Morphological changes of the workers developing into NRs
The differentiation of NRs was induced by isolated 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 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[12]. 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).
Fig. 1 The workers of Reticulitermes labralis created new nests (colonies) in the field and laboratory. (a) In the field, the workers created a new colony in new food site (dead wood). The workers exploited new food
resources outside the natal nest and moved to a new nest site. The differentiation
of female workers into NRs was triggered in the new nest site. (b) In the laboratory, A group of workers from their natal colony was isolated in a Petri dish with moist
sawdust. After four weeks, the NRs appeared in the group of isolated workers, indicating
that a new colony was created. W(blue), workers; PR(red), primary reproductives; NR(yellow),
female neotenic reproductives.
Fig.2 Morphological changes of the workers developing into neotenic reproductives
(NRs). (a) The primary reproductives with darkened pigmentation. Scale bar =1.0 mm. (b) NRs (queens) in colony. Scale bar =1.0 mm. (c) A female workers. Scale bar =1.0 mm.
(d) The NRs derived from female workers had brown pigmentation on their heads. Scale
bar =1.0mm. (e) The ovary of NRs. Scale bar = 100 µm. (f) Each oocyte in workers was
surrounded by degenerated follicle cells. Scale bar = 20 µm. (g) Each oocyte in the
NRs was surrounded by a thick layer of cube-shaped follicle cells. Scale bar = 25
µm.
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).
A total of 112,954 unigenes were used for functional annotation. 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. 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 first nine pathways were
shown in Additional files 7. 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.
Fig. 3 A total of 112,954 unigenes were used for functional annotation. (a) Venn diagram of the distribution of unigene and database matching results. The
numbers of unique sequence-based annotations is the sum of the unique best BLASTX
hits from the Nr, Swiss-Prot, KOG and KEGG databases. The overlapping regions between
the four circles contain the numbers of unigenes that shared BLASTX similarities with
the respectively databases. (b) Histogram presentation of the clusters of the KOG
function classifications. A total of 25,453 unigenes were grouped into 24 KOG classifications.
The y-axis indicates the number of unigenes in a specific functional cluster. The
legend presents the 24 functional categories.
Differentially expressed genes (DEGs) in workers, IWs and NRs
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 least 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.
Fig.4 Differentially expressed genes in workers, IWs and NRs. (a) The number of the
differentially expressed genes in workers, IWs and NRs. The x-axis indicates the three
stages (worker: female workers in natal colonies; IW: isolated female workers from
their natal colonies; NR: female neotenic reproductives differentiated from isolated
workers). Red represent transcripts that were significant up-regulated, and green
indicate that those transcripts were significantly downregulated. The parameters FDR≤0.001
and log2Ratio≥1 were used as the thresholds to judge the significance of gene expression
differences. (b) The scatter plot of the differentially expressed genes in workers,
IWs and NRs. Red and green scatter plot represent upregulated genes and downregulated,
respectively. Black scatter plot represent that the genes were not differential expressed.
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.
Activity transposable elements (TEs) can “jump” within the genome, thereby disturbing
the regulation and expression of other genes by transposing into another gene. In
IWs, 2 of 739 annotated TE genes were upregulated, one PogoTE gene and one PiggyBac
TE gene, and then were downregulated in NRs.
Fig. 5 Trend analysis of DEGs enrichment across the developmental stages of workers,
IW into NRs. (a)Trend analysis of differentially expressed genes. The 38,070 DEGs
from workers, IWs and NRs of R. labralis were clustered into 8 profiles. (b) Significantly enriched trend analysis of profile
3 and profile 5.(c)Thop 20 of KEGG pathway enrichment of DEGs in Profile 5.The bubble
size represents the number of DEGs, and the bubble color represents the Q value.(d)Heat
map of the DEGs for the 41genes involved in Ras signaling pathway.
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
transduction15. 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).
Fig. 6 Predicted model for Ras-MAPK pathway controlling the differentiation of workers
into neotenic reproductives. When cell surface receptor accepts extracellular signaling
from the absence of queen pheromone, a increase in intracellular Ca2+ is triggered via Ca2+ channels. Ca2+ enters the cytoplasm where it binds and activates calmodulin(CaM). Ca2+-CaM can act on Ras in the cytoplasmic surface. Ras binds GTP. Activated Ras initiates
the Raf-1, and through the MEK/ERK MAPK pathway, triggers signaling to the nucleus.
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 IWs were higher compared with the expression
levels of the workers and 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.
Fig. 7 The Comparison of expression of five genes between RNA-Seq and qRT-PCR during the differentiation of workers into
NRs. (A) Ras protein; (B) Calmodulin ; (C) calcium-binding protein; (D) Ras GTPase-activating
protein 4; (E) protein kinase C. Worker, the female workers in natal colonies; IWs,
isolated female workers for three weeks; NRs, female neotenic reproductives. The blue
color represents the qRT-PCR analyses of the gene expression levels. The red color
represents the transcriptome sequencing analyses (RPKM) of the gene expression levels
.
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[16,17].
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 IW1, IW2 and IWs compared with those of workers. However, after the IWs moulted into NR, the
expression levels of the CAT gene increased significantly, and the expression levels
of the CAT gene in NR (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
in an uncomfortable environment (Additional files 11).
Fig. 8 Dynamic changes in the expression levels of catalase gene related to longevity
during the differentiation of workers into NRs by qRT-PCR analysis. Error bars represent
the standard deviation of the mean. Different letters over the bars denote significant
differences atp< 0.05. Worker, the female workers in their natal colonies; IW1, isolated female workers for one week; IW2, isolated female workers for two weeks; IW, isolated female workers for three weeks;
NRs, female neotenic reproductives.