Long live the queen, the king and the commoner? Transcript expression differences between old and young in the termite Cryptotermes secundus

Social insects provide promising new avenues for aging research. Within a colony, individuals that share the same genetic background can differ in lifespan by up to two orders of magnitude. Reproducing queens (and in termites also kings) can live for more than 20 years, extraordinary lifespans for insects. We studied aging in a termite species, Cryptotermes secundus, which lives in less socially complex societies with a few hundred colony members. Reproductives develop from workers which are totipotent immatures. Comparing transcriptomes of young and old individuals, we found evidence for aging in reproductives that was especially associated with DNA and protein damage and the activity of transposable elements. By contrast, workers seemed to be better protected against aging. Thus our results differed from those obtained for social insects that live in more complex societies. Yet, they are in agreement with lifespan estimates for the study species. Our data are also in line with expectations from evolutionary theory. For individuals that are able to reproduce, it predicts that aging should only start after reaching maturity. As C. secundus workers are immatures with full reproductive options we expect them to invest into anti-aging processes. Our study illustrates that the degree of aging can differ between social insects and that it may be associated with caste-specific opportunities for reproduction.


Introduction
"Mors certa-hora incerta" (Anselm of Canterbury): all living organisms die at some point due to extrinsic factors (extrinsic mortality) and/or intrinsic aging processes (intrinsic mortality, senescence). Aging is associated with time passing by and for most living organisms with a progressive deterioration of the body, which becomes apparent by a decline in motor and cognitive performance, a collapse of physiological processes and a decline in the ability to deal with various sources of stress [1][2][3]. Ultimately, evolutionary theory explains aging by a decline of the force of natural selection with age [4][5][6][7]

Sample collection, preparation, sequencing and library construction
For this study, C. secundus colonies were collected near Palmerston-Channel Island from a mangrove ecosystem (12˚30' S, 131˚00' E; Northern Territory, Australia) [38]. The colonies were transferred to wooden blocks of Pinus radiata. After transport to the laboratory, the stock colonies were kept under standardized conditions in a climate room with a temperature of 28˚C, 70% humidity and a 12h day/night rhythm. Under these conditions colonies develop as in the field (see [34]). Total RNA was extracted from twelve different individuals belonging to nine different colonies: four primary queens, four primary kings and four workers (two young and two old individuals of each caste, S1 Table). Young queens and kings were collected one year after the nuptial flight (yearling reproductives). The age was determined by the color of the cuticle, tunnel size, and colony size, which was approximately 20 workers and one soldier. The old primary reproductives had been reproducing for at least seven years. Young workers were selected from the youngest worker instars, no more than six months old. The age of the old workers was at least three years, based on the instar (last larval worker instar) and the minimum developmental time it takes to reach this instar [34]. They might be older due to further stationary molts (size remains the same and no morphological changes take place) or regressive molts (decrease in size and wing-bud development). The sex of C. secundus workers can currently not be determined; there are no reliable morphological differences between male and female workers and attempts with traditional karyological methods failed.
An in-house protocol was followed for RNA extraction. Briefly, individuals were placed on ice and the gut was removed and discarded. Whole bodies were then used for RNA extraction. Samples were transferred into peqGOLD Tri Fast (PEQLAB) and homogenized in a Tissue Lyser II (QIAGEN). Chloroform was used for protein precipitation. From the aqueous phase, RNA was precipitated using Ambion isopropyl alcohol and then washed with 75% ethanol. Obtained pellets were solved in nuclease-free water. DNA was subsequently digested using the DNase I Amplification Grade kit (Sigma Aldrich, Cat. No. AMPD1). We performed an RNA Integrity Number Analysis (RIN Analysis) measuring the RNA concentration with the Agilent RNA 6000 Nano Kit using an Agilent 2100 Bioanalyzer (Agilent Technologies) for quality control. Samples with total RNA were sent on dry ice to BGI Tech Solutions (HONGKONG) Co. and then to Shenzhen, PR China for sequencing. The preparation of the cDNA libraries was performed by BGI according to their internal and proprietary standard operating procedure. The company performed paired-end sequencing (not-strand specific) on Illumina HiSeq platforms (S1 Table).

RNA-seq data analysis
The raw reads' quality was evaluated with FastQC (v0.11. 4, [44]). Common Illumina-and BGI in-house adapter sequences were trimmed from the raw reads using Trimmomatic (v0.33, [45]) and we kept only reads with a minimum length of 70 bp (for a rationale see [46]). Expression differences were assessed at the level of the transcripts. Reads were pseudo-aligned with Kallisto (v0.43.0, [47]) against C. secundus transcriptome obtained from a draft version of the C. secundus genome (provided as supplementary material on Dryad) [48]. The counts estimated by Kallisto were rounded to the nearest integer and used to check the completeness of our samples with BUSCO (v. 3.0. 2 [49]).The protein coding sequences of the transcripts that had at least one positive match (n>0 in count table) were used as input for BUSCO. We ran BUSCO (at the protein level -m prot) against the insect gene set of ortholog groups (insecta odb 9) with default settings. This gene set includes 1,658 single copy orthologs.
The Kallisto counts were used as input for DESeq2 (v1. 10.1, [50]) in R (v3.2.3, R Core Team 2015) to determine transcript expression differences between old and young individuals of queens, kings and workers. In DESeq2 the p-values are adjusted for multiple testing using the false discovery rate (FDR) approach [51]. In order to correct for unaccounted sources of variation the 'Surrogate Variable Analysis' (sva) package (v 3.20.0) as implemented in R was used [52,53]. The software identifies and estimates surrogate variables for unknown sources of variation (for instance, batch or colony effects). For data visualization, a principal component analysis (PCA) was performed with DESeq2 using transformed count data (variance stabilization). The age-specific differentially expressed transcripts (DETs) were compared between castes and the overlaps were visually represented with Venn diagrams generated using the online tool Venny (v2. 1 [54]) and graphically processed with Inkscape (v0.91, www.inkscape. org). ). For all species, we only kept the longest isoform per gene and additionally removed sequences with internal stop codons and/or selenocysteines (U). The latter causes problems in downstream analyses and sometimes sequences with selenocysteine are considered to be potential pseudo-genes (see [48]). COGs were inferred among the four reference species using the software OrthoFinder (v. 0.4.0 [64]).

Functional annotation and enrichment
A functional enrichment analysis (GO enrichment) was done with DAVID (v6. 8, [65]). For this we made a second set of COGs (see above), but only with C. secundus and D. melanogaster, and additionally a BLASTP search of C. secundus sequences against the protein coding sequences (longest isoforms only) of D. melanogaster with a threshold E-value of 1e-05. The set of homologs used for the GO enrichment consisted of (i) single copy 1:1 orthologs and (ii) homologs (not single copy 1:1 orthologs) found via the BLASTP search and filtered using the best bit scores. germanica (BGER) (sequences provided on Dryad for each COG). COGs were found using the text search option available in OrthoDB (search for gene name or IDs e.g. "AGO2", "FBgn0087035"). We aligned the sequences separately for each COG with MAFFT (v7.294b), with the G-INS-i algorithm default settings [67]. We built a profile hidden Markov model (HMM) from each multiple sequence alignment using hmmbuild (HMMER v3.1b2). We then used the HMM to search with hmmsearch against C. secundus and M. natalensis protein coding sequences to identify candidate sequences for each COG. For inference of gene trees, we only kept sequences with a threshold E-value of �1e-40. All candidate sequences were searched with hmmscan against the Pfam database (Pfam A, release 30).
Phylogenetic trees of closely related genes were inferred for a total of 26 genes of interest provided on Dryad. We included sequences of all species mentioned above that passed the threshold (�1e-40). In all sequences, selenocysteine (U) was replaced by "X". Sequences were aligned as described above. We subsequently identified ambiguously aligned sequence regions with Aliscore (v. 2, [68,69]) (settings: -r: maximal number of pairwise comparisons) and subsequently removed those sections with Alicut (v. 2.3, https://www.zfmk.de/en/research/researchcentres-and-groups/utilities, masked alignments provided on Dryad). We used IQ-Tree (v. 1.5.5, [70]) to infer phylogenetic gene trees (single genes or groups of closely related genes) using the Maximum-Likelihood approach separately for genes (gene groups) of interest. Using parsimony start trees, we estimated the best model with the implemented ModelFinder [71] for available standard models, including the free rate models LG4M and LG4X [72]. We used default settings for rates, number of rate categories, and the Bayesian Information Criterion (BIC) to estimate the best substitution model. Statistical support was inferred from 2,000 nonparametric bootstrap replicates. We visualized the unrooted trees with bootstrap support using Seaview (v4.5.4, [73]) and graphically processed the trees with Inskape (v0.91, www.inkscape. org) (provided in Dryad).
All raw sequencing reads are deposited on EMBL (Primary Accession PRJEB27153, for sample accession numbers see S1 Table). Additional supplementary data can be accessed from the Dryad repository (dryad link can be provided upon acceptance).

Comparison with other species
The results obtained in this study were compared with similar aging studies for the fruit fly D. melanogaster In addition, C. secundus DETs were compared with D. melanogaster 'aging genes' available in the ageing gene database (GenAge Build 18, [75]).

Differential expression between age-classes within castes in C. secundus
The results of the BUSCO analysis showed for each sample a completeness of more than 90% of the insect gene set of OGs (insecta odb 9) (see S1 Table). Using transcript transformed count data, the PCA separated the castes on the first component, which explained 35% of the variance (Fig 2). The queens clustered together and were clearly separated from the other castes, whereas the kings clustered closer to the workers. The second component explained 18% of the variance and separated the age classes. The age class pattern was reversed in the workers in comparison to the reproductives (Fig 2). Two surrogate variables were estimated with the package sva [53], and they were included in the formula for differential expression testing to correct for unaccounted sources of variation. A total of 815 transcripts were significantly differentially expressed (DETs) between the age-classes within castes: 193 for the queens, 248 for the kings, and 374 for the workers (padj < 0.05, S2 Table and Fig 3). These caste-specific age-related DETs were divided into two groups: transcripts that were more highly expressed in younger individuals than in older ones (HY) and transcripts that were more highly expressed in old individuals compared to young ones (HO) (Fig 4). Only one HO DET, annotated as a carboxylesterase like gene, was shared across all castes (see Fig 3); 10% of the DETs were shared in a pairwise manner. Hence, most DETs were caste specific. DETs in the workers were involved in a plethora of different functions including metabolic processes, growth, development and morphogenesis, regulation of transcription, alternative splicing and chromatin remodeling (S2 and S3 Tables).

Differential expression of genes linked to aging
The DETs that were considered important for aging were split into the following categories: (i) DNA damage response, genome stability and telomeres; (ii) transposable elements (TEs); (iii) oxidative stress; (iv) neural aging; (v) reproduction, (vi) immunity and (vii) other age-related genes.
Transposable elements. Seven HY and thirteen HO DETs in queens were linked to TEs. Out of these, two HY and three HO were shared with the kings (Fig 4A and 4B). One of the shared HY DETs was identified as 'pelota' (pelo), a gene that encodes a conserved protein involved in transposon silencing [87]. Kings had in total eight HY and seven HO DETs related to TEs. One of the HY DETs in kings was identified as a homolog of vasa, a gene involved in transposon silencing in the germline by participating in piRNA biogenesis and amplification [88], and another one was identified as 'Heat shock protein cognate 4' (Hsc70-4), which is involved in the RNA interference pathway (RNAi) and in heterochromatin formation [89,90]. Strikingly, 21 DETs in workers were associated with TEs and contrary to the reproductives most of these DETs were more highly expressed in young compared to old individuals: 19 HY and two HO (S2 Table). One of the HY DETs in workers was identified as 'eggless', which is an H3K9 methyltransferase that interacts with heterochromatin protein 1 (HP1) to spread heterochromatin, and it is necessary for TE silencing via piRNA pathway [91,92]. Of the two HO DETs in workers, one was annotated as 'Vacuolar H+ ATPase 16kD subunit 1' (Vha16-1). Vha16-1 participates in the uptake of dsRNA to the cell, which allows the activation of the RNAi pathway [93].
Oxidative stress. In queens and kings, one HY DET was identified as darkener of apricot (Doa) (Fig 4A), which among other functions can negatively regulate the expression of superoxide dismutases (SOD1 and SOD2) in D. melanogaster [94]. Additionally, in kings the GO enrichment revealed for the HY DETs multiple GO BP terms related to catabolism of xenobiotic compounds and insecticides (S3 Table). The signal for these BP terms originated from seven HY DETs annotated as cytochrome P450 family 6 (CYP6). Another HY DET in kings was identified as a catalase-like gene. Other DETs in kings potentially involved in oxidative stress were identified as cytochrome P450s family 4 (CYP4, one HY and two HO), cytochrome P450 of unknown family (HY, homologue to CYP303a1), and one mitochondrial uncoupling protein (HY), which can modulate and reduce ROS production [95]. In workers multiple transcripts were associated to oxidative stress: two catalase-like proteins (both HO, one of them shared in opposite direction with kings, see S1 Fig), a homolog of 'inactivation no afterpotential E' (HO, inaE in D. melanogaster)", a carbonic anhydrase (HO), isocitrate dehydrogenase (HO) and five cytochrome P450s from different families (one HY and four HO, see S2 Table). P450 genes, however, can also play an important role in caste differentiation in termites and hence need not be linked with oxidative stress.
Neural aging. In queens three HY DETs were potentially linked to aging of the brain: Palmitoyl-protein thioesterase 1 (Ppt1) and two Down syndrome cell adhesion molecules (Dscam) [96,97]. One of the Dscam proteins was shared with workers but the change in expression with age was in opposite direction (HY in the queens but HO in the workers, S1 Fig). In kings, one HY DET was identified as a homolog of apolipoprotein D (ApoD). ApoD is mainly expressed in the brain and nervous system and its over expression in D. melanogaster increases lifespan and stress resistance [98,99]. Additionally in kings, the 1:1 ortholog of refractory to sigma p (ref (2) S3 Table). This signal came from multiple transcripts including pelo, hsc-4 and nap1 that are involved in transposon silencing, heterochromatin formation and genome stability (see Transposable elements).  Immunity. Several DETs between age-classes were related to immunity. Two DETs were shared between queens and workers in the same direction (HY, Fig 4), and were identified as a 'thaumatin-like' protein and a 'termicin family' protein. These proteins are involved in immune defense against fungal infections [104-106]. Queens and kings shared the differential expression of 'modular serine protease', which is involved in innate immunity [107,108], but in opposite direction (HY in kings and HO in queens). In kings, two more HO DETs were linked to immunity: an antimicrobial peptide Attacin, and a protein containing a 'single domain Von Willebrand factor type C' (SVWC). Proteins with the SVWC domain are involved in response to bacterial infections and viral challenges (see [109]). In workers two genes involved in innate immune response coding for a transglutaminase and a peptidoglycan recognition protein [110,111] were significantly differentially expressed (both HO).
Other age-related genes. The 1:1 ortholog of ATP citrate lyase (ATPCL), a gene that in D. melanogaster males has a pro-aging effect [112], was upregulated in old reproductives (HO, Fig 4). In contrast to ATPCL, the 1:1 ortholog of 'ATP synthase, subunit D' (ATPsynD), whose downregulation in D. melanogaster with low carbohydrate to protein diet has been associated to a lifespan extension [113], was downregulated in old kings (HY).
In old workers, the 1:1 ortholog of 'S-adenosylmethionine Synthetase' (Sam-S) was downregulated. The decrease of S-adenosylmethionine has been linked to a lifespan extension in C. elegans and D. melanogaster [114,115].

Comparison with other species
The results for the age-associated expression differences in C. secundus were compared with the results of similar aging studies for reproducing females of a solitary insect, D. melanogaster,  (Table 1). In workers ten GO terms were shared with D. melanogaster females in opposite direction and nine GO terms were shared in the same direction. One GO term of interest was 'oxidationreduction process' (oxidative theory of aging) (see Table 1 Table). In the reproductives, four HY DETs (four transcripts in queens and two of these shared by kings, same gene) were identified as homologs of Ecdysone-inducible gene L2 (Imp-l2), which is a gene contained in the GenAge database.

Comparison with the termite M. bellicosus
The GO enrichment results of C. secundus can only be compared with the enrichment results of M. bellicosus major workers as very few genes were differentially expressed between young and old M. bellicosus reproductives [18]. C. secundus queens shared three enriched GO terms with M. bellicosus major workers two in opposite direction and one in the same direction (Table 1). C. secundus workers shared with M. bellicosus major workers ten GO terms in opposite direction and 32 in the same direction ( Table 1).Most of the GO terms were related to development, morphogenesis and metabolism.
The comparison using C. secundus DETs (without considering GO terms) was done against M. bellicosus DE genes (old vs young) of queens, kings, minor workers and major workers. Elsner et al. (2018) [18] found only two DE genes in M. bellicosus kings, and no overlap was found with C. secundus DETs. C. secundus queens shared a 1:1 orthologue with M. bellicosus queens, a cytochrome b5 (HO) (S4 Table). Additionally, the queens shared 17 genes with M. bellicosus major workers (1:1 orthologs to Drosophila in both species, S4 Table). Five genes showed an opposite expression pattern (HY in C. secundus queens and HO in M. bellicosus major worker) of which one was identified as an aging marker in Drosophila: cora [40,41] (S4 and S6 Tables). Cora was also upregulated in old minor workers of M. bellicosus (HO). The function of cora is not known but it interacts with forkhead box, sub-group O (foxo, IIS pathway) and rictor (mTOR pathway) [116]. Three genes were upregulated in C. secundus queens and M. bellicosus major workers (HO in both species) and nine genes were downregulated (HY in both species). Workers shared 46 genes (more than half functionally annotated) with M. bellicosus major workers, 24 in the same direction (15 HY and 9 HO) and 22 in opposite direction (all 1:1 orthologs to Drosophila in C. secundus, S4 Table).

Discussion
We identified many age-related expression changes in C. secundus reproductives, which live in less complex colonies. This differs from the higher termite M. bellicosus, for which only very few genes were differentially expressed in the heads of young and old reproductives [18]. Our results imply that C. secundus reproductives have an earlier onset of aging than those of M. bellicosus, as supported by longevity data. Macrotermes reproductives can live for 20 years and more, while the maximum lifespan of C. secundus reproductives is around 10-13 years [22]. In contrast, the totipotent C. secundus workers showed fewer signs of aging in line with the hypothesis that in organisms with full reproductive options, aging is expected to start only after reaching maturity. For the present study we assumed that transcript expression differences reflected changes at the protein level (as most transcriptome studies do). This is known to be the case for certain proteins in D. melanogaster (e.g., heat shock proteins and prophenoloxidase) [117] and in bees (e.g. vitellogenin [118]), but results need to be taken with caution because transcript and protein abundances do not necessarily correspond [119][120][121].

Aging in C. secundus compared to other social insect species
Age-associated expression changes in C. secundus can be related to senescence processes, such as TE-activity and its resulting DNA and protein damage. In queens and kings, DETs related to DNA damage repair and genome stability were downregulated in old compared to young individuals (S2 Table). This also applied to old workers but in addition they had two ATR homologs upregulated in old individuals (S2 Table). ATR coordinates DNA damage responses (repair and cell cycle checkpoint signaling), participates in telomere maintenance and contributes to genome stability [76][77][78]. This implies that the protection against damage decreased with age in the reproductives but less so in workers.
This view is also supported by our results for TE activity. TEs have been associated with aging in a broad range of animals, from C. elegans, D. melanogaster, and mice to humans [17, [122][123][124][125][126] and most recently in termites [18]. Old queens had more TE-related transcripts than young queens and associated with this was a downregulation of a 1:1 ortholog (D. melanogaster) encoding pelo, which is involved in TE silencing [87]. In kings the number of TE related transcripts was similar in old compared to young kings (see S2 Table) but three genes involved in TE silencing were downregulated with age (HY): pelo and the homologs of vasa and Hsc70-4 [87][88][89][90]. The reverse was found for workers in which the very low TE activity (two DETs) in old workers was accompanied with an upregulation of Vha16-1, which is necessary for the activation of the RNAi pathway and a systemic response against TEs [93]. In line with TE activity, was the downregulation of ATPCL (Acetyl-CoA metabolism) in the reproductives (HY). This gene has been connected to a disregulation of histone acetylation levels and the disruption of heterochromatin formation (high TE content) [112,127].
The TE-related transcripts upregulated in young workers might not be related to aging but developmental processes. Some upregulated TE-related transcripts in young workers were LINE (long interspersed nuclear elements) retrotransposons (Dfam results, S7 Table). LINE-1 retrotransposons have been shown to regulate global chromatin accessibility in the early embryo in mice [128]. Though it still requires further studies, we propose that the TE related activity in young C. secundus workers reflects regulation of complex postembryonic development of these immature instars, in line with other development-related worker DETs and the GO enrichment results (S2 and S3 Tables).
Despite the differences in tissue specificity (C. secundus: whole body without gut vs M. bellicosus: head only), which might increase the number of false negatives (see [129]), our current results compared with those for the termite M. bellicosus [18] reveal some striking similarities and differences. While the same aging-related mechanisms were detected (TE-activity and piRNA mediated TE defense), the aging pattern seems reversed between castes. Old major M. bellicosus workers were characterized by a high TE activity whereas reproductives in this 'higher' termite seem to be protected against TEs by the piRNA pathway [18]. This difference can be explained by the fact that the major workers of M. bellicosus are sterile, while C. secundus workers are totipotent immatures from which all reproductives derive (Fig 1B). Evolutionary theory predicts that in organisms with full reproductive options aging starts only after the onset of maturity [8]. Hence, we expect that the totipotent C. secundus workers are selected to invest in anti-aging processes, while this is less the case for sterile M. bellicosus workers in which all reproduction is channeled through a separate sexual line (Fig 1A).
Many other factors differ between the two termite species, most importantly lifestyle and associated with it, colony size and pathogen load. C. secundus is a drywood termite that nests in a single, non-decomposed piece of wood which the workers never leave to forage outside. Hence, like other drywood termites [130], C. secundus probably has a low pathogen load for workers and reproductives, and it has small colony sizes and workers (and reproductives) that experience low extrinsic mortality risks [131]. By contrast, M. bellicosus belongs to the foraging termites in which workers leave the nest to collect food and bring it back. Associated with this lifestyle are-besides reduced reproductive options for workers-larger colony sizes, and increased pathogen loads and extrinsic mortality risks, especially for workers [22]. Extrinsic mortality risks are the most important factors to influence aging in non-social organisms where all individuals can reproduce [4,8]. Hence, we also would predict faster aging of workers due to increased extrinsic mortality risks in M. bellicosus. However, these models do not consider social organisms and how things change with sterility. Models including all these factors are warranted to determine the separate contribution of different factors (e.g. worker sterility, colony size, and extrinsic mortality) for the evolution of lifespans in social insects. In termites, these factors form 'syndromes' associated with life type. Foraging termite species are always socially more complex, and have workers with reduced reproductive capacities and shorter lifespans than wood-dwelling species [22]. Hence, these traits seem to co-evolve and it will be difficult to separate and test the contribution of single factors.
An upregulation of genes, such as catalases and SODs, involved in the protection against oxidative stress has been associated with the long lifespan of queens of the termite Reticulitermes speratus in studies that compared neotenic queens and workers of unknown age [32,33]. This species lives in colonies with a degree of social complexity and worker reproductive options intermediate between M. bellicosus and C. secundus. Its developmental system is complicated and lifespan data have not been published. As a foraging species, it has bifurcated development, hence workers have reduced reproductive options but they seem to be not fully sterile. Similar to C. secundus, founding queens (primary queens) can be replaced in R. speratus. However, they seem to be rather short lived and replacements are produced from nymphal instars via parthenogenesis [132]. We also found evidence that oxidative stress plays a role in aging in C. secundus, but it supports the view that non-sterile workers, in contrast to the reproductives, increase their protection with age: catalase (two DETs), inaE, carbonic anhydrase and idh-all genes important in oxidative stress defense [133][134][135][136]-were upregulated in old compared to young workers. The expression of one catalase (upregulated in old workers) was downregulated in old kings and the 1:1 ortholog of doa was downregulated in old reproductives. Doa negatively regulates the expression of SOD1 and SOD2 [94]. Although cytochrome P450s from the families CYP6 and CYP4 were differentially expressed in workers, these genes might be connected to caste differentiation rather than oxidative stress [137,138], and for kings, other functions for the differentially expressed CYP4 can also not be excluded.

Neural aging
Neurodegeneration and the decline in cognitive functions is characteristic of aging in many species [119,139,140]. However, such a decline with age was not observed in the ant Pheidole dentata: old workers neither experienced a deterioration of cognitive functions nor changes in the brain characteristic of aging [141]. These results question whether neural deterioration is a hallmark of aging in social insects. Although in our study no tests were performed regarding the decline of cognitive performance and sensory perception in old workers, the higher expression of odorant receptors (four DETs, see S2 Table) and the expression of idh involved in the protection against oxidative stress damage in dopaminergic neurons [136], might be a sign of no aging in C. secundus workers. In contrast, we consider the upregulation with age (HO) of ref (2)p in kings and the downregulation of Ppt1 (HY) in queens as clear signs of nervous system degeneration and aging. Ref (2)p is a conserved marker of neuronal aging [100,142], and the deficiency of Ppt1 in flies is associated with an abnormal accumulation of lipofuscin in the nervous system (lipofuscin accumulates with age) and a reduced lifespan [96]. Some neurodegenerative diseases in humans are caused by mutations in Ppt1 [143].

The role of immunity in aging
The loss of immunocompetence and the dysregulation of the immune system are considered a hallmark of aging [144][145][146], but how aging and immunity influence each other is still not completely understood [147]. In D. melanogaster, an overexpression of immunity genes is characteristic of aging [148][149][150]. An immunity response clears pathogens but it can also cause tissue damage via inflammation. Thus a hyper-activation of the immune system can represent a state of immunopathology and promote aging [146]. Overall, all castes in C. secundus showed an increase in expression of immunity genes with age, but to properly interpret these results the caste-specific interplay between immunity and oxidative stress response should be considered. Immune defense is associated with an upregulation of oxidative stress response genes [151][152][153][154]. Hence, the overexpression of immunity genes and the downregulation of stress and oxidative stress response genes in C. secundus reproductives (but not workers) could reflect the lack of homeostasis between pathways, and be eventually interpreted as a sign of aging.

Typical aging pathways in C. secundus
We found few but relevant DETs connected with typical aging pathways like the IIS or the TOR pathway: Imp-12, prmt1 and cora, which were all over-expressed in young compared to old queens (Imp-12 also in young kings). These genes have been described to directly or indirectly interact with foxo (IIS pathway). Imp-l2 binds insulin like peptides (ILPs) and its overexpression in D. melanogaster leads to a lifespan extension [155,156]. Thus a decline in expression of Imp-l2 in old C. secundus reproductives might again be a sign of aging. Prmt1 encodes an arginine methyltransferase, which, in mammals, methylates foxo [157,158]. The function of cora is unknown, yet for D. melanogaster, experimental evidence exists that it physically interacts with foxo (IIS pathway) and rictor (TOR pathway) [116]. ATPsynD, differentially expressed in kings (HY), was another gene connected to the TOR pathway and involved in aging. The downregulation of this gene in D. melanogaster conferred a lifespan extension to females fed with a low carbohydrate-to-protein diet [113]. The effect of this gene in the context of termites and their diet should be explored in future studies. The expression patterns of Imp-l2, prmt1, ATPsynD and cora in C. secundus suggest that, as in other insects/animals lifespan determination and aging processes might be modulated by the typical aging pathways IIS and TOR.

Conclusions
Our results demonstrate that reproductives of the lower termite C. secundus show signs of senescence. The age related changes in expression suggest that aging might be linked to TE activity, oxidative stress and 'wear and tear', and that it may partly be modulated by the IIS pathway, immunity response and epigenetic modifications. For the totipotent workers, which can become reproductives, we did not find evidence of aging but rather a strong signal of metabolism, growth and development as indicated by the identity of DETs and the GO enrichment analysis. Hence our results contrast strongly with those for the higher termite M. bellicosus where reproductives do not show signs of aging while the major workers did [18]. This is in line with the general prediction of life history theory for organisms with reproductive options that aging should only start after having reached maturity [8]: the totipotent C. secundus workers are immatures (they are larval instars) that can develop into reproductives and hence should not age. M. bellicosus workers are also non-adults, but in contrast they are irrevocably sterile and experience much higher extrinsic mortality risks which should both favor faster aging. The degree to which workers can reproduce differs across social insect species and generally correlates with sociality [159]. Therefore, we suggest that aging in reproductives and workers depends, among other factors, on the degree of worker's reproductive options which should be tested in future studies.

S1 Fig. Shared differentially expressed transcripts (DETs) with age (old vs young) in opposite direction between C. secundus queens, kings and workers.
The heatmap depicts the log2fold changes in expression; each row corresponds to a caste and each column to a transcript. NS: differential expression not significant (padj>0.05). (EPS) S1