Morphological changes of the workers developing into neotenic reproductives
The differentiation of neotenic reproductives (NRs) was induced by isolating late instar workers of R. labralis (Fig.1). The difference in morphology between the female workers and NRs was significant (Fig.2a-d).The abdomen lengths of the workers, isolated workers and NRs were 1.89±0.06,1.95±0.03 and 2.70±0.11mm, respectively. The abdomen length significantly increased after the female workers moulted into NRs. The heads of the female NRs had brown pigmentation stripes. No male NRs emerged in experimental colonies. Moreover, the ovaries of the workers produced degenerate follicles. In the ovaries of the workers, a conspicuous feature was that each oocyte was surrounded by a thin layer of follicle cells in which the elongated follicle cells were degenerated (Fig.2f). In NRs, each oocyte was surrounded by a thicker layer of follicle cells (Fig.2e and g).
RNA sequencing and de novo transcriptome assembly
A total of 460,864,858 clean-sequence reads were obtained by transcriptome sequencing (Additional files 1). Based on the clean reads, 112,954 unigenes ranging from 201 bp to 38,952 bp were assembled. The mean length was 848 bp, and the N50 length was 1,536 bp (Additional files 2 and 3).
We annotated 112,954 unigenes. The Venn diagram illustrates that the number of unique sequence-based annotations is the sum of the best unique BLASTx hits from the Nr, Swiss-Prot, KEGG and KOG databases (Fig.3). A total of 40,972 unigenes (36.27%) were successfully annotated to at least one database. 17,535 unigenes (15.52％) were successfully annotated in the four databases. In total, 40,073 (35.48％) unigenes had significant matches in the Nr database, followed by 29,540 unigenes (26.15％) in the Swiss-Prot database. 25,453 (27.99％) and 20,116 unigenes (22.53％) had specific matches in the KOG and KEGG databases, respectively (Fig. 2a and Additional files 4). As expected, the top hits were found in the insect genomes, especially termite Zootermopsis nevadensis (25.82％) (Additional files 5 ).
Based on the Nr protein database annotation, we classified the function of the predicted unigenes using GO, KOG and KEGG (Fig.3). In the GO functional classification, the annotated unigenes belonged to clusters of biological process (30,717 unigenes, 36.68%), molecular function (29,711 unigenes, 35.48%) and cellular components (23,312 unigenes, 27.84%), which were distributed into 55 categories (Additional files 6). For the KOG database, 25,453 unigenes were classified into 25 functional categories (Fig. 3). The largest group was general function prediction(10,042 unigenes, 39.45%), followed by signal transduction mechanisms (10,007 unigenes, 39.31%), posttranslational modification, protein turnover, and chaperones (8,124 unigenes, 31.92%). A total of 20,116 unigenes were annotated into 233 pathways in the KEGG database. The most enriched pathway was “Ribosome” (1,533 unigenes, 6.26%), followed by “Protein processing in the endoplasmic reticulum” and “Endocytosis”, for which the unigene percentages were 4.48% and 3.36%, respectively (Additional files 7).
Differentially expressed genes (DEGs) in workers, isolated workers and neotenic reproductives
We performed a quantitative comparison of the gene expression between workers and IWs (worker vs IW), between IWs and NRs (IW vs NR), and between workers and NRs (worker vs NR) (Fig. 4). We identified 38,070 DEGs (upregulated and downregulated genes) among workers, IWs and NRs. There were 17,405 DEGs in “workers vs IWs”, 30,332 DEGs in “IWs vs NRs” and 7,016 DEGs in “workers vs NRs”. The results showed that “worker vs NR” had the fewest number of DEGs compared with “worker vs IW” and “IW vs NR”. The number of DEGs from “IW vs NR” and “worker vs IW” was 2.5-fold and 4.3-fold, respectively, compared with the number of DEGs of “worker vs NR”. Of the 17,405 DEGs, 16,910 (97.2%) genes were upregulated in IWs compared with workers. After the IWs developed into NRs, most of the DEGs (27,741, 91.5%) were downregulated.
Trend analysis of DEGs enrichment across the developmental stages of workers, IW into NRs
The 38,070 DEGs from workers, IWs and NRs of R. labralis were clustered into 8 profiles (Fig. 5a), in which 32,622 (85.69%) genes were in profile 3 and profile 5 (Fig. 5b). Only in profile 5, all 12,543 genes were specifically upregulated in IWs and were downregulated in workers and in NRs. The expression levels of the genes in profile 5 were significantly higher in IWs than those in the workers and NRs, which indicated that these genes were specifically overexpressed during the stage of IWs developing into NRs, and then were rapidly downregulated after the IWs became NRs. The genes in profile 5 were involved in the differentiation of IWs into NRs. Therefore, we focused on the gene expression in profile 5.
To determine the biological function of the DEGs in profile 5, the GO classification and the KEGG pathway analysis were carried out. All of the DEGs were annotated in the GO database and were categorized into 51 functional groups, including the three main GO ontologies: biological process, cellular component and molecular function (Additional files 8). Of these DEGs, a large number were dominant in catalytic activity, binding, metabolic processes, cellular processes, the cell part and the cell.
Our KEGG pathways analysis showed that 2,272 DEGs in profile 5 were mapped in 201 pathways, in which 31pathways were significantly related to IWs differentiation into NRs (Q-value <0.05) , including ribosomes, protein processing, biosynthesis of amino acids, calcium signalling pathway, etc.(Additional files 9). Five signal transduction pathways were involved in profile 5, and the signal transduction pathways in the top 20 pathways were the phosphatidylinositol signalling pathway, the calcium signalling pathway, the Ras signalling pathway and the PI3K-Akt signalling pathway (Fig.5c). The number of genes enriched in the above pathways were 111 (4.89%), 67 (2.95%), 41 (1.8%) and 27 (1.19%) unigenes, respectively. All of the five signal transduction pathways in profile 5 were classified into environmental information processing, which indicated that these signal transduction pathways were related to the transition of workers into NRs after the workers were isolated from their natal colonies.
Active transposable elements (TEs) can “jump” within the genome, thereby disturbing the regulation and expression of other genes by transposing into another gene. In profile 5, 2 of 739 annotated TE genes were upregulated in IWs, one PogoTE gene and one PiggyBac TE gene, and then were downregulated in NRs.
Signal transduction along the Ras-MAPK pathway axis from plasma membrane to nucleus Ras proteins transmit signals from cell surface receptors to a variety of effectors and thereby regulate pathways that govern cell proliferation and differentiation. Ras proteins are binary switches that cycle between ON and OFF states during signal transduction. In profile 5, we found that 41 genes involved in the Ras signalling pathways were significantly upregulated in IWs compared with workers, and then were significantly downregulated in NRs (Fig.5d). All expression levels of genes in the IWs were significantly higher than those in the workers (natal colonies) and in the NRs.
The transcriptome analysis showed that a Ras (N-Ras protein) gene was expressed only in IWs and was not detected in workers and NRs, and Ras-extracellular signal regulated kinase (ERK) was significantly upregulated in IWs. The signalling pathways in profile 5 showed that Ca2+- Calmodulin (Ca2+-CaM) involved the activation of the Ras signalling pathway, and then activated the downstream effector pathway MAPK during IW development into NRs (Additional files 10).We identified that 52 genes in the MAPK signalling pathways was significantly upregulated in IWs and was significantly downregulated in workers and in NRs. 61 genes involved in the calcium signalling pathway were significantly upregulated in IWs and were significantly downregulated in workers and in NRs. Predicted model for Ras-MAPK pathway controlling the differentiation of workers into neotenic reproductives (queens) was shown (Fig. 6).
Dynamic changes in Ras signalling pathway-related gene expression during the differentiation of workers into NRs
We performed qRT-PCR analysis on the expression of five Ras signalling pathway-related genes. The expression levels of the five genes in the isolated female workers (IWs) were higher compared with the expression levels of the workers and neotenic reproductives(NRs). The relative expression level of Ras in the IWs was 131-fold and 20-fold higher than the expression in the workers and NRs, respectively (Fig. 7A). The relative expression level of Calmodulin (CaM) in the IWs was 333-fold and 47-fold higher than the levels in the workers and NRs, respectively (Fig.7B). The relative expression level of the calcium-binding protein in the IWs was 462-fold and 205-fold higher than the levels in the workers and NRs, respectively (Fig.7C). The relative expression level of the Ras GTPase-activating protein 4 (Ras-GDP) in the IWs was 48-fold and 7-fold higher than the levels in the workers and NRs, respectively (Fig. 7D). The relative expression level of protein kinase C (PKC) was 18-fold and 11-fold higher in the IWs than the levels in the workers and NRs, respectively (p<0.05) (Fig. 7E). The five genes in the IWs were extremely upregulated. After the IWs differentiated into NRs, the expression levels of the five genes in the NRs derived from IWs extremely decreased. There were no significant differences in the expression levels of the five genes between the workers and the NRs.
The data from transcriptome sequencing also showed that the five genes were significantly upregulated in IWs compared with workers and NRs, and there were no significant differences in the expression levels of the five genes between the workers and the NRs. The expression of Ras, Calmodulin and calcium-binding genes in the workers and NRs was not detected. Our results showed that the expression profiles of the candidate unigenes revealed by the qRT-PCR data were consistent with those derived from transcriptome sequencing (Fig. 7A-E), which indicate that the RNA-Seq analysis was reliable and provided a valuable gene sequence for biological analysis. Ras signalling pathway-related genes exhibited a significant IW-specific overexpression and have putative associations with the reproductive plasticity of workers.
Dynamic changes in the expression levels of catalase gene related to longevity during the differentiation of workers into NRs.
The oxidative stress theory of ageing states that the accumulation of oxidative damage causes ageing. In the model insect Drosophila melanogaster, overexpression of catalase (CAT) resulted in reduced levels of oxidative stress and an extended lifespan, which indicated that CAT positively effects longevity[18,19]. In this study, the CAT of R. labralis closely matched the CAT gene sequence of the termite Zootermopsis nevadensis. We found a dynamic change in the expression of the CAT gene during the differentiation of workers into reproductives. Our qRT-PCR analysis showed that the expression levels of the CAT gene decrease significantly in isolated female workers for one week (IW1), isolated female workers for two weeks (IW2) and isolated female workers for three weeks (IW) compared with those of workers. However, after the IW moulted into NR, the expression levels of the CAT gene increased significantly, and the expression levels of the CAT gene in NRs (1.670 ± 0.065) was approximately 5-fold higher compared with those of IW(0.316 ± 0.012) (Fig. 8). Our KEGG pathway analysis showed that the CAT gene, as a downstream gene in the longevity-regulating pathway, directly caused longevity. However, the expression of CAT was usually suppressed by the Ras-P13k-Akt-FOXO pathway (Additional files 11).