- Open Access
Disentangling the aging gene expression network of termite queens
BMC Genomics volume 22, Article number: 339 (2021)
Most insects are relatively short-lived, with a maximum lifespan of a few weeks, like the aging model organism, the fruit-fly Drosophila melanogaster. By contrast, the queens of many social insects (termites, ants and some bees) can live from a few years to decades. This makes social insects promising models in aging research providing insights into how a long reproductive life can be achieved. Yet, aging studies on social insect reproductives are hampered by a lack of quantitative data on age-dependent survival and time series analyses that cover the whole lifespan of such long-lived individuals. We studied aging in queens of the drywood termite Cryptotermes secundus by determining survival probabilities over a period of 15 years and performed transcriptome analyses for queens of known age that covered their whole lifespan.
The maximum lifespan of C. secundus queens was 13 years, with a median maximum longevity of 11.0 years. Time course and co-expression network analyses of gene expression patterns over time indicated a non-gradual aging pattern. It was characterized by networks of genes that became differentially expressed only late in life, namely after ten years, which associates well with the median maximum lifespan for queens. These old-age gene networks reflect processes of physiological upheaval. We detected strong signs of stress, decline, defense and repair at the transcriptional level of epigenetic control as well as at the post-transcriptional level with changes in transposable element activity and the proteostasis network. The latter depicts an upregulation of protein degradation, together with protein synthesis and protein folding, processes which are often down-regulated in old animals. The simultaneous upregulation of protein synthesis and autophagy is indicative of a stress-response mediated by the transcription factor cnc, a homolog of human nrf genes.
Our results show non-linear senescence with a rather sudden physiological upheaval at old-age. Most importantly, they point to a re-wiring in the proteostasis network and stress as part of the aging process of social insect queens, shortly before queens die.
Almost all animals age, but at a different pace . The fruit fly Drosophila melanogaster lives only for around seven weeks , while the clam Ocean Quahog, Arctica islandica, can have a lifespan of more than 400 years , and the giant barrel sponge Xestopongia muta can live more than two millennia . Generally, organisms with large differences in aging rates are found between widely divergent species [1, 5], which makes controlled comparisons of the underlying aging mechanisms difficult.
Classical model organisms typically have a short lifespan and can be characterized by r-life history strategies (‘live fast, have many offspring and die young’) as exactly these traits make them good model organisms. Social insects such as termites, ants, or the honeybee, offer promising new insights into aging because individuals with the same genetic background can differ by orders of magnitudes in lifespan. Within a social insect colony, which is generally a large family, the reproducing queen (and in termites, also the king) can reach lifespans of more than 20 years, while non-reproducing workers have a lifespan of a few months only [6,7,8,9]. However, quantitative demographic data covering the whole lifespan of queens are inherently rare (for ants: [10,11,12] and references therein; for termites: [13, 14]) and many reports on queen-longevity are more anecdotal. Thus, it is largely unknown for long-lived queens whether they age gradually or whether aging is a more sudden event.
During recent years, several pioneering studies, especially on the honeybee, revealed exciting new insights into the mechanisms of how queens can live so long. Generally, the TOR (target of rapamycin) and the IIS (insulin/insulin-like growth factor1 signalling) pathway have been associated with longevity in model organisms from D. melanogaster to mice and humans [15,16,17]. They are the most intensively studied aging-related pathways and they have also been associated with caste differences in social Hymenoptera (e.g.,[18,19,20,21,22,23]). Additionally, in the honeybee, juvenile hormone (JH) seems to have lost its direct gonadotropic function in adults so that queens have a high expression of vitellogenin (Vg), which encodes yolk precursors, without requiring high JH titers (e.g., [24, 25]). This result has led to the hypothesis that an uncoupling between JH and Vg expression might account for the long life of honeybee queens , as well as social insect queens more generally , because the life-time shortening consequences of high JH titers are absent. However, this re-wiring along the JH-Vg axis is not universal for all social Hymenoptera since the queens of many ant and bee species require JH for vitellogenesis (e.g.  and references therein). For termites, fewer studies exist but JH is required for vitellogenesis [28, 29] and a recent study revealed that no re-wiring exists along the JH-Vg axis . Hence, other mechanisms must exist to explain the long life of termite queens. Studies of the subterranean termite Reticulitermes speratus implicated the involvement of a breast cancer type 1 susceptibility (BRCA1) homolog , which is involved in DNA repair , and better protection against oxidative stress by superoxide dismutases and catalases [33, 34]. The latter has also been discussed for other social Hymenoptera, including ants and the honeybee. Yet, the overall evidence of the role of oxidative stress is less clear (reviewed in [35,36,37,38]). Furthermore, regulation of the activity of transposable elements (TEs)  and changes in the insulin/insulin-like growth factor1 signalling (IIS) and target of rapamycin (TOR) pathways  have been linked with caste-specific aging differences in termites Yet, all studies on social insects suffer from a lack of time-series data to investigate molecular changes across the lifespan of long-lived queens. Like the demographic life history data, such data are inherently difficult to obtain due to the long lifespan of queens. However, they are necessary (i) to understand the aging process, (ii) to work out potential changes compared to solitary insects, and (iii) to identify the relevant age-classes for detailed studies. The latter is a completely overlooked issue but highly relevant. Differences across studies might be consequences of non-comparable age-classes between studies, if, for example, aging is a non-linear process.
We studied aging in termite queens of known age across their entire lifespan to measure at the ultimate, eco-evolutionary level age-dependent survival and at the proximate, mechanistic level age-specific changes in gene expression. For the latter, we generated head / thorax transcriptomes of queens of different ages (for an outline of the workflow, see Additional File 1, Figure S1; for rational of tissue choice, see Methods). We used field collected, newly established colonies of the wood-dwelling (i.e.., one-piece nesting) termite Cryptotermes secundus (Hill, 1925) (Blattodea, Isoptera, Kalotermitidae) that were kept under identical conditions in the laboratory for a time span of up to 15 years. Laboratory conditions have been adjusted to C. secundus so that colony development (e.g. fecundity, growth rates, colony composition, molting types) does not differ compared to that of field colonies as has been shown in earlier studies [40, 41]. Keeping colonies under constant good conditions without external mortality allowed us to study intrinsic aging, disentangled from causes of extrinsic mortality such as predation, food shortage, or disease. As typical for wood-dwelling termite species, C. secundus colonies are founded in a piece of wood, which serves as food and shelter and which workers never leave to forage outside. Such species have a low social complexity with small colonies and totipotent workers that develop into sexuals.
Our survival analysis covered a time period of up to 15 years; yet none of the queens that would have had an age of 14 and 15 years was alive and the maximal age we recorded was 13 years (Additional File 2, Table S1). In fact, most queens with an expected age of ≥ 12 years were dead. Out of eight queens in this ‘old-age’ class, only a single queen (13 years) had survived (Fig. 1; Additional File 2, Table S1). Kaplan Meier survival analysis estimated the median longevity of the queens to be 12.0 years (SE: ± 0.54) (mean longevity: 11.1 years, SE: ± 0.66) in the laboratory after successful colony foundation (Fig. 1).
Identifying transcripts that change their expression with age: age‐related DETs
To study gene expression changes over the life-time of queens, we generated transcriptomes of head / thorax from twelve queens with different chronological age since the onset of reproduction, from two until 13 years, covering the complete lifespan of C. secundus queens: 2, 3, 4, 5, 6, 7, 8, 9, 10 (two samples), 11, and 13 years (Additional File 2, Table S2). The queens used for gene expression analyses came from the same data set as those for the survival analysis; they were alive queens that entered the survival analysis as censored data (for more details see Methods).
A total of 169 transcripts were significantly differentially expressed (DETs) over time as revealed by Iso-MaSigPro time series analysis (Additional File 2, Table S3). According to their expression pattern, DETs were grouped into six Iso-MaSigPro clusters (hereafter, ‘cluster’) (Fig. 2). Cluster 1 represented 44 DETs, which were slightly expressed in young queens followed by a decline at middle ages and a strong increase when queens became older. The 32 DETs of cluster 2 characterized young queens with a declining expression with age. Clusters 3 and 5 comprised 31 and 37 DETs, respectively, that were highly expressed in middle-aged queens, while cluster 4 and cluster 6 (15 and 10 DETs) characterized old queens with no expression in young ones. Thus, in the following text, we refer to the DETs as young (cluster two), middle-aged (clusters three and five) and old DETs (clusters one, four and six). Details for all clusters are provided in Additional File 2, Table S3.
Identifying modules of co‐expressed transcripts
To identify modules of co-expressed transcripts, we performed a weighted gene co-expression network analysis (WGCNA). It revealed a total of 254 modules of co-expressed transcripts. Based on eigengene values, 13 modules correlated significantly positively with age and 13 negatively (see Additional File 1, Figures S2 and S3; Additional File 3 (WGCNA module-age association, shown are eigengene values for all modules).
Identifying transcript co‐expression modules with age‐related DETs
Within the age-correlated WGCNA modules, we identified age-related DETs. The negatively age-correlated module ‘seashell4’ had the highest number of young DETs (10 DETs). No gene ontology (GO) term was enriched for this module. The highest number of old DETs was found in the positively age-correlated modules ‘cyan’ (89 DETs) and ‘tan’ (79 DETs) (Additional File 2, Table S4 and S5). Only broad categories were enriched in the ‘cyan’ module (e.g., RNA metabolic process and gene expression). The ‘tan’ module was enriched for ribosomal and tRNA related functions (Additional File 1, Figure S4).
Extracting age‐related subnetworks based on age‐related DETs
To generate subnetworks related to the age-related DETs, we located them in the WGCNA co-expression network. These DETs and their one- and two-step neighbors (i.e., the ‘second level neighborhood’) were then extracted from the co-expression network, which resulted in 50 subnetworks of different sizes (for more details, see Methods) (Additional File 1, Figure S5). Note, DETs might be located at the boundaries of multiple WGCNA modules, which means the subnetworks obtained consist of fragments of multiple WGCNA modules. The resulting subnetworks either contained young and middle-aged DETs or old DETs, with a single exception where a middle-aged DET was in the periphery of the largest subnetwork containing old DETs. The largest subnetwork containing either young and middle-aged DETs (hereafter, young subnetwork) or old DETs (plus a single middle-age DET; hereafter old subnetwork) were further analyzed.
The largest young subnetwork comprised 164 transcripts (out of these 24 Iso-MaSigPro DETs), of which only 12 (7 %) were one-to-one orthologs to D. melanogaster genes (Additional File 2, Table S4). The GO enrichment analysis of the young subnetwork showed multiple Biological Process (BP) terms related to RNA catabolism, but these GO terms were not significant after correcting for multiple testing (Additional File 2, Figure S6).
TE activity and genome instability
53 DETs (32 %) of the young subnetwork were related to TEs (Fig. 3 and Additional File 2, Table S4), comprising TEs and genes from TE defense pathways. This included one homolog of the gene argonaute 2 (ago2) (two transcripts), an essential gene of the endo-siRNA pathway which silences TEs , and arsenite 2 (ars-2), which is required for siRNA and miRNA-mediated TE silencing . Additionally, we found two genes connected to DNA damage response and genome instability: kin17 and PIF1-like gene.
From well-known aging pathways, we identified (i) inositol polyphosphate phosphatase 2 (mipp2) and (ii) adenylyl cyclase 76E (ac76E). The former is part of the TOR pathway and has been associated with longevity , and the latter is activated by the transcription factor ‘Forkhead box O’ (foxo). Additionally, we found several fecundity-related DETs. They included two transcripts of the gene hu li tai shao (hts) (one a DET of IsoMaSigPro cluster two) and one homolog of the gene bällchen (ball) (two transcripts) (one a DET of cluster five, Fig. 3).
The largest ‘old subnetwork’ comprised 1,098 transcripts (out of these 42 Iso-MaSigPro DETs). 521 transcripts (47 %) were identified as one-to-one orthologs of D. melanogaster genes (Additional File 2, Table S4). Iso-MaSigPro DETs in the old subnetwork belonged mainly to Iso-MaSigPro clusters 1 and 4. The second level neighborhoods of these DETs were connected in the network, and a GO enrichment analysis revealed multiple GO terms associated with protein-related functions, including translation, protein folding, unfolded protein binding, proteolysis involved in cellular protein catabolic process, protein targeting to ER, ribosome, and proteasome complex (Additional File 1, Figures S7 and S8). 198 transcripts of the old subnetwork (18 %) were involved in protein translation, protein folding, and protein catabolism and proteolysis [45, 46] (Figs. 4 and 5, and Additional File 2, Table S5). Additionally, 61 transcripts (~ 6 %) were related to TEs (Additional File 2, Table S5).
Epigenetic modifications, transcriptional regulation, and TE activity
Many old subnetwork genes are involved in de/acetylation and methylation of DNA, which are important epigenetic modifications that regulate gene expression and genome stability [47,48,49] (Additional File 2, Table S5).
Most strikingly, two crucial histone acetylation modifying complexes, the Tip60 acetyltransferase complex and the male specific lethal (msl) complex were represented in the old subnetwork. The former included the genes dom, ing3, mrg15, pont and rept, and the latter msl-1, msl-3 and mof. Genes involved in deacetylation of DNA were, for instance, sirtuin 1 (sirt1), histone deacetylase 3 (HDAC3), and histone deacetylase 6 (HDAC6). Genes linked to epigenetic histone methylation included, for instance, ash-1 and lid. Another well-represented group of genes connected to expression regulation in the old subnetwork were spliceosome components and splicing factors. Additionally, we found in the old subnetwork important transcripts related to TE silencing: dicer-2, Hsc70-4, Hsc70-3, Hsp83, trsn, armi, Rm62, Gasz, Tudor-SN, and Hel25E. Details are given in Additional File 2, Table S5.
Proteostasis and oxidative stress
Related to proteostasis, we detected a strong signal for protein synthesis and degradation. Regarding protein synthesis, the old subnetwork comprised many transcripts coding for initiation, elongation and termination factors, as well as many ribosomal proteins and aminoacyl-tRNA synthetases (Fig. 4; Additional File 2, Table S5). Regarding protein degradation, almost all subunits of the ubiquitin proteasome system (UPS) were present (Fig. 5), as well as autophagy genes, heat shock proteins, and the transcription factor xbp1. Xbp1 is involved in the ‘unfolded protein response’ (UPR) and the ER-associated protein degradation (ERAD) pathway [50, 51].
Additionally, BRCA1 was also present in the old subnetwork. This gene is involved in oxidative stress response, and in the transcriptional activation of proteasomal genes by stabilizing the transcription factor cnc/nrf-2 (cap-n-collar/nuclear factor erythroid 2–related factor 2) . Other genes in the old subnetwork involved in oxidative stress response and transcriptionally activated by nrf-2 were thioredoxin and S-glutathione transferase.
Additionally, the old subnetwork was characterized by a signature of ecdysone biosynthesis with ecdysone receptor (EcR), ecdysone-induced protein 75B (Eip75B), phantom and disembodied. The presence of Phosphatidylinositol 3 kinase 68D (Pi3K92E) links to the IIS pathway.
Locating age‐related co‐expression modules in the age‐related subnetworks
Finally, we inspected those WGCNA modules with a large fraction of transcripts in the young and old subnetworks and those modules that were significantly associated with age.
In the young subnetwork, WGCNA modules with a large fraction of transcripts were ‘saddlebrown’ and ‘skyblue4’, which both did not significantly correlate with age. Significantly age-correlated co-expression modules were firebrick2, indianred1 and seashell4 (Additional File 1, Figure S2). No GO terms were significantly enriched for any of these modules.
In the old subnetwork, the modules with a large fraction of transcripts were ‘green’ and ‘paleturquoise’, which did not significantly correlate with age. The old subnetwork contained transcripts of 13 significantly age-correlated co-expression modules (Additional File 1, Figure S3). The GO enrichment analysis of these modules revealed several terms involved in protein-related functions, including ribosome biogenesis, rRNA processing, protein folding, translation, unfolded protein binding, protein catabolic process, protein transport, tRNA aminoacylation for protein translation, and proteasome core complex (Additional File 1, Figure S4, S9, S10 and Additional File 4, Table S6).
Our study revealed a median maximum reproductive longevity of C. secundus queens of 12 years with a maximum lifespan of 13 years when excluding all causes of extrinsic mortality in the laboratory (Fig. 1, Additional File 2, Table S1). The small difference between the median and the maximum recorded lifespan reflects the rather sudden decline in life expectancy after an age of around 11–12 years. Of eight queens with a potential age ≥ 12 years, all had died, except one 13-year old queen. The survival curve indicates a type I survivorship with high age-specific survival probabilities until midage and a rapid decline in survival later in life, after queens successfully founded a colony and without extrinsic mortality. Including the early colony founding stages, which are characterized by very high mortalities with more than 99 % of the dispersing sexuals dying in C. secundus , this suggests a bathtube curve of mortality for queens. A high early failure period is followed by a stable failure period and a final wear-out failure period .
Our transcriptome study identified six clusters of transcripts that were significantly differentially expressed with age (DETs) (Fig. 2): one cluster for young queens (cluster 2), two for medium-aged queens (cluster 3 and 5) and three for old queens (cluster 1, 4, and 6). This implies that three ‘molecular’ life stages can be distinguished in C. secundus queens, with the third corresponding to old, aged queens that will probably die soon as no queen reached a lifespan beyond 13 years. The co-expression network analysis, which extracted subnetworks based on age-related DETs, resulted in two main subnetworks, a young and an old subnetwork. This implies that there are two age-related ‘molecular’ life stages, as DETs/genes of young and medium ages belonged to the same young subnetwork.
Compared to the old subnetwork, the young subnetwork contained relatives few DETs (164 vs. 1,098 in the old subnetwork) and they were characteristic for young and medium ages. This shows the similarity in expression of these two age stages. Not unexpectedly, the young subnetwork indicates an upregulation of transcripts linked with fecundity (e.g., hts, ball) and of the TOR pathway (mipp2) which has been associated with longevity . The upregulation of Ac76E may imply that the IIS pathway is down because this gene is activated by the transcription factor foxo, which is inhibited by an upregulated IIS pathway. However, other evidence suggests that, like in other social insects (e.g.  and references therein), queens are characterized by an upregulated IIS pathway [30, 39]. Additionally, we detected several upregulated TEs-related transcripts associated with signs of an upregulation of the endo-siRNA pathway (e.g., ago2, ars), which is a transcriptional and post-transcriptional TE-defence mechanism of the soma [42, 43, 55].
In an earlier study on C. secundus, we compared one-year old reproducing queens with queens that were at least seven years old, though the exact age of the later was unknown . We identified 193 DETs between both age classes with no signs of physiological upheaval as revealed in the current old subnetwork. This implies that the old queens used in the earlier study were not physiological old and rather classify as queens of medium age. This is in line with the gene functions identified in the earlier study . Similar to the current study, we had found an IIS- and TE-related signal as well as some fecundity-related genes, the latter were higher expressed in the medium- than young-age queens.
The old subnetwork contained many more transcripts (1,098 vs. 164 in the young subnetwork). Our results imply a physiological stage of upheaval shortly before queens die. There were strong signs of decline and repair at the upstream transcriptional level of epigenetic control as well as at the post-transcriptional level of TE-activity and the proteostasis network.
An upregulation of genes modifying histone marks implied considerable epigenetic changes that lead to altered gene expression as is typical for aging organisms:
First, our results indicate dynamic changes of ‘active’ histone marks of euchromatin because many genes related to H3K4 and H3K36 de/methylation and H4K16 de/acetylation were present in the old subnetwork (Additional File 2, Table S5). For instance, the Tip60- as well as the msl-complex were well represented. Both complexes are involved in the acetylation of histones, including H4K16, which, for instance, is indicative for replicative aging in yeast . The upregulation of Sirt1 may function as an antagonist as it deacetylates H4K16. The same applies to the deacetylases HDAC6 and HDAC3, which can also deacetylase histones .
Second, there is also evidence for changes of repressive histone mark of heterochromatin (e.g., linked to H3K9 and H3K27 acetylation) (Additional File 2, Table S5). For instance, the old age transcript ALP1 is an antagonist of HP1, the latter is involved in the maintenance of heterochromatin. HP1 generally decreases with age, and its overexpression can lead to an increased lifespan in D. melanogaster [57, 58].
In line with a deregulation of repressive histone mark of heterochromatin, several TE-related transcripts were also connected with the old subnetwork (Additional File 2, Table S5). TEs often are accumulated in heterochromatin, which silences their activity . Yet, dysregulation of heterochromatin can allow TEs to become active and this has been associated with aging [60, 61].
Also, the upregulation of several genes from two TE-defense pathways - the endo-siRNA pathway (e.g., dicer2, Hsc-70-4, trsn) as well as the piRNA pathways (e.g., piRNA biosynthesis. armi, gasz, Hel25E, Rm62; ping-pong cycle: Tudor-SN, qin) - support the notion of active TEs. Both pathways silence TEs posttranscriptionally [62,63,64]. TE activity and the piRNA pathway have also been associated with aging and termite reproductives’ longevity in another termite species .
Loss of proteostasis
Our results revealed a very strong proteostasis signal indicative of an upheaval in protein synthesis, protein folding and protein degradation in old queens. Many genes from the proteostasis network were detected in the old subnetwork (Additional File 2, Table S5). They indicate an upregulation of protein degradation, together with protein synthesis and protein folding (Figs. 4 and 5). This is unusual because old organisms are typically characterized by a downregulation of all these processes (Fig. 6).
Many transcripts coding for ribosomal proteins and aminoacyl-tRNA synthetases were found in the old subnetwork, indicative of upregulated protein synthesis (Additional File 2, Table S5). This is further supported by a strong signal of an active TORC1 system which promotes protein synthesis [65,66,67] (Fig. 7). Thus, for instance, many downstream eukaryotic initiation factors (eIF4A, eIF4B, eIF4E, eIF4G) and eukaryotic elongation factor (eEF2) transcripts were found in the old subnetwork (Fig. 4). They activate ribosome biogenesis, translational elongation, and cap-dependent translocation (Fig. 7).
Normally, an active TORC1 system is associated with a downregulation of protein degradation as it inhibits proteolytic systems [66,67,68] and autophagy (e.g., upregulated TORC1 inhibits ATG1, which is necessary for autophagy activation; Fig. 7). Surprisingly, however, we found strong evidence of upregulated protein degradation in the old subnetwork. Several transcripts linked to autophagy, almost all subunits of the UPS, the UPR-, and the ERAD pathway as well as heat shock proteins characterized the old subnetwork (Figs. 5 and 6; Additional File 2, Table S5).
Linking protein synthesis and degradation
The simultaneous upregulation of protein synthesis and autophagy might be explained by a stress response. In D. melanogaster and humans, upregulated TORC1 enhances an oxidative stress response under stress conditions, controlled by the transcription factor cnc, a homolog of human nrf genes (Fig. 7). Ubiquitinated proteins and damaged mitochondria activate cnc/ nrf-2 via p62, supported by upregulated TORC1, which then activates oxidative stress response genes [69, 70]. Additionally, cnc is known to activate chaperones (protein folding) and the proteasome , and this has been associated with lifespan expansion in D. melanogaster and Caenorhabditis elegans [72, 73].
Support for the conclusion of a stress-associated, active cnc transcription factor in old queens comes from several transcripts in the old subnetwork (Fig. 7): (i) BRCA1, which indirectly actives cnc by inhibiting cnc inhibitor Keap1, and (ii) the Tip 60 complex as well as genes that are transcriptionally activated by cnc, such as thioredoxin, S-Glutathione transferases and UPS genes (Fig. 7).
Hence, our results imply that old C. secundus queens are in a stage of stress. They have mounted stress response systems mediated by cnc, including protein degradation and protein folding. However, it is unusual that old queens can do this. In D. melanogaster, only young individuals can mount this stress response . The constant activation of the proteasome in these very old queens may lead to their death (note, the studied queens had reached their maximum lifespan, we never had older queens) as the proteasome’s constant activation in transgenic flies was detrimental for survival .
Oxidative stress in other social insects
There has been a debate about the role of oxidative stress to explain the long lifespan of social Hymenoptera queens, yet the evidence is inconclusive (e.g., reviewed by [35,36,37,38]). For instance, markers of oxidative stress in honey bee workers’ brains do not increase with age, although they live shorter than queens . In the ant Lasius niger, both workers and males show higher activity of the antioxidant superoxide dismutase than queens, but both live shorter than queens . These studies have shown that higher expression of oxidative stress response genes does not necessarily correlate with longer lifespans. A recent comparative analysis including one termite, one ant and several bee species, also did not obtain consistent results of an association between oxidative stress and longevity . For the termite Reticulitermes speratus, it has been suggested that queens are better protected against oxidative stress as qRT-PCRs studies showed a higher expression of the antioxidants catalases and peroxiredoxins in queens compared to workers , while kings were characterized by a high expression of BRCA1 (in the fat body) compared to workers . Unfortunately, the age of the studied reproductives is not known. It would be of interest to study expression of these genes with age, as this would contribute to a better understanding of the aging process, but such studies are rare.
Our results imply that C. secundus queens do not age gradually. At old age, there is a physiological stage of upheaval, characterized by signs of stress (activity of TEs, active crc), defence (piRNA pathways) and repair (protein degradation and synthesis) before the animals die.
This sudden decline is in line with our life history records for C. secundus (Fig. 1), which match survival curves for the ant Cardiocondyla obscurior . Most of the few other studies that surveyed survival of social insect colonies and queens did not cover the final stage of death, although some were long-term studies [10, 12]. These studies recorded a low, steady rate of mortality, which is consistent with the stable failure period of a bathtube curve of mortality, before a short, final wear-out failure period starts . These non-linear aging curves stress the importance of using queens of known age for aging studies as processes revealed for middle-aged versus old queens are expected to differ considerably.
Figure S1 (Additional File 1) provides a schematized workflow of the analyses described in the following sections.
Collection and colony maintenance
From 2002 until 2016, C. secundus colonies were collected from mangroves near Palmerston - Channel Island (12°30’ S, 131°00’ E; Northern Territory, Australia) when they were less than one year old . Colonies of an age of less than one year can be unambiguously identified by the size and slightly lighter sclerotization of the founding reproductives (primary reproductives), the presence of fewer than 20 workers and short tunnel systems of a few centimeters. All collected colonies were transferred to Pinus radiata wood blocks and transported to the laboratory in Germany, where they were kept under standardized conditions in a climate room with a temperature of 28 °C, 70 % humidity and a 12 h day/night rhythm. Under these conditions, colonies develop like in the field (see [40, 41]). Only colonies that survived the transfer from the field to the laboratory in Germany and the re-establishment in the new wood block were used for survival analysis.
The survival of primary queens (and kings) was determined by their presence /absence in a colony.
Founding (primary) queens can be identified by their dark brown color, compound eyes and wing abscission scars. If the primary reproductives had died and the colony was still alive, they had been replaced by neotenic replacement reproductives that lack these traits.
Survival census of colonies was restricted to the times when colonies had to be transferred to new wood blocks to provide consistently abundant food conditions (low food conditions change colony trajectories and interfere with normal development, ). We restrained from splitting wood blocks more often as splitting and colony transfer is stressful for colonies. Thus, censuses were restricted to around every third year. To ensure that time of death of the queen – and hence her age - was accurately determined (i.e., that they have died just within a few months prior to splitting), we only included colonies / queens, which clearly had newly developed neotenic replacement reproductives. Newly developed neotenic replacement reproductives can be unambiguously identified, for instance, by their very slight sclerotization.
We did not have a single cohort that we followed through time, as it is impossible to collect sufficient sample sizes of colonies with a known age of one year in the field. Hence, we used all one –year old colonies that we had collected over time to maximize sample size and which we had left undisturbed in the laboratory (except for necessary translocations to new wood blocks). These colonies were split in either 2017 or in 2018, which were the endpoints of this study. As sample size was still relatively small in 2017, we used another set of colonies that we split in 2018 to obtain more reliable survival analysis. Queens that were alive in 2017, respectively 2018, entered the survival analysis as censored data; representative samples of these queens were used for the transcriptome study. In total, we had over 100 one-year old, established colonies that entered this experiment. Out of these colonies, 41 colonies / queens could be used as the age of death could be reliably determined. Thus, the sample size was 41 queens (Additional File 2, Table S1).
The median longevity of queens was determined using Kaplan Meier survival analysis in SPSS 23  with alive queens entering the data set as censored data. Confidence intervals were based on a transformation of intervals for the log-minus-log survival function according to recommendations by Hosmer et al. . Overall, we had surviving primary queens from one year up to a maximum of 13 years; the three expected 14- and 15-year old queens were all dead (Additional File 2, Table S1).
Rational of tissue choice
We choose brain and thorax (including legs) to obtain a brain and endocrine signature, including the corpora allata and some hemolymph and fat body from the thorax. Recent work on C. secundus has pointed out that transcriptomes of the head and (pro)thorax conveys important information of queens and life history traits .
RNA extraction and sequencing
RNA was extracted from twelve queens with different chronological ages since onset of reproduction from two years until 13 years: 2, 3, 4, 5, 6, 7, 8, 9, 10 (two samples), 11, and 13 years. In colonies older than 13 years, the original queen was always replaced by a neotenic replacement queen.
An in-house protocol was followed for RNA extraction (see ). Individuals were placed on ice and the gut was removed and discarded. The head together with the thorax were 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 Beijing Genomics Institute (BGI) Tech Solutions (HONGKONG) Co. and then to the BGI-Shenzhen (PR China) for sequencing. BGI performed the preparation of the cDNA libraries according to their internal and proprietary standard operating procedure. The cDNA libraries were paired-end sequenced (not-strand specific) on Illumina HiSeq 2500 and 4000 platforms (100 base pairs read length and about 4 Giga bases per sample). Index sequences from the machine reads were demultiplexed and a quality -check and filtering of raw reads was done using the package soapuke (-n 0.05 -l 20 -q 0.2 -p 1 -i -Q 2 -G --seqType 1 and -A 0.5, http://soap.genomics.org.cn/).
Processing of RNASeq raw reads
FastQC (v0.11.4)  was used to evaluate the quality of the cleaned raw reads. To obtain a transcript count table, the cleaned raw reads were pseudo-aligned with Kallisto (default settings, v0.43.0)  against a C. secundus transcriptome obtained from a draft version of the C. secundus genome (with estimated gene and transcript models, see [39, 80]). The counts estimated by Kallisto were normalized using DESeq2 (v1.18.1, count/size factor)  in R (v3.4.4) .
Time course analysis to identify age‐associated differentially expressed transcripts (DETs)
The normalized counts were used as input for the R package Iso-MaSigPro (v1.50.0)  to test for differentially expressed transcripts (DETs) across time. Iso-MaSigPro is designed for the analysis of multiple time-course transcriptome data. It implements negative binomial generalized linear models (GLMs) [83, 84]. Significantly differentially expressed transcripts (FDR corrected p-value set to 0.05) were clustered using the clustering algorithm mclust in Iso-MaSigPro and the k.mclust option set as ‘true’. With this option, mclust computes an optimal number of clusters (k) . We obtained six Iso-MaSigPro clusters (Additional File 2, Table S2).
Weighted gene co‐expression Network analysis (WGCNA) to identify networks of co‐expressed transcripts
To obtain networks of co-expressed transcripts that were categorized as modules, we performed a Weighted Gene Co-expression Network Analysis (WGCNA). The counts obtained with Kallisto (v0.43.0)  were transformed using variance stabilizing transformation (vst) as implemented in DESeq2 (v1.18.1) . The vst transformed counts were used to perform a co-expression network analysis with the R package WGCNA (v1.63) . For more details on the methodology, see [85,86,87]. In short, (Additional File 1, Figure S1, workflow, right side), a similarity matrix was built by calculating Pearson correlations between the expression values of all pairs of transcripts. Using the similarity matrix, a signed weighted adjacency matrix was obtained as described by the formula:
Where corij is the Pearson correlation between the expression pattern of transcript ‘i’ and transcript ‘j’ (the similarity value). The value of β was chosen based on the soft-thresholding approach . With this value of β, we obtained a weighted network with an approximate scale-free topology (β = 14, scale-free topology R2 = 0.84). In a signed weighted adjacency matrix negative and small positive correlations get negligibly small adjacency values shifting the focus on strong positive correlations. Seeing the adjacency matrix as a network, the nodes correspond to the transcripts and the connections between nodes correspond to the adjacency values (transformed correlation coefficients). A topological overlap matrix (TOM), which in addition to the adjacency matrix considers topological similarity (shared neighbors reinforce the connection strength between two nodes), was constructed using the adjacency matrix . To define transcript modules, a hierarchical clustering tree was constructed using the dissimilarity measure (1-TOM). Transcript modules were defined by cutting the branches of the tree using the Dynamic Hybrid Tree Cut algorithm  and the minimum module size was set to 30 transcripts. Transcript modules with similar expression profiles were merged by hierarchical clustering of the eigengene correlation values. Briefly, a hierarchical clustering tree was created with the eigengene dissimilarity measure (1-correlation coefficient of eigengenes) and a tree height cut of 0.20 was used (corresponds to an eigengene cor ≥ 0.80). Eigengenes were calculated with the function moduleEigengenes (default settings) . A module eigengene corresponds to the module’s first principal component and can be seen as a weighted average expression profile . To find significantly associated modules with age, correlations between age and eigengenes of the merged modules were calculated. Each module was named after a color by WGCNA.
The adjacency matrix of the WGCNA was visualized using Cytoscape (v3.7.1) , only including pairs of nodes with a corij ≥ 0.90. Each module’s color corresponds with the respective module name (e.g., saddlebrown color for the saddlebrown module).
To identify co-expression modules containing age-related DETs, we looked for age-related DETs from the Iso-MaSigPro analysis in the WGCNA modules. Those modules that were significantly correlated with age and held the highest number of Iso-MaSigPro DETs were further inspected.
Identifying/Extracting age‐related subnetworks based on age‐related DETs
To identify age-related subnetworks within the co-expression network, we combined the results of the Iso-MaSigPro analysis with those from the WGCNA and extracted subnetworks that were based on age-related Iso-MaSigPro DETs. Therefore, we extracted 1st and 2nd neighbors of DETs based on the WGCNA co-expression network (i.e., the visual representation of the adjacency matrix). To do this, we used the ‘First neighbors’ function of Cytoscape. We selected an age-related DET from Iso-MaSigPro as transcript of interest. By calling the function, the neighboring transcripts were selected, which were then extracted to form a subnetwork. By calling the function twice, one obtains the one- and two-step neighbors (called ‘second level neighborhood’) of the transcript of interest. This was done for each DET identified in IsoMaSigPro.
The obtained subnetworks were clearly separated in those containing young Iso-MaSigPro DETs (young subnetworks) and those containing old Iso-MaSigPro DETs (old subnetworks). The largest subnetwork obtained for each group was used for further analysis paying attention to both transcript identity and WGCNA module content. Thus, for instance, we looked for WGCNA modules that had been identified to be age-related within the global WGCNA in these subnetworks.
The AutoAnnotate Cytoscape plug-in (v1.3)  was used to annotate the subnetworks using the clustering algorithm ‘Markov Cluster’ (MCL)  to define and visualize sub-clusters, and the labeling algorithm ‘Adjacent Words’ to label the sub-clusters. The Cytoscape plug-in BiNGO (v3.0.3)  was used for GO enrichment analysis. The p-values of the GO enrichment analysis were adjusted for multiple testing using the FDR approach . Subnetworks were graphically processed with Inkscape (v0.91, www.inkscape.org).
Transcript (functional) annotation
Nucleotide and protein sequences were obtained from the draft version of C. secundus genome [39, 80]. For annotation, the translated transcripts were searched against the Pfam database (Pfam A, release 30)  with the software hmmscan (option --cut_ga, HMMer v.3.1b2)  and against the InterPro database with the software InterProScan (v5.17-56.0) . Additionally, we did a BLASTX search (NCBI BLAST suite v. 2.3.0)  with an e-value of 1e− 05 as cut-off against the protein coding sequences of the termite Zootermopsis nevadensis (official gene set version 2.2) . To further assist the annotation, we inferred a set of clusters of orthologous sequence groups (COGs) from the official gene sets at the amino acid level of C. secundus (draft version) and D. melanogaster, and 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.
To detect possible TEs, transcripts were searched against the Dfam database (v2.0)  using nhmmer . A transcript was considered TE related if there was a hit against the Dfam database and the other annotation sources (Pfam, Interpro and BLAST) were not pointing in the direction of a known gene.
Gene identification and construction of gene trees
In addition to the functional annotation, we inferred phylogenetic trees for selected transcripts (Supplementary Archive 1, DRYAD, doi available upon acceptance, see below for a private reviewer link). Following the procedure described in , we retrieved protein coding sequences of the respective cluster of orthologous sequence groups (COGs) from OrthoDB v.9.1  for the following species: D. melanogaster (DMEL), Apis mellifera (AMEL), Cardiocondyla obscurior (COBS), Polistes canadensis (PCAN), Tribolium castaneum (TCAS), Z. nevadensis (ZNEV) and Blattella germanica (BGER). COGs were identified using text search by searching for the gene name or IDs of D. melanogaster. In case Selenocysteine (U) was included in sequences, “U” was replaced by “X“ to avoid problems in downstream analyses since many programs cannot handle Selenocysteine properly. Protein sequences of COGs of the above-listed species were separately aligned with MAFFT (v.7.294b) applying the G-INS-i algorithm and otherwise default options . For each multiple sequence alignment, a profile hidden Markov model (pHMM) was built with hmmbuild (HMMER v.3.1b2) . Then the pHMM was used to search (hmmsearch) for corresponding protein coding sequences in C. secundus and Macrotermes natalensis to identify orthologous candidate sequences for each COG in both species. For gene (transcript) tree inference, we only kept sequences with a threshold e-value of ≤ 1e− 40. In addition, we annotated all candidate sequences identified in C. secundus and M. natalensis against the Pfam database (Pfam A, release 30) using hmmscan (HMMER v.3.1b2).
To infer phylogenetic gene trees, we merged for each COG the COGs (amino-acid level) of the seven species listed above with the putatively orthologous amino-acid sequences of C. secundus and M. natalensis. We generated multiple sequence alignments for a total of 29 genes of interest applying MAFFT (G-INS-i, see above). Ambiguously aligned sequence sections were identified with Aliscore (v. 2 [104, 105]; settings: -r: 10,000,000 and otherwise defaults) and removed with Alicut (v. 2.3, https://www.zfmk.de/en/research/research-centres-and-groups/utilities); masked alignments are provided as Supplementary Archive S1 (deposited at DRYAD, doi available upon acceptance, see below for a private reviewer link). Phylogenetic trees were inferred with IQ-TREE (1.7-beta12 ) for each gene. The best model was selected with the implemented ModelFinder  from all available nuclear models implemented in IQ-TREE plus the two protein mixture models LG4M and LG4X  based on the Bayesian Information Criterion (BIC). We applied default settings for rates and the number of rate categories. Statistical support was inferred from 2,000 non-parametric bootstrap replicates. Unrooted trees with the bootstrap support mapped were visualized with Seaview (v4.5.4 ) and provided in Newick Format with Supplementary Archive S1 at DRYAD (doi available upon acceptance, see below for a private reviewer link).
Availability of data and materials
Raw sequence reads are deposited on NCBI (Bioproject PRJNA691762, BioSample Accession numbers:
Jones OR, Scheuerlein A, Salguero-Gómez R, Camarda CG, Schaible R, Casper BB et al. Diversity of ageing across the tree of life. Nature. 2014;505:169–73. doi:https://doi.org/10.1038/nature12789.
Linford NJ, Bilgir C, Ro J, Pletcher SD. Measurement of lifespan in Drosophila melanogaster. J Vis Exp. 2013;71. doi:https://doi.org/10.3791/50068.
Ungvari Z, Ridgway I, Philipp EER, Campbell CM, McQuary P, Chow T, et al. Extreme longevity is associated with increased resistance to oxidative stress in Arctica islandica, the longest-living non-colonial animal. J Gerontol A Biol Sci Med Sci. 2011;66:741–50. doi: https://doi.org/10.1093/gerona/glr044.
McMurray SE, Blum JE, Pawlik JR 2008 Redwood of the reef: growth and age of the giant barrel sponge Xestospongia muta in the Forida Keys. Marine Biol. 2008;155:159–171. doi:https://doi.org/10.1007/s00227-008-1014-z.
Rose MR. Evolutionary biology of aging. Oxford: Oxford University Press; 1991. ill.; 25 cm.
Keller L, Genoud M. Extraordinary lifespans in ants: a test of evolutionary theories of ageing. Nature. 1997;389:958–60. doi:https://doi.org/10.1038/40130.
Keller L. Queen lifespan and colony characteristics in ants and termites. Insectes Soc. 1998;45:235–46. doi:https://doi.org/10.1007/s000400050084.
Toth AL, Sumner S, Jeanne RL. Patterns of longevity across a sociality gradient in vespid wasps. Curr Opin insect Sci. 2016;16:28–35. doi:https://doi.org/10.1016/j.cois.2016.05.006.
Korb J, Thorne B. Sociality in Termites. In: Rubenstein DR, Abbot P, editors. Comparative Social Evolution. Cambridge: Cambridge University Press; 2017. p. 124–53. doi: https://doi.org/10.1017/9781107338319.006.
Cole BJ The ecological setting of social evolution: the demography of ant populations. In: Gadau J, Fewell J, editors. Organization of Insect Societies. Cambridge: Harvard Univ Press; 2009. p. 74–104.
Heinze J, Schrempf A. Terminal investment: individual reproduction of ant queens increases with age. PLoS One. 2012;7:e35201. doi:https://doi.org/10.1371/journal.pone.0035201.
Tschinkel WR Lifespan, age, size-specific mortality and dispersion of colonies of the Florida harvester ant, Pogonomyrmex badius. Insect Soc. 2017;64:285–296. doi:https://doi.org/10.1007/s00040-017-0544-0
Thorne BL, Breisch NL, Haverty MI Longevity of kings and queens and first time of production of fertile progeny in dampwood termite (Isoptera; Termopsidae; Zootermopsis) colonies with different reproductive structures. J Anim Biol. 2002;71:1030-41. doi:j.1365-2656.2002.00666.x.
Elsner D, Meusemann K., Korb J. Longevity and transposon defense, the case of termite reproductives. Proc Natl Acad Sci USA. 2018;115:5504–9. doi:https://doi.org/10.1073/pnas.1804046115
Evans DS, Kapahi P, Hsueh W-C, Kockel L. TOR signaling never gets old: aging, longevity and TORC1 activity. Ageing Res Rev. 2011;10:225–37. doi: https://doi.org/10.1016/j.arr.2010.04.001.
Partridge L, Alic N, Bjedov I, Piper MD. Ageing in Drosophila: the role of the insulin/Igf and TOR signalling network. Exp Gerontol. 2011;46:376–81. doi: https://doi.org/10.1016/j.exger.2010.09.003.
Antikainen H, Driscoll M, Haspel G, Dobrowolski R. TOR-mediated regulation of metabolism in aging. Aging Cell. 2017;16:1219–33. doi: https://doi.org/10.1111/acel.12689.
Libbrecht R, Corona M, Wende F, Azevedo DO, Serrão, JR, Keller L. Interplay between insulin signaling, juvenile hormone, and vitellogenin regulates maternal effects on polyphenism in ants. Proc Natl Acad Sci U S A. 2013;110:11050–55. doi: https://doi.org/10.1073/pnas.1221781110.
Mutti NS, Dolezal AG, Wolschin F, Mutti JS, Gill KS, Amdam GV. IRS and TOR nutrient-signaling pathways act via uvenile hormone to influence honey bee caste fate. J Experim Biol. 2011;214:3977–84. https://doi.org/10.1242/jeb.061499.
Amdam GV, Norberg K, Fondrk MK, Page RE Jr. Reproductive ground plan may mediate colony-level selection effects on individual foraging behavior in honey bees. Proc Natl Acad Sci U S A. 2004;101:11350–55. doi: https://doi.org/10.1073/pnas.0403073101.
Warner MR, Qiu L, Holmes MJ, Mikheyev AS, Linksvayer TA. Convergent eusocial evolution is based on a shared reproductive groundplan plus lineage-specific plastic genes. Nat. Commun. 2019;10:2651. doi: https://doi.org/10.1038/s41467-019-10546-w.
Chandra V, Fetter-Pruneda I, Oxley PR, Ritger A, McKenzie S, Libbrecht R, Kronauer DJC.Social regulation of insulin signaling and the evolution of eusociality in ants. Science 2018;361:398–402. doi: https://doi.org/10.1126/science.aar5723.
Korb J, Meusemann K, Aumer D, Bernadou A, Elsner D, Feldmeyer et al. Comparative transcriptomic analysis of the mechanisms underpinning ageing and fecundity in social insects. Phil Trans R Soc B. 2021;376:20190728. doi:https://doi.org/10.1098/rstb.2019.0728.
Corona M, Velarde RA, Remolina S, Moran-Lauter A, Wang Y, Hughes KA, et al. Vitellogenin, juvenile hormone, insulin signaling, and queen honey bee longevity. Proc Natl Acad Sci U S A. 2007;104:7128–33. doi:https://doi.org/10.1073/pnas.0701909104.
Rascon B, Mutti NS, Tolfsen C, Amdam G V. Honey bee life history plasticity: Development, behavior, and aging. In: Flatt T, Heyland A, editors. Mechanisms of Life History Evolution: The Genetics and Physiology of Life History Traits and Trade-Offs. New York: Oxford Univ Press; 2011. p. 253–66.
Rodrigues MA, Flatt T. Endocrine uncoupling of the trade-off between reproduction and somatic maintenance in eusocial insects. Curr Opin Insect Sci. 2016;16:1–8.
Weitekamp CA, Libbrecht R, Keller L. Genetics and evolution of social behavior in insects. Annu Rev Genet. 2017;51:219–39. https://doi.org/10.1146/annurev-genet-120116-024515.
Maekawa K, Ishitani K, Gotoh H, CornetteR, Miura T. Juvenile hormone titre and vitellogenin gene expression related to ovarian development in primary reproductives compared with nymphs and nymphoid reproductives of the termite Reticulitermes speratus. Physiol Entomol. 2010;35:52–8. doi:https://doi.org/10.1111/j.1365-3032.2009.00711.x.
Korb J. Juvenile hormone: a central regulator of termite caste polyphenism. In: Zayed A, Kent C, editors. Genomics, Physiology and Behaviour of Social Insects. Academic Press; 2015. p. 131 – 61. doi:https://doi.org/10.1016/bs.aiip.2014.12.004.
Lin S, Werle J, Korb J Transcriptomic analyses of the termite, Cryptotermes secundus, reveal a gene network underlying a long lifespan and high fecundity. Comm Biol. accepted
Tasaki E, Mitaka Y, Nozaki T, Kobayashi K, Matsuura K, Iuchi Y. High expression of the breast cancer susceptibility gene BRCA1 in long-lived termite kings. Aging (Albany NY). 2018;10:2668–83. doi:https://doi.org/10.18632/aging.101578.
Wu J, Lu L-Y, Yu X. The role of BRCA1 in DNA damage response. Protein Cell. 2010;1:117–23. doi:https://doi.org/10.1007/s13238-010-0010-5.
Tasaki E, Kobayashi K, Matsuura K, Iuchi Y. An efficient antioxidant system in a long-lived termite queen. PLoS One. 2017;12:e0167412. doi:https://doi.org/10.1371/journal.pone.0167412.
Tasaki E, Kobayashi K, Matsuura K, Iuchi Y. Long-lived termite queens exhibit high Cu/Zn-superoxide dismutase activity. Oxid Med Cell Longev. 2018;2018:5127251. doi:https://doi.org/10.1155/2018/5127251.
Münch D, Amdam G V, Wolschin F. Ageing in a eusocial insect: molecular and physiological characteristics of life span plasticity in the honey bee. Funct Ecol. 2008;22:407–21. doi:https://doi.org/10.1111/j.1365-2435.2008.01419.x.
Lucas ER, Keller L. Ageing and somatic maintenance in social insects. Curr Opin Insect Sci. 2014;5:31–6. doi:https://doi.org/10.1016/j.cois.2014.09.009.
de Verges J, Nehring V. A critical look at proximate causes of social insect senescence: damage accumulation or hyperfunction? Curr Opin Insect Sci. 2016;16:69–75. doi:https://doi.org/10.1016/j.cois.2016.05.003.
Kramer BH, Nehring V, Buttstedt A, Heinze J, Korb J, Libbrecht R et al. Oxidative stress and senescence in social insects: a significant but inconsistent link. Phil Trans R Soc B. 2021;76:20190732. doi:https://doi.org/10.1098/rstb.2019.0732.
Monroy Kuhn JM, Meusemann K, Korb J. Long live the queen, the king and the commoner? Transcript expression differences between old and young in the termite Cryptotermes secundus. PLoS One. 2019;14:e0210371. doi:https://doi.org/10.1371/journal.pone.0210371.
Korb J, Lenz M. Reproductive decision-making in the termite, Cryptotermes secundus (Kalotermitidae), under variable food conditions. Behav Ecol. 2004;15:390–5. doi: https://doi.org/10.1093/beheco/arh033.
Korb J, Katrantzis S. Influence of environmental conditions on the expression of the sexual dispersal phenotype in a lower termite: implications for the evolution of workers in termites. Evol Dev. 2004;6:342–52. doi: https://doi.org/10.1111/j.1525-142X.2004.04042.x.
Kim K, Lee YS, Carthew RW. Conversion of pre-RISC to holo-RISC by Ago2 during assembly of RNAi complexes. RNA. 2007;13:22–9. doi:https://doi.org/10.1261/rna.283207.
Sabin LR, Zhou R, Gruber JJ, Lukinova N, Bambina S, Berman A, et al. Ars2 regulates both miRNA- and siRNA- dependent silencing and suppresses RNA virus infection in Drosophila. Cell. 2009;138:340–51. doi: https://doi.org/10.1016/j.cell.2009.04.045.
Ivanov DK, Escott-Price V, Ziehm M, Magwire MM, Mackay TFC, Partridge L, et al. Longevity GWAS Using the Drosophila Genetic Reference Panel. J Gerontol A Biol Sci Med Sci. 2015;70:1470–8. doi: https://doi.org/10.1093/gerona/glv047.
Cooper, GM. Protein Synthesis, Processing, and Regulation (Chap. 7). In: The Cell: A Molecular Approach. 2nd. Sunderland: Sinauer Associates; 2000.
Wong E, Cuervo AM. Integration of clearance mechanisms: the proteasome and autophagy. Cold Spring Harb Perspect Biol. 2010;2:a006734. doi: https://doi.org/10.1101/cshperspect.a006734.
Putiri EL, Robertson KD. Epigenetic mechanisms and genome stability. Clin Epigenetics. 2010;2:299–314. doi:https://doi.org/10.1007/s13148-010-0017-z.
Gozani O, Shi Y. In: Workman JL, Abmayr SM, editors. Histone Methylation in Chromatin Signaling BT - Fundamentals of Chromatin. New York: Springer; 2014. p. 213–56. https://doi.org/10.1007/978-1-4614-8624-4_5.
Steunou A-L, Rossetto D, Côté J. Regulating Chromatin by Histone Acetylation BT - Fundamentals of Chromatin. In: Workman JL, Abmayr SM, editors. New York, NY: Springer New York; 2014. p. 147–212. doi:https://doi.org/10.1007/978-1-4614-8624-4_4.
Plongthongkum N, Kullawong N, Panyim S, Tirasophon W. Ire1 regulated XBP1 mRNA splicing is essential for the unfolded protein response (UPR) in Drosophila melanogaster. Biochem Biophys Res Commun. 2007;354:789–94.
Remondelli P, Renna M. The Endoplasmic Reticulum Unfolded Protein Response in Neurodegenerative Disorders and Its Potential Therapeutic Significance. Front Mol Neurosci. 2017;10:187. doi:https://doi.org/10.3389/fnmol.2017.00187.
Gorrini C, Baniasadi PS, Harris IS, Silvester J, Inoue S, Snow B, et al. BRCA1 interacts with Nrf2 to regulate antioxidant signaling and cell survival. J Exp Med. 2013;210:1529–44. doi: https://doi.org/10.1084/jem.20121337.
Korb J, Schneider K. Does kin structure explain the occurrence of workers in the lower termite Cryptotermes secundus? Evol Ecol. 2007;21:817–28.
Lienig J, Bruemmer H. Fundamentals of Electronic Systems Design. Berlin: Springer International Publishing; 2017. https://doi.org/10.1007/978-3-319-55840-0.
Piatek MJ, Werner A. Endogenous siRNAs: regulators of internal affairs. Biochem Soc Trans. 2014;42:1174–9. doi:https://doi.org/10.1042/BST20140068.
Dang W, Steffen KK, Perry R, Dorsey JA, Johnson FB, Shilatifard A, et al. Histone H4 lysine 16 acetylation regulates cellular lifespan. Nature. 2009;459:802–7. doi: https://doi.org/10.1038/nature08085.
Wood JG, Hillenmeyer S, Lawrence C, Chang C, Hosier S, Lightfoot W, et al. Chromatin remodeling in the aging genome of Drosophila. Aging Cell. 2010;9:971–8. doi: https://doi.org/10.1111/j.1474-9726.2010.00624.x.
Larson K, Yan S-J, Tsurumi A, Liu J, Zhou J, Gaur K, et al. Heterochromatin formation promotes longevity and represses ribosomal RNA synthesis. PLoS Genet. 2012;8:e1002473. doi: https://doi.org/10.1371/journal.pgen.1002473.
Wood JG, Helfand SL. Chromatin structure and transposable elements in organismal aging. Front Genet. 2013;4:274. doi: https://doi.org/10.3389/fgene.2013.00274.
Chen H, Zheng X, Xiao D, Zheng Y. Age-associated de‐repression of retrotransposons in the Drosophila fat body, its potential cause and consequence. Aging Cell. 2016;15:542–52. doi:https://doi.org/10.1111/acel.12465.
Wood JG, Jones BC, Jiang N, Chang C, Hosier S, Wickremesinghe P, et al. Chromatin-modifying genetic interventions suppress age-associated transposable element activation and extend life span in Drosophila. Proc Natl Acad Sci U S A. 2016;113:11277–82 doi: https://doi.org/10.1073/pnas.1604621113.
Liu Y, Ye X, Jiang F, Liang C, Chen D, Peng J, et al. C3PO, an endoribonuclease that promotes RNAi by facilitating RISC activation. Science. 2009;325:750–3. doi: https://doi.org/10.1126/science.1176325.
Kandasamy SK, Zhu L, Fukunaga R. The C-terminal dsRNA-binding domain of Drosophila Dicer-2 is crucial for efficient and high-fidelity production of siRNA and loading of siRNA to Argonaute2. RNA. 2017;23:1139–53. doi: https://doi.org/10.1261/rna.059915.116.
Tsuboyama K, Tadakuma H, Tomari Y. Conformational activation of Argonaute by distinct yet coordinated actions of the Hsp70 and Hsp90 chaperone systems. Mol Cell. 2018;70:722–729.e4. doi: https://doi.org/10.1016/j.molcel.2018.04.010.
Liu Y, Beyer A, Aebersold R. On the dependency of cellular protein levels on mRNA abundance. Cell. 2016;165:535–50. doi: https://doi.org/10.1016/j.cell.2016.03.014.
Cao JQ, Tong WS, Yu HY, Tobe SS, Bendena WG, Hui JHL. Chapter Three - The Role of MicroRNAs in Drosophila Regulation of Insulin-Like Peptides and Ecdysteroid Signalling: Where Are We Now? In: Verlinden HBT-A in IP, editor. Insect Epigenetics. Academic Press; 2017; p. 55–85. doi: https://doi.org/10.1016/bs.aiip.2017.02.002.
Laplante M, Sabatini DM. mTOR signaling at a glance. J Cell Sci. 2009;122 Pt 20:3589–94. doi: https://doi.org/10.1242/jcs.051011.
Zhao J, Zhai B, Gygi SP, Goldberg AL. mTOR inhibition activates overall protein degradation by the ubiquitin proteasome system as well as by autophagy. Proc Natl Acad Sci U S A. 2015;112:15790–7. doi: https://doi.org/10.1073/pnas.1521919112.
Ichimura Y, Waguri S, Sou Y-S, Kageyama S, Hasegawa J, Ishimura R, et al. Phosphorylation of p62 activates the Keap1-Nrf2 pathway during selective autophagy. Mol Cell. 2013;51:618–31. doi: https://doi.org/10.1016/j.molcel.2013.08.003.
Aramburu J, Ortells MC, Tejedor S, Buxade M, Lopez-Rodriguez C. Transcriptional regulation of the stress response by mTOR. Sci Signal. 2014;7:re2. doi: https://doi.org/10.1126/scisignal.2005326.
Tsakiri EN, Sykiotis GP, Papassideri IS, Terpos E, Dimopoulos MA, Gorgoulis VG, et al. Proteasome dysfunction in Drosophila signals to an Nrf2-dependent regulatory circuit aiming to restore proteostasis and prevent premature aging. Aging Cell. 2013;12:802–13. doi: https://doi.org/10.1111/acel.12111.
Sykiotis GP, Bohmann D. Keap1/Nrf2 signaling regulates oxidative stress tolerance and lifespan in Drosophila. Dev Cell. 2008;14:76–85. doi: https://doi.org/10.1016/j.devcel.2007.12.002.
Tullet JMA, Hertweck M, An JH, Baker J, Hwang JY, Liu S, et al. Direct inhibition of the longevity-promoting factor SKN-1 by insulin-like signaling in C. elegans. Cell. 2008;132:1025–38. doi: https://doi.org/10.1016/j.cell.2008.01.030.
Tolfsen CC, Baker N, Kreibisch C, Amdam GV. Flight restriction prevents associative learning deficits but no changes in brain protein-adduct formation during honeybee ageing. J Exp Biol. 2011;214:1322–32.
Parker JD, Parker KM, Sohal BH, Sohal RS, Keller L. Decreased expression of Cu-Zn superoxide dismutase 1 in ants with extreme lifespan. Proc Natl Acad Sci U S A. 2004;101:3486–9. doi:https://doi.org/10.1073/pnas.0400222101.
IBM Corp. Released 2015 IBM SPSS Statistics for Windows, Version 23.0. 2015. Armonk: IBM Corp.
Hosmer DW, Lemeshow S, May S. Applied Survival Analysis 2nd. 2008. Hoboken: Wiley.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. http://www.bioinformatics.babraham.ac.uk/projects/fastqc.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotech. 2016;34:525–7. doi:https://doi.org/10.1038/nbt.3519.
Harrison MC, Jongepier E, Robertson HM, Arning N, Bitard-Feildel T, Chao H, et al. Hemimetabolous genomes reveal molecular basis of termite eusociality. Nat Ecol Evol. 2018;2:557–66. doi: https://doi.org/10.1038/s41559-017-0459-1.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550. doi: https://doi.org/10.1186/s13059-014-0550-8.
R Core Team. 2018 R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna). https://www.R-project.org/.
Nueda MJ, Martorell-Marugan J, Marti C, Tarazona S, Conesa A. Identification and visualization of differential isoform expression in RNA-seq time series. Bioinformatics. 2018;34:524–6. doi: https://doi.org/10.1093/bioinformatics/btx578.
Nueda MJ, Tarazona S, Conesa A. Next maSigPro: updating maSigPro bioconductor package for RNA-seq time series. Bioinformatics. 2014;30:2598–602. doi:https://doi.org/10.1093/bioinformatics/btu333.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. doi:https://doi.org/10.1186/1471-2105-9-559.
Zhang B, Horvath S. A general framework for weighted gene coexpression network analysis. Stat Appl Genet Mol Biol. 2005;4:17. doi:https://doi.org/10.2202/1544-6115.1128.
Horvath S, Dong J. Geometric interpretation of gene coexpression network analysis. PLoS Comput Biol. 2008;4:e1000117. doi:https://doi.org/10.1371/journal.pcbi.1000117.
Yip AM, Horvath S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinformatics. 2007;8:22. doi:https://doi.org/10.1186/1471-2105-8-22.
Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics. 2008;24:719–20. doi: https://doi.org/10.1093/bioinformatics/btm563.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504. doi: https://doi.org/10.1101/gr.1239303.
Kucera M, Isserlin R, Arkhangorodsky A, Bader GD. AutoAnnotate: A Cytoscape app for summarizing networks with semantic annotations. F1000Research. 2016;5:1717. doi:https://doi.org/10.12688/f1000research.9090.1.
Enright AJ, Van Dongen S, Ouzounis CA. An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res. 2002;30:1575–84.
Maere S, Heymans K, Kuiper M. BiNGO: a Cytoscape plugin to assess overrepresentation of Gene Ontology categories in Biological Networks. Bioinformatics. 2005;21:3448–9. doi: https://doi.org/10.1093/bioinformatics/bti551.
Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Ser B. 1995;57:289–300. http://www.jstor.org/stable/2346101.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44:D279-85. doi: https://doi.org/10.1093/nar/gkv1344.
Eddy SR. Accelerated Profile HMM Searches. PLoS Comput Biol. 2011;7:e1002195. doi: https://doi.org/10.1371/journal.pcbi.1002195.
Jones P, Binns D, Chang H-Y, Fraser M, Li W, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30:1236–40. doi: https://doi.org/10.1093/bioinformatics/btu031.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10. doi: https://doi.org/10.1016/S0022-2836(05)80360-2.
Terrapon N, Li C, Robertson HM, Ji L, Meng X, Booth W, et al. Molecular traces of alternative social organization in a termite genome. Nat Commun. 2014;5:3636. doi: https://doi.org/10.1038/ncomms4636.
Hubley R, Finn RD, Clements J, Eddy SR, Jones TA, Bao W, et al. The Dfam database of repetitive DNA families. Nucleic Acids Res. 2016;44:D81-9. doi: https://doi.org/10.1093/nar/gkv1272.
Wheeler TJ, Eddy SR. nhmmer: DNA homology search with profile HMMs. Bioinformatics. 2013;29:2487–9. doi: https://doi.org/10.1093/bioinformatics/btt403.
Zdobnov EM, Tegenfeldt F, Kuznetsov D, Waterhouse RM, Simao FA, Ioannidis P, et al. OrthoDB v9.1: cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs. Nucleic Acids Res. 2017;45:D744-9. doi: https://doi.org/10.1093/nar/gkw1119.
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80. doi:https://doi.org/10.1093/molbev/mst010-
Misof B, Misof K. A Monte Carlo approach successfully identifies randomness in multiple sequence alignments: a more objective means of data exclusion. Syst Biol. 2009;58:21–34. doi: https://doi.org/10.1093/sysbio/syp006.
Kück P, Meusemann K, Dambach J, Thormann B, von Reumont BM, Wägele JW et al. Parametric and non-parametric masking of randomness in sequence alignments can be improved and leads to better resolved trees. Front Zool. 2010;7:10. doi: https://doi.org/10.1186/1742-9994-7-10.
Nguyen L-T, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74. doi: https://doi.org/10.1093/molbev/msu300.
Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9. doi: https://doi.org/10.1038/nmeth.4285.
Le SQ, Dang CC, Gascuel O. Modeling protein evolution with several amino acid replacement matrices depending on site rates. Mol Biol Evol. 2012;29:2921–36. doi: https://doi.org/10.1093/molbev/mss112
Gouy M, Guindon S, Gascuel O. SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010;27:221–4. doi:https://doi.org/10.1093/molbev/msp259.
We thank Florentine Schaub for assistance in the field and wet lab, Daniela Schnaiter for termite colony maintenance, Charles Darwin University (Australia), and especially S. Garnett and the Horticulture and Aquaculture team, provided logistic support to collect C. secundus. The Parks and Wildlife Commission, Northern Territory, the Department of the Environment, Water, Heritage and the Arts gave permission to collect. Permit numbers 2002 until 2016, export permit numbers: PWS2001_1508, PWS2003_39852, PWS2004_5769, PWS2007_4154, PWS2010_6997, PWS2014_001342, PWS2016_001559; collection permit numbers: 15656, 18310, 26851, 30073, 36401, 51402, 59044. The study was conducted according to the Nagoya protocol. KM thanks Ondrej Hlinka and the CSIRO IM&T HPC Cluster team.
This research was supported by the Deutsche Forschungsgemeinschaft by grants to JK (DFG; KO1895/16 − 1; KO1895/20 − 1), one within the Research Unit FOR2281. Open Access funding enabled and organized by Projekt DEAL.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Schematic workflow. Details are explained in the Methods. Figure S2. WGCNA modules of co-expressed transcripts that negatively correlate with age. Modules are named after colors by WGCNA. Eigengenes of all these modules showed a negative correlation with age while * indicates that the age-correlation was significant. Modules marked with $ were found in the young subnetwork. Figure S3. WGCNA modules of co-expressed transcripts that positively correlate with age. Modules are named after colors by WGCNA. Eigengenes of all these modules showed a positive correlation with age; * indicates that the age-correlation was significant.Additionally, modules marked with $ were found in the old subnetwork. Figure S4. GO enrichment for the WGCNAmodule ‘tan’ that positively correlated with age and which had many old-age DETs. Details are shown for BP (Biological Process), which revealed an enrichment of transcripts for ribosomal and tRNA related functions. Figure S5. Young and old transcript subnetworks corresponding to the second level neighborhood of Iso-MaSigPro DETs. Age-related DETs were located in the WGCNA co-expression network and these DETs and their one- and two-step neighbors (i.e., the ‘second level neighborhood’) were then extracted from the co-expression network to provide the shown networks. Figure S6. BiNGO GO enrichment(Biological Process) for the young subnetwork. No terms were significantly enriched after correcting for multiple testing (FDR). Figure S7. BiNGO GO enrichment (Biological Process) for the old subnetwork. Colored nodes are GO terms that were significantly enriched after correcting for multiple testing (FDR). FigureS8. BiNGO GO enrichment of (Molecular Function; Cellular Component) for the old subnetwork. Colored nodes are GO terms that were significantly enriched after correcting for multiple testing (FDR). Figure S9. BiNGO GO enrichment (Biological Process; MolecularFunction; Cellular Component) for the ‘green’ WGCNA module, which is part of the old subnetwork. Colored nodes are GO terms that were significantly enriched after correcting for multiple testing (FDR). Figure S10. BiNGO GO enrichment (Biological Process, Molecular Function; Cellular Component) for the ‘paleturquoise’ WGCNA module, which is part of the old subnetwork. Nodes in color are GO terms significantly enriched after correcting for multiple testing (FDR).
Data for survival analysis. Table S2. Sample information of samples included in this study. Table S3. Differentially expressed transcripts of the Iso-MaSigPro analysis, separately for cluster 1-6. Table S4. Transcripts in the young subnetwork (SNW). Transcripts in the young subnetwork (SNW) classified into major categories; TE-related transcripts of the young subnetwork (SNW). Table S5. Transcripts in the old subnetwork (SNW). Transcripts in the old subnetwork (SNW) classified into major categories; TE-related transcripts of the old subnetwork (SNW).
WGCNA module-age associations. Listed are the eigengenes and the correlation coefficients with respect to age (p-values in parenthesis).
GO terms for enriched differentially expressed transcripts (DETs) included in BiNGO module green; BiNGO module paleturquoise; BiNGO module tan; BiNGO module Thistle2; BiNGO module snow;BiNGO module cyan; BiNGO module deeppink1; BiNGO module navajowhite; BiNGO module blue2_NS; BiNGO module cornflowerblue_NS; BiNGO module pink3_NS; BiNGO module steelblue4_NS.
About this article
Cite this article
Monroy Kuhn, J.M., Meusemann, K. & Korb, J. Disentangling the aging gene expression network of termite queens. BMC Genomics 22, 339 (2021). https://doi.org/10.1186/s12864-021-07649-4
- Social insects
- Weighted gene co‐expression networks
- Time series