De Novo Transcriptome Assembly and Analysis of Longevity Genes Using Subterranean Termite (Reticulitermes chinensis) Castes

The longevity phenomenon is entirely controlled by the insulin signaling pathway (IIS-pathway). Both vertebrates and invertebrates have IIS-pathways that are comparable to one another, though no one has previously described de novo transcriptome assembly of IIS-pathway-associated genes in termites. In this research, we analyzed the transcriptomes of both reproductive (primary kings “PK” and queens “PQ”, secondary worker reproductive kings “SWRK” and queens “SWRQ”) and non-reproductive (male “WM” and female “WF” workers) castes of the subterranean termite Reticulitermes chinensis. The goal was to identify the genes responsible for longevity in the reproductive and non-reproductive castes. Through transcriptome analysis, we annotated 103,589,264 sequence reads and 184,436 (7G) unigenes were assembled, GC performance was measured at 43.02%, and 64,046 sequences were reported as CDs sequences. Of which 35 IIS-pathway-associated genes were identified, among 35 genes, we focused on the phosphoinositide-dependent kinase-1 (Pdk1), protein kinase B2 (akt2-a), tuberous sclerosis-2 (Tsc2), mammalian target of rapamycin (mTOR), eukaryotic translation initiation factor 4E (EIF4E) and ribosomal protein S6 (RPS6) genes. Previously these genes (Pdk1, akt2-a, mTOR, EIF4E, and RPS6) were investigated in various organisms, that regulate physiological effects, growth factors, protein translation, cell survival, proliferation, protein synthesis, cell metabolism and survival, autophagy, fecundity rate, egg size, and follicle number, although the critical reason for longevity is still unclear in the termite castes. However, based on transcriptome profiling, the IIS-pathway-associated genes could prolong the reproductive caste lifespan and health span. Therefore, the transcriptomic shreds of evidence related to IIS-pathway genes provide new insights into the maintenance and relationships between biomolecular homeostasis and remarkable longevity. Finally, we propose a strategy for future research to decrypt the hidden costs associated with termite aging in reproductive and non-reproductive castes.


Introduction
The insulin signaling pathway (IIS-pathway) comprises all the proteins and components involved in the action of insulin within the body [1][2][3]. Insulin is the most potent physiological anabolic agent that has ever been discovered (1921) [4,5]. It is responsible for storing and synthesizing lipids, protein, and carbohydrates, and it prevents the breakdown of these macromolecules and their release into circulation [6,7]. The IIS-pathway has been thoroughly investigated in different organisms to regulate their life span [8,9]. The IIS-pathway is responsible for mediating and transducing signals across cell membranes and regulating the physiological aspects of reproduction [10][11][12].
The inhibition of upstream IIS-pathway components (daf-2 and age-1) significantly prolongs the life and health span of Caenorhabditis elegans [13]. On the other hand, IGF-I and IGF-II monitor growth and growth factors (GF) [13]; nutrient sensing, stress replication [14], cell functions, and metabolic regulation is gradually deteriorating due to aging [14,15]. Similarly, IIS-related genes are also responsible for activating hormones and the motility of cells [16,17].
Researchers previously reported that IIS-associated genes that prolong the lifespan of Drosophila melanogaster [18,19], dauer state in C. elegans [20,21], and reproductives diapause in Culex pipiens mosquito [22]. It also indicated that squirrels live for twenty-five years, rats live for three years, and C. elegans, D. longispina, and D. melanogaster live only for a few weeks [23,24]. Long-lived queens of insects, including Ephemera simulans [23], Pogonomyrmex owyheei [25], and Lasius niger can live up to ten years [26]. Polistes canadensis (a social insect) has a lifespan and genetic variation up to the genus level within a single colony [27]. However, there is still a significant gap between what we know and what we do not know regarding the enormous variety of aging rates among different animals and the mechanisms that might be responsible for this diversity [28]. Similarly, the royal caste of termites lives for 18-30 years, ultimately making termites an emerging model for longevity [29] and therefore requires further attention to explore its molecular pathways. However, there are no appropriate mechanisms available for determining the life span of a termite reproductive group [30,31].
We followed the IIS-pathway (PI3K-Akt are the main effector pathways in IIS-signalling) in insects [12,[32][33][34] to determine the IIS-pathway in primary and secondary worker reproductives (many subterranean termite colonies consist of former worker termites or immatures that have developed into larvae and eventually become wingless reproductives to supplement the colony), and non-reproductive subterranean termites [29,35,36]. Longevity in Reticulitermes chinensis was determined through next-generation sequencing (NGS) or transcriptomic [36]. A detailed evaluation of transcriptome sequences was examined with developmental processes, biological, and physiological changes in R. chinensis castes [29]. The de novo transcriptome assembly of Reticulitermes chinensis, IIS-pathway genes are enriched, contributing to termite longevity and health. In addition, the reported IIS-pathway and longevity genes in secondary worker reproductive kings (SWRK) and queens (SWRQ), primary kings (PK) and queens (PK), and male (WM) and female (WF) workers were analyzed to determine their life extension mechanism and co-evolutionary process. To better understand the genetic basis for longevity, healthy lifespan, and growth regulation in social insects and other insects, this study presents a comprehensive analysis of the IIS-pathway-associated genes implicated in R. chinensis longevity.

COG, KEGG, and GO Ontology Classifications
In the functional COG classification with 25 categories, additionally annotated 36,681 unigenes, the general functional prediction had 5707 unigenes (the colossal group), followed by cell motility (97) and nuclear structure (93), the minuscule group (Table 1). We have mapped the unique sequences into KEGG to understand the biological pathways of R. chinensis, including the mitogen-activated protein kinase (MAPK) pathway and the PI3K-Akt pathway, the two pathways that are involved in insulin signaling. However, the PI3K pathway is the most prominent insect insulin signaling pathway. From total castes, 93,078 unigenes were analyzed in IIS-pathway using transcriptome sequencing, and 343 genes were allocated to KEGG pathways ( Figure 1A-E). All DEGs in this study were mapped to terms in the GO database and determined the functions of the differentially expressed genes. Among 59 GO terms, the level2 GO enrichment analysis classified genes according to their molecular function (MF), cellular components (CC), and biological process (PB). The results indicated that PK-vs-SWRK castes exhibited significantly up-regulated unigenes. The categorical presentation suggests that a maximum number of up-regulated genes are reported in BP (cellular processes), followed by CC (cell and cell part), and MF (catalytic activity) ( Figure 1A). Among PQvs-SWRQ castes, the level2 GO terms were reported significantly as down-regulated. The maximum number of down-regulated are cellular processes at BP, followed by cell and cell part in CC and binding at MF ( Figure 1B).
Furthermore, level 2 GO terms of SWRK-vs-SWRQ individuals are additionally considerably reported as down-regulated genes shown at BP ( Figure 1C). The level2 GO terms for WF-vs-SWRQ and WM-vs-SWRK were significantly down-regulated. The categorical cluster indicates that the BP possessed cellular processes, followed by CC (cell and cell part) and catalytic activity at MF ( Figure 1D,E). from PQ-vs-PK (2277 up-regulated and 832 downregulated), followed by SWRQ-vs-SWRK (55 up-regulated and 272 downregulated) and WM (13 up-regulated and 124 downregulated) ( Figure 3). In DEGs with a p-value < 0.05, the expression levels of GO and KEGG pathways showed a significantly increased. The total number of CDs was 64,046 sequences ( Figure S5). The transcriptome sequence analysis and annotations of R. chinensis caste provided valuable evidence for evaluating all unigenes ( Figure 4A-E). Transcriptome review showed high expression genes cognate to the IIS-pathway in SWRK, SWRQ, PQ, PK, WM, and WF castes.

Discussion
The transcriptome sequences of R. chinensis castes were compiled to determine the longevity-related genes that contribute to extending the lifespan of SWRK, SWRQ, PK, and PQ caste's. Transcriptomically, we confirmed that several genes (Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6) in R. chinensis had extended the lifespan like D. melanogaster [18,19] and C. elegans [20]. The SWRK and SWRQ live longer (18-30 years) than WM and WF (a few weeks to months) [19,33]. It is challenging to sort out and estimate the lifespan of termites. Therefore, we evaluated transcriptomic sequences with next-generation sequencing (NGS) associated with developmental, biological, and physiological changes in cells or tissues [29,38]. The transcriptome research was determined through quantitative RT-qPCR. However, the current study is the second endeavor to longevitycognate genes of R. chinensis castes [29,38]. According to the transcriptome data, 35 DEGs involved in the longevity-related genes; Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6 are involved in the longevity of invertebrates and vertebrates. The insulin/insulin-like growth factor signaling pathway is a hormonally mediated cell-signaling pathway, that involves insulin-like peptides, transmembrane receptors of their cognate cell surface, and downstream effectors [19,33], which shreds of evidence annexing genetic and biochemical Figure 6. Log2 fold changes in the insulin signaling pathway in R. chinensis castes SWRK, SWRQ, PK, PQ, WM, and WF. Differential expression of genes in each group is shown as log2 fold changes compared to the reference group. The y-axis represents the log2 fold changes, and the x-axis shows the gene name and individual caste name.

Discussion
The transcriptome sequences of R. chinensis castes were compiled to determine the longevity-related genes that contribute to extending the lifespan of SWRK, SWRQ, PK, and PQ caste's. Transcriptomically, we confirmed that several genes (Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6) in R. chinensis had extended the lifespan like D. melanogaster [18,19] and C. elegans [20]. The SWRK and SWRQ live longer (18-30 years) than WM and WF (a few weeks to months) [19,33]. It is challenging to sort out and estimate the lifespan of termites. Therefore, we evaluated transcriptomic sequences with next-generation sequencing (NGS) associated with developmental, biological, and physiological changes in cells or tissues [29,38]. The transcriptome research was determined through quantitative RT-qPCR. However, the current study is the second endeavor to longevity-cognate genes of R. chinensis castes [29,38]. According to the transcriptome data, 35 DEGs involved in the longevity-related genes; Pdk1, akt2-a, Tsc2, mTOR, EIF4E, and RPS6 are involved in the longevity of invertebrates and vertebrates. The insulin/insulin-like growth factor signaling pathway is a hormonally mediated cell-signaling pathway, that involves insulin-like peptides, transmembrane receptors of their cognate cell surface, and downstream effectors [19,33], which shreds of evidence annexing genetic and biochemical changes [39]. The PI3K-Akt signaling pathway plays a principal role in insects' longevity [12,33,34].
Downstream Akt recruited and activated phosphorylation at S473 and T308 at the plasma membrane [40]. Permitted Akt controls cell survival, proliferation, and protein synthesis [41], and activates mTORC2 and Pdk1 [42,43]. In many organisms like yeast, C. elegans, and D. melanogaster, Pdk1 is important in cell survival and evolution [44,45]. Additionally, for murine embryonic development, Pdk1 is essential; mice without a Pdk1 gene died at 9.5 days of the embryo and showed abnormalities in different tissues [46]. Pdk1 hypomorphic mice have smaller cardiovascular organ volumes, and Pdk1 conditional abstraction in muscle cells leads to heart failure and minimizes lifespan [47]. Pdk1 has direct effectors on Akt, S6K, and RSK, causing embryonic stem cells to activate all three of these kinases [46]. However, Akt/PKB activations are involved in several cellular replications to magnification factor signalings, such as apoptosis bulwark, glucose transporter (GLT), and glycogen synthesis [10], while in D. melanogaster, it regulates lifespan, reproductive status, growth, and metabolism [48][49][50]. Moreover, Pdk1 in neuronal IIS-pathway can influence chemotaxis and learning [50]. PI3K-Akt and FOXO proteins are also essential in Drosophila, Nilaparvata lugens, and Sogatella furcifera flight muscle genes [51,52]. In C. elegans, Pdk1 controls insulin physiological effects, and growth factor (GF) enters the staggering dauer phase and elongates their life span when Pdk1 becomes dormant [53]. The expression of the Pdk1 gene through transcriptome and RT-PCRs analysis was significantly higher in SWRK than in any other castes. As a consequence of these findings, our results indicated that Pdk1 serves the same active function in the life span regulation of R. chinensis termites as it does in C. elegans, S. furcifera, Drosophila, and N. lugens.
Tuberous sclerosis (TSC) and its domains (Tsc1 and Tsc2) are responsible for autosomaldominating disorders, including epilepsy, skin, retina, heart, kidney, and the central nervous system [54,55]. A germ-line mutation in the rat homologous human Tsc2 gene makes the Eker rat prone to many tumors and death in the intermediate stage. Tetracycline-dependent conditions and overexpression of Tsc2 inhibited the proliferation of an Eker rat-derived kidney tumor cell line [55]. Tsc1 in the Drosophila homolog is associated with the cell mutation study (mosaic screens) and polyploidy intended to trigger a cell switch [56]. Current RT-qPCR findings in this study showed that Tsc2 is abundantly conserved among all the castes.
Protein kinase B and its three paralogs, akt1, akt2, and akt3, exhibit intensively deliberated cancer and metabolism [41], and also regulate the plasma membrane [57,58] and cell size mutations in Drosophila [56]. Akt activates considerable tissue growth/size and is essential for the PTEN facility to act as a tumor suppressor in Drosophila [10,59]. The expression of akt2-a may also help longevity and remove tumor cell incursion and metastasis in reproductive castes of R. chinensis. mTOR was first identified and cloned within the budding yeast Saccharomyces cerevisiae [60,61]. Extracellular signaling activates the mTORC1 and mTORC2 network in innate immune cells [11,62,63] and regulates translation, protein synthesis, mRNA translation, anabolic cell growth, and metabolism [64]. Knock-down of mTOR regulates low fertility in N. lugens and activates vitellogenin gene (Vg) expression in A. aegypti [65,66], also essential for transduction alimentation during mosquito egg development [67][68][69]. mTOR also controls ovary size in Drosophila [70], which promotes the number of follicles and owns diphenic development in Apis mellifera [71]. The previously reported results indicate that mTOR increases life span, fecundity, and GF in reproductive termite castes. The same strategy (RT-PCR analysis) is used for mTOR in WM, WF, PK, PQ, SWRK, and SWRQ, where it shows a possibility of prolonging the lifespan of termite castes.
The RPS6 deficiency cells did not significantly reduce global protein translation or 5 -TOP mRNAs [9]. However, functional analysis and mutations in eIF4E-1 result in embryonic defects and lethality [72,73]. Our results suggest that IIS-related pathway genes increase reproductive castes' life expectancy and control a metabolic pathway in D. melanogaster, R. chinensis, C. elegans, and fly. In addition, their nutritional reactions and sensory compensation for IIS-pathways and insulin peptides are similar [13,74,75].
In the current investigation, the hypothesized findings (RT-qPCR results) were obtained from the reproductive castes (PK, PQ, SWRK, and SWRQ), with much higher expression than in the non-reproductive castes (WM and WF). Our findings ultimately shed new insight into the maintenance and significance of biomolecular homeostasis in reproductive and nonproductive R. chinensis castes. As previously described by a number of researchers, IIS pathway-related genes desperately promote growth factors, physiological processes, cell metabolism and survival, autophagy, fecundity rate, egg size, and follicle number. These findings imply that a relatively conserved protein in the insulin signaling system significantly prolongs or avoids age-related diseases, processes, structures, and associated pathways. Further research is required to determine the biological function of genes connected to the insulin signaling pathway using genetic tools such as RNA interference and transgenic structures in order to determine how these genes contribute to the continued existence of termites. In the end, the findings we were able to obtain new shreds of insight into the processes involved in preserving biomolecule homeostasis and its connection to remarkable longevity.

Collection and Rearing of Reticulitermes Chinensis
Termite Reticulitermes chinensis colonies were dug from Chengdu in April-May 2014 ( Table 2) and transferred to the laboratory via plastic boxes (25 cm × 18 cm × 15 cm). Isolated colonies were reared for 6 years under generous (munificent) and crowded conditions (the lab was 24 h open) at room temperature (25-28 ± 1 • C) at Northwest University (http://english.nwu.edu.cn/ (accessed on 7 June 2021), Xian, Shaanxi, China. A total of 250 monohybrid colonies were established; each colony corresponds to male × female alates 78% (195/250) and female × female alates 22% (55/250) using pine sawdust (50-60% humidity) in specially designed plastic boxes (80 mm × 65 mm × 40 mm) [76]. At an early stage, active and mature colonies were selected for the experiment, from where we accumulated R. chinensis reproductive (PQ and PK; SWRQ and SWRK) and non-reproductive (WF and WM) castes [29]. These collections were done randomly from different colonies (not specific to the quantities and targeted colonies), where we found the primary and secondary reproductives and non-reproductives castes. Anyhow, many factors influence the division of labour in an insect community, including the size of the colony, castes and instar demography, food availability, and competition [77]. For example, the size of the Rhytidoponera metallica (Smith) ant and Coptotermes formosanus [78] colony influences the existence of age polyethism [79]. Colonies of R. chinensis grow from a mated pair (30-40 days start egg laying) to an immature colony (7-25 months), a juvenile colony (3-6 years), and a sexually mature colony (after 6 years). In a single colony, each caste performs its own role (e.g., food foraging and colony carrying). Due to the fact that, during the rearing of the colonies, a variety of colonies died due to fungus attacks (paralyzing the termites till death), while some colonies laid eggs very late (an average period for eggs laying was observed from 30-40 days). According to these fundamental reasons, the termites were collected randomly.

Experimental Samples
During the study, R. chinensis colonies were subjected to synchronization procedures immediately after six years of rearing of these colonies because we were looking for selected castes, i.e., primary king (PK) and queen (PQ), secondary worker reproductive king (SWRK) and queen (SWRQ), worker male (WM) and female (WF). Total RNA was extracted from the whole body (for RT-PCR, we used heads of individuals free from protozoans) of PQ, PK, SWRQ, SWRK, WF, and WM for RNA Illumina sequencing [29]. Three technical (e.g., PK-1, PK-2, and PK-3) and five biological replicates (e.g., PK-1 had five individuals) (randomly collected from 250 colonies) of each caste (six in number "PQ, PK, SWRQ, SWRK, WF, and WM") were used.

Total RNA Extraction, cDNA Synthesis, and Illumina Sequencing
The whole body of R. chinensis castes was micro-dissected (into the head, legs, thorax, and wing (of an adult) to extract DNA and RNA), and it was promptly stored in liquid nitrogen.
We used TRIzol reagent and an Agilent 2100 Bioanalyzer in order to collect sufficient RNA for Illumina sequencing (Agilent Technologies, Palo Alto, CA, USA). The Lysate RZ (TRIzol reagent) and RNAsimple Total RNA-Kit (Tiangen Biotech "Beijing" Co., Ltd., Beijing, China) were applied to the tissue to homogenize and prevent the tissue from degrading [76]. Tissue homogenates were heated to temperatures ranging from 15-30 • C to separate the nucleic acid-protein complex. Afterward, the tissue homogenate was transferred to RNase-free centrifuge tubes (for centrifuge spinning) to separate the aqueous phase, precipitation, deproteinized, and abstracted residual liquid. An adequate quantity (for Illumina sequencing and RT-PCR experiments) of total RNA was extracted and stored at −80 • C for further experiments. Finally, with the help of a spectrophotometer (NanoReady "Model: F-1100", Shanghai, China), the total amount of RNA (protein: 260/280 and salt: 260/230) was checked to quantify and verify the entire amount of RNA integrity [29].
After that, the NEB-Next Prep-Kit for Illumina (NEB) sequences [80] was used to revert cDNA in accordance with the manufacturer protocol. Furthermore, we used the QiaQuick PCR extraction kit for cDNA fragments emasculated. The poly(A) end-paired were tailed and combined with the sequence of the Illumina adapter. We obtained sequence readings of an average length (561 bp) and a minimum length (201 bp) by using the Illumina HiSeqTM 4000 platform ( Figure S1). The ligation product was selected by size from the Gene Denovo Biotechnology Co (Guangzhou, China). In addition, the Trinity system (trinityrnaseq r2012-04-27) was used for each test simultaneously in order to obtain an adequate number of clear unigenes [81,82].

Transcriptome Assembly and Reads Mapping
Raw data evaluation affects the quality and screening; therefore, we used fastp (version 0.18.0) before and after filtering and using the following parameters (¿10% of unknown and ¡50% inference of reads) for readings that contain adapters of low-quality q-value (≤20 bases) [83]. Quality clean reads of the transcriptome de novo was assembled with the short-read assembly through the trinityrnaseq r2012-04-27 assembly program [81]. Further data on the development of transcriptomic assemblies using a de Bruijn graph algorithm from short-read sequences have been provided in supplementary data [84].

Read Alignments, Normalization, and the Gene Expression Level
Short sequences and readings have been mapped into a reference SOAPaligner/soap2 tool [85,86]. The edgeR (Bioconductor package version 2.4.0) was used to generate read counts for genomic features and summarize short reads (edgeR is a method that is based on the weighted mean of log ratios). This method normalizes data in the beginning by calculating size and normalization factors from raw reads of the Illumina TM sequence. The program implements precisely established statistical approaches for multigroup experiments [87]. These genes were tenacious by the quantifications allocated exclusively to the exon region per million mapped reads (RPKM) [76]. The R-package (http://www.rproject.org (accessed on 15 March 2021)) was used to calculate all statistical data expression and visualization [37].

Differentially Expressed Genes (DEGs), Identification, Validation, and Functional Enrichment Analyses
After reading alignments and normalization of genes, finally, we obtained DEGs, using DEseq2 (because each caste has three replicates); the relative expression level of each gene was calculated [88]. These statistical methods distinguish between digital gene expression data (electronic data) of a negative binomial distribution model, fold changes (FC > 2), and the false discovery rate (FDR < 0.01) to justify the threshold level (p-value). The value of log2ratio ≥ 1 is an absolute value in calculations of significant DEGs between pooling samples. Annotations from the Gene Ontology (GO), the Kyoto Genes and Genomes Encyclopedia (KEGG) (http://www.genome.jp/kegg (accessed on 15 March 2021), the Clusters of Orthologous Groups (COG) (https://www.ncbi.nlm.nih.gov/COG/ (accessed on 15 March 2021) and the Non-redundant (NR) (http://www.ncbi.nlm.nih.gov (accessed on 15 March 2021) database were all relegated [86]. Blast2GO software assigned the GO annotation (obtained from the Nr database profile) [89]. In addition, WEGO software was used to classify unigenes functions [90].

Statistical Data Analysis
A standard method 2 −∆∆Ct (formulae are given below) was used to calculate the relative gene expression of the mRNA genes used with each replicate [92,93]. The RT-qPCR experiment quantifying and counting results were statistically analyzed using IBM Corp. Released 2019 (IBM SPSS Statistics for Windows, Version 26.0. Armonk, NY, USA: IBM Corp.) to determine the relationship between groups (reproductive and non-reproductives) using a t-test [29]. All values were expressed as mean ± standard deviation (p-values < 0.05) and statistical figures constructed with Microsoft Excel (0365) and OriginPro (2018).