Transcriptomic evidence that insulin signalling pathway regulates the ageing of subterranean termite castes

Insulin is a protein hormone that controls the metabolism of sugar, fat and protein via signal transduction in cells, influencing growth and developmental processes such as reproduction and ageing. From nematodes to fruit flies, rodents and other animals, glucose signalling mechanisms are highly conserved. Reproductive termites (queens and kings) exhibit an extraordinarily long lifespan relative to non-reproductive individuals such as workers, despite being generated from the same genome, thus providing a unique model for the investigation of longevity. The key reason for this molecular mechanism, however, remains unclear. To clarify the molecular mechanism underlying this phenomenon, we sequenced the transcriptomes of the primary kings (PKs), primary queens (PQs), male (WMs) and female (WFs) workers of the lower subterranean termite Reticulitermes chinensis. We performed RNA sequencing and identified 33 insulin signalling pathway-related genes in R. chinensis. RT-qPCR analyses revealed that EIF4E and RPS6 genes were highly expressed in WMs and WFs workers, while mTOR expression was lower in PKs and PQs than in WMs and WFs. PQs and PKs exhibited lower expression of akt2-a than female workers. As the highly conserved insulin signalling pathway can significantly prolong the healthspan and lifespan, so we infer that the insulin signalling pathway regulates ageing in the subterranean termite R. chinensis. Further studies are recommended to reveal the biological function of insulin signalling pathway-related genes in the survival of termites to provide new insights into biomolecular homeostasis maintenance and its relationship to remarkable longevity.

www.nature.com/scientificreports www.nature.com/scientificreports/ Growth factors such as growth factor-1, insulin and rapamycin have been extensively investigated, and sirtuin genes have been identified as genes that are upregulated in association with longevity in animals whose lifespan is extended due to caloric restriction 10 . Therefore, animals (e.g., D. melanogaster, C. elegans, D. longispina and A. islandica) that produce insulin-like peptides, mutations that reduce the levels or actions of these peptides may lead to significant increase in lifespan. The genetic control of aging and longevity of C. elegans results in more than 30 insulin-like peptides signalling the species via a common Daf2 receptor. It has been shown that the Daf2 gene exhibits significant homology to mammalian genes for the insulin receptor and the insulin-like growth factor 1 (IGF-1) receptor 11,12 . Consequently, insulin/IGF-1-like pathways also mediate metabolic effects relying on the winged-helix transcription factor, FOXO. In C. elegans FOXO is first identified as Daf-16, which elevates lipid level and longevity induced by Daf-2 loss 13 . Similarly, D. melanogaster has a homologous insulin/IGF-1-like pathway consists of insulin receptors, including mammalian insulin receptor substrates (IRS) and kinase homologs, and a forkhead transcription factor. These insulin receptors reduce signalling stress resistance and extend the longevity 14 . Senescence represents a decline fitness mechanism in age and life span, a highly plastic life-history trait that can be strongly influenced by the biotic and abiotic characteristics of the ambient environment 15 . As a result, the evolutionary and ecological importance of such declines in free-living populations and the intensively studied model laboratory systems is mainly unknown among individuals surviving into old age 16 . Molecular genetic studies have shown that insulin transduction regulates many physiological phenomena, such as the growth, development and longevity of insects 17 . The mutant insulin receptor substrate gene chico can significantly prolong the lifespan of female fruit flies 18 . In addition, the PI3K-Akt is the main pathway of the two-central insulin signalling pathways (ERK/MAPK and PI3K-Akt) in insects 19 , we, therefore, follow PI3K-Akt pathway to determine the insulin signalling pathway in the social insect of termite on the first time.
Previous studies of ageing in termites have focused mainly on antioxidant enzymes, DNA damage, and transposable element activity [20][21][22] . For example, the expression levels of DNA repair genes in the non-reproductive and primary reproductive states of Reticulitermes speratus were determined 19 , and several antioxidant enzyme genes were reported as potentially being associated with fertility and longevity in different castes 20 . Similar studies relating termite insulin signalling pathway-related genes and longevity have not yet been reported. The present study compared the expression levels of insulin signalling pathway-related genes among different castes of the subterranean termite Reticulitermes chinensis, including WMs, WFs, PQs and PKs, to reveal the mechanisms underlying different lifespans in these termite castes.

Results
illumina data sequence and assembly. We established RNA-seq libraries using mRNA isolated from WMs, WFs, PQs and PKs of R. chinensis. A total of 109,126,456 clean sequencing readings were obtained via Hi-SeqTM 4000 (Illumina) paired-end sequencing. Each test provided more than 7 G of transcriptome data based on the clean reads, and the Trinity system assembled a total of 161,932 unigenes ranging from 201 bp to 19,428 bp ( Supplementary Fig. 1). The mean length was 673 bp, with an N50 value of 933 bp. The performance evaluation based on GC showed an accurate approximation (42.90%). These findings suggest that the development and quality of the sequencing results were sufficient for further research. functional annotation of the R. chinensis transcriptome. For unigene annotation, we applied BLASTx (http://www.ncbi.nlm.nih.gov/BLAST/) with an e-value cut-off of 1e-5 to the NCBI nonredundant (Nr) protein database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (http://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg) [23][24][25] , and the KOG database (ftp://ftp.ncbi.nih.gov/pub/COG/KOG/kyva). In total, 161,933 unigenes were annotated. The Venn diagram showed that the number of unique sequence-based annotations was the sum of the unique best BLASTX similarities from the Nr, Swiss-Prot, KOG and KEGG databases (Fig. 1). In the Nr database, 60,736 unigenes (37,51% of R. chinensis unigenes) presented significant matches, whereas in the Swiss-Prot database, 31,699 unigenes (19,58%), in the KOG database, 27,181 unigenes (16,79%), and in the KEGG database, 28,003 unigenes (17,30%) exhibited significant matches (Supplementary Table 1). Most of the sequences (100,843) were compiled, and 62.27% of the genes could not be identified due to the relatively short length of the distinct gene sequences and the lack of annotation of the R. chinensis genome information. These unigenes may also be R. chinensis-specific genes or short fragments that mainly come from untranslated regions (e.g., 5′ and 3′ UTRs) or nonconserved regions of protein-coding transcripts (Fig. 2) 26 . The experiment was designed for the whole termite, and the Nr full library was used for unigenes, where data on the intestinal protist could not be ruled out because it was impossible to determine whether the unigene annotated as a protist species was a termite or an intestinal protist gene, so it was kept during analysis.

Gene Ontology (GO), Clusters of euKaryotic Orthologous Groups (KOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) ontology classifications.
We classified the functions of the predicted unigenes using GO, KEGG, and KOG analyses based on the protein annotation results of Nr database homology searches. A total of 52,430 unigenes were annotated into 25 classes (Fig. 3) in the KOG functional classification.
To determine the biological functions of the DEGs among PKs, PQs, WMs and WFs, GO classification was carried out. A total of 59 categories (Fig. 4) of functional groups were analysed. Biological processes showed the highest probability density distribution of gene expression, followed by cellular components and molecular functions. All DEGs identified in this study were mapped to terms in the GO database to determine the functions of the differentially expressed genes. These included cellular processes, catalytic activity, metabolic processes, cell components, single-organism processes, and binding ( Fig. 5A-D). Through sequencing, we found a total of 568 related genes involved in the insulin signalling pathway (Supplementary Fig. 4). There are 33 genes exhibiting significant differences in the primary kings, primary queens, and female and male workers of the lower www.nature.com/scientificreports www.nature.com/scientificreports/ subterranean termites. There are two insulin signalling pathways: the mitogen-active kinase (MAPK) pathway and the PI3K-Akt pathway. However, the PI3K pathway is the most prominent insulin signalling pathway in insects. Therefore, we selected 6 genes (Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6) in the PI3K-Akt pathway showing significant differences.

Differentially expressed gene (DEG) analysis in
The RPKM and RT-qPCR validation was performed for 6 genes found to be differentially expressed in the scatter plot correlation among the WMs, WFs, PQs and PKs castes. The correlation coefficient between the RT-qPCRs and transcriptome data validation results was 0.721 (Fig. 10).

Discussion
Age verification or determination is a requirement for the identification of individuals in living organisms and for the prediction of death. Drosophila melanogaster and Caenorhabditis elegans are well-established research organisms in the field of ageing research. These organisms may not only exhibit a prolonged life but may also be healthier throughout their life span. Studies among social insects, including termites, show that there is diverse caste differentiation (workers, soldiers, king and queen) within a single colony. The queen and king can live for 20-30 years, while non-reproductive castes only live for a few weeks to months. The specific associated developmental, biological and physiological changes in a given cell or tissue type have been examined through the evaluation of transcriptome sequences 19,27 . Transcriptomics refers to the collection of all RNA molecules, from protein-coding (mRNA) to non-coding RNAs, including rRNAs, tRNAs, lncRNAs, and pri-miRNAs among others, for which specialized library preparation methods and appropriate bioinformatic data processing and abundance quantification techniques are required for functional analysis 26 .
To understand global changes in the expression of thousands of genes, transcriptome analysis requires an effective statistical method with multiple comparison tests. The accuracy of quantitative RT-qPCR and microarray analyses has become evident in recent years and is highly dependent on the selection of genes for standardization. However, the present study is the first on termite insulin signalling pathway-related genes. From the transcriptome data, we obtained a total of 33 DEGs involved in the insulin signalling pathway by producing a high-quality dataset from R. chinensis (Fig. 6). These genes (Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6) are involved in the insulin signalling pathway in different insects (i.e., Drosophila) 28 . Rintelen 29 investigated whether in vivo Pdk1 presents more than one target and functions in PI(3)K signalling at the downstream level in C. elegans. In addition, many in vitro invertebrate cell culture studies have indicated that Pdk1, in the AGC family of kinases (consisting of Akt, S6K, RSK, PKN and all protein kinase C isoforms), is the critical regulator of T-loop phosphorylation. Pdk1 has two active domains, an amino-terminally bound serine-threonine kinase domain and a pleckstrin homology domain with excessive PIP3 affinity. Pdk1 is the direct effector of Akt, S6K, and RSK, and the activation of all three of these kinases is blocked in Pdk1-deficient embryonic stem cells. However, previous genetic evaluations of C. elegans and Drosophila were primarily based on the finding that Akt plays a significant role in Pdk1 activity as a central regulator of cell growth through AGC kinases Akt and S6K, respectively 28,29 . Recently, Gao et al. 30 examined the role of PI3K-Akt and FOXO proteins in insulin signalling pathways in Sogatella furcifera. They reported physiological and genetic changes in the wing patterning of S. furcifera.
A similar study was reported in C. elegans indicating that Pdk1 regulates the physiological effects of insulin and growth factors by activating a series of kinases controlling cell growth, differentiation, survival, protein www.nature.com/scientificreports www.nature.com/scientificreports/ translation and glucose metabolism. If inactivation of Pdk1 occurs, C. elegans will enter the stagnant dauer stage and extend its lifespan 31 . This phenomenon (inactivation of Pdk1) was also reported in primary reproductives (king and queen). The results of the present study regarding the expression levels of Pdk1 in different termite castes showed that there was no statistically significant difference between reproductive and worker termites. Consequently, our results indicated that the function of Pdk1 in regulating termite lifespan showed the same changes as in nematode and other insects 28,30 . Therefore, the low expression 31 (C. elegans) of Pdk1 had little effect on the longevity or physical morphology of the termites. The Pdk1 induced bent wing phenotype in S. furcifera depends on average levels of dS6K and Akt, because null mutations in either of the corresponding genes dominantly suppress the longevity. These findings, together with biochemical evidence from cultured cells showing that Pdk1 regulates the activity of Akt and dS6K, provide functional evidence that Pdk1 is a key regulator of growth and cell size by controlling the activity of two AGC kinases, Akt and dS6K 29 . The akt2-a belongs to a subfamily of serine/threonine protein kinases consisting of three members, akt1, akt2 and akt3. The activation of these genes involves a combination of numerous stimuli and the activation of hormones, metabolism, growth factors and cell motility 32,33 . In the classical control model of Akt1, the PIP3 phospholipid recruits akt1 to the www.nature.com/scientificreports www.nature.com/scientificreports/ plasma membrane, where two protein kinases, mTOR and Pdk1, phosphorylate Akt1 at its C-terminus and activation loop, respectively. The dual phosphorylation of Akt1 results in increased kinase activity targeting protein substrates such as GSK3 and FOXO. The activity of Akt1 may promote cell growth, block apoptosis, and mediate the insulin response, and clinically produced akt inhibitors include binding compounds targeting the ATP site and allosteric site 34,35 . Our findings suggest that the low expression of akt2-a in the primary reproductive castes may also prevent tumour cell invasion and metastasis. A related cell mutation study was conducted for the Tsc1 Drosophila homologue via mosaic screens, and the size of Tsc1 mutant cells was significantly increased. The body volume was also increased, specifically by the tissues that contained the most mutant cells. Mutations in the Drosophila, Tsc2 gene have previously been shown to induce similar phenotypes, and it has been proposed that polyploidy triggers the switch in the cell. Clones of Tsc1 mutant cells undergo certain divisions within imaginal disks; however, they retain normal ploidy. The ectopic overexpression of Tsc1 or Tsc2 in Drosophila tissue does not inhibit but competes with the expression of Tsc1 and Tsc2, leading to a great decline in cell growth and proliferation 36 . Tsc1 and Tsc2 are genes responsible for the suppression of tumours that contain a protein called tuberin 37 . Tsc2 is widely distributed across cell types and organ systems, and these genes and proteins are highly conserved in interspecies sequences from Drosophila to humans. Tsc2 binds to a third protein, TBC1D7, as part of a heteromeric protein complex to control cell growth, cell size, the cell cycle, and mTOR pathway proliferation. In contrast, large deletions, indels, nonsense and missense mutations and splicing errors are included in Tsc2 mutations 38 . The hypothesized findings (RT-qPCR results) were obtained in the present study from the reproductive castes, in which higher expression of Tsc2 was identified than in the non-reproductive castes of this termite.
Ageing is described as an accumulation of cellular damage over time promoting disease and death 39 . The germline genome is vulnerable to the collection of deleterious mutations throughout meiotic DNA replication. If a mutation is not eliminated from reproductive cells, it can be passed on to offspring, which is linked to an elevated threat of diseases in future generations. When mutations occur in somatic cells, they cannot be transmitted to the next generation. However, if mutations occur in the gametes of an organism, they may affect the offspring 40 . The mammalian target of rapamycin (mTOR) has been reported in many organisms, such as yeast and mammals. mTOR belongs to the serine/threonine kinase family of phosphatidylinositol kinase-related kinases (PIKKs) 41 and contributes to increases in growth factors, physiological processes, cell metabolism and survival, and autophagy. The basal amount of mTOR leads to a low fecundity rate, egg size, and follicle numbers and is correlated with low vestigial (Vg) expression in Aedes aegypti, Nilaparvata lugens, Drosophila and Apis mellifera 42,43 . Accumulating evidence suggests that aberrant regulation of both cell growth and metabolism substantially contributes to cancer improvement and progression 44 . Therefore, earlier research has indicated that the high expression of akt2-a, mTOR, EIF4E, and RPS6 is closely associated with cancer, growth factors, reproduction, physiological processes (wing patterning of S. furcifera) 30 , the inflammatory response, cell survival, weight problems and autophagy [45][46][47][48] . Furthermore, previous researchers discovered that the high expression of the breast cancer susceptibility gene BRCA1 leads to a long life in insects 49,50 . Hence, an efficient antioxidant system can prevent the accumulation of detrimental DNA changes and contribute to the longevity of termite kings. Long-lived reproductive individuals exhibit a stable defence system against transposons as a possible source of DNA damage, in contrast to short-lived workers of the termite Macrotermes bellicosus 51 . The present experimental results indicated that the expression of cancer-related genes is low in the primary reproductive castes of termites. They are unlikely to experience DNA damage because cancer incidence and longevity are associated with multi-gene regulation processes [52][53][54][55][56][57][58][59] . Consistent with these studies, our results suggest that extreme evolutionary pressures (body changes over time) 4 potentially led to the low expression of insulin signalling pathway-related genes in primary reproductive castes of R. chinensis.
In D. melanogaster and C. elegans, the most prominent role of the insulin signalling pathway is to control longevity, puberty, growth and body size. The ability to generate observable aberrant phenotypes through the www.nature.com/scientificreports www.nature.com/scientificreports/ perturbation of cell growth has allowed rapid progress in growth regulation studies. These species show strong similarities in their food reactions and sensory compensation associated with insulin signalling pathways and insulin-like peptides. The insulin signalling pathway also controls the metabolic pathway in the fly. Interestingly, flies exhibit female sterility, increased triglyceride levels and an extended diapause life stage as physiological responses to harsh environmental conditions such as a low nutrient supply or low temperature. Increased stress tolerance often occurs during diapause, which together with greater energy reserves, increases longevity and therefore the probability of reproduction [60][61][62] . Similarly, the results provide a valuable resource for the study of ageing mechanisms, structures and related pathways. These findings suggest that relatively conserved proteins alter the insulin signalling pathway and considerably prolong the lifespan or even prevent the development of diseases related to ageing. Further studies are needed to reveal the biological function of insulin signalling pathway-related genes in the survival of termites using genetic tools such as RNA interference and transgenic constructs [63][64][65] . Ultimately, our results provide new insights into biomolecular homeostasis maintenance and its relationship to remarkable longevity.

Materials and Methods
Sample collection. Alates of R. chinensis were collected from colonies in Chengdu in April 2014. The initial founder colony (each colony consisted of one male and one female alate) was by rearing a randomly selected male and female alate in a plastic box (80× 65× 40 mm) with damp chips of sliced pine wood at 25 °C in the laboratory (Termite House of Northwest University, China) (Fig. 11D). Individuals were extracted from each early-stage colony after 4 years, and they were regarded as PKs, PQs, and WMs and WFs workers (Fig. 11). As biological replicates for the experimental procedure, the whole bodies of the termites were temporarily preserved at −80 °C in liquid nitrogen. The WMs, WFs, PQs and PKs were separated (as a male and female) based on the characterization of the seventh sterna under an anatomical microscope. Reticulitermes chinensis is not an endangered or protected species; thus, no specific permits were required for the sampling and experiments.

RnA extraction and illumina sequencing.
To obtain enough RNA, from PQs, PKs, WMs and WFs samples using TRIzol reagent and checked the total RNA quality assessed an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Three replicates of each caste WMs, WFs, PKs and PQs, were pooled to obtain enough RNA and build cDNA libraries using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB). In brief, using Oligo(dT) beads to extract the total RNA, enriched from mRNA and fragmentation buffer were used to transcribed into short fragments into cDNA with random primers. Second-strand cDNA was synthesized by using DNA polymerase I, RNase H, dNTPs and buffer. The cDNA fragments were then purified with the QiaQuick PCR extraction kit, poly(A) end-repaired was tailed, and attached to the Illumina sequencing adapters. Ligation products were selected in size by agarose gel electrophoresis, PCR amplified and sequenced using the Illumina HiSeqTM 4000 platform by Gene Denovo Biotechnology Co. (Guangzhou, China) 66,67 . De novo transcriptome assembly. The remaining reads from all samples were assembled using the Trinity version of trinityrnaseq r2012-04-27 68,69 , which generates transcriptomic assemblies from short-read sequences using the de Bruijn graph algorithm. Further description of the methodology was previously provided in detail 48 , as was a summary of the assembly statistics. Assemblies were generated for WMs, WFs, PKs and PQs (each biological replicate consisted of 5 individuals) to obtain appropriate RNA for RNA-seq.
Read alignments and normalization of gene expression levels. SOAPaligner/soap2 70 was used for short sequence alignment to read the sequences and map them to the reference sequences. The read coverage of one gene was used to calculate the level of expression of that gene. We obtained the expression levels of all the detected genes using this process. A read that was uniquely mapped to a gene was used to calculate the degree of expression. The level of gene expression was measured from the number of reads per kilobase of the exon region per million mapped reads (RPKM). The R package (http://www.r-project.org/) was used for statistical data expression and visualization.
Differentially expressed genes (DEGs) and functional enrichment analyses. After calculating the rate of the expression of each gene, differential expression analysis was performed using edgeR 71 . After multiple tests, the false discovery rate (FDR) was used to determine the threshold for the P-value, and an FDR threshold of 0.01 and an absolute value of log2Ratio ≥ 1 were used to measure the significance of differences in gene expression in the sample. Using a method similar to that described by Zhang 72 , the differentially expressed genes were subjected to the analysis of GO and KEGG enrichment. Operational annotations for unigenes can be extracted from the analysis of Nr annotations. The GO annotation profile was checked with Blast2GO software 73 . The functional classification of unigenes was performed using WEGO software 74 . For DEGs, a Q-value <0.05 for GO terms and KEGG pathways were considered to indicate significant enrichment. Quantitative real-time PCR (RT-qPCR). We designed primer pairs for each gene associated with the insulin signalling pathway with Primer3 v1.1.4 (Supplementary Table 2). Total RNA was extracted from the entire body of individuals of each caste (PQs, PKs, and WFs and WMs workers) using the RNAsimple Total RNA Kit (Tiangen). After the extraction of RNA, its quality and quantity (purity of protein and salt) were checked with a NanoReady spectrophotometer (Model: F-1100 made in China). To build the cDNA (the cDNA was held at -20 °C for further experiments) library, the FastKing RT Kit (Tiangen) was used. RT-qPCR was used to amplify the cDNA by using a CFX 96 instrument (Bio-Rad) with SuperReal PreMix Plus (Tiangen). All procedures were carried out in compliance with the manufacturer's protocol. Since beta-actin (RsACT) 75,76 was evaluated as the most reliable reference gene for Reticulitermes termites, it was selected as the reference 77 . The standard 2 −∆∆Ct method was used to calculate relative gene expression 78 . We conducted three biological and three experimental replicates (each replicate consisted of 5 individuals) for all RT-qPCR experiments.
www.nature.com/scientificreports www.nature.com/scientificreports/ Statistical analysis. For the statistical analysis of the RT-qPCR experiments, SPSS Statistics 17.0 was used.
Significance between groups was measured through a one-way variance analysis (ANOVA) followed by the Duncan post hoc test. All data in graphs are shown as the mean ± standard mean error, and all estimated P values <0.05 are given.

Data availability
The datasets generated and analyzed during the current study are deposited under BioProject accession number PRJNA592596 at the NCBI. Any reasonable requests will be answered by the corresponding author.