Longevity and transposon defense, the case of termite reproductives.

Social insects are promising new models in aging research. Within single colonies, longevity differences of several magnitudes exist that can be found elsewhere only between different species. Reproducing queens (and, in termites, also kings) can live for several decades, whereas sterile workers often have a lifespan of a few weeks only. We studied aging in the wild in a highly social insect, the termite Macrotermes bellicosus, which has one of the most pronounced longevity differences between reproductives and workers. We show that gene-expression patterns differed little between young and old reproductives, implying negligible aging. By contrast, old major workers had many genes up-regulated that are related to transposable elements (TEs), which can cause aging. Strikingly, genes from the PIWI-interacting RNA (piRNA) pathway, which are generally known to silence TEs in the germline of multicellular animals, were down-regulated only in old major workers but not in reproductives. Continued up-regulation of the piRNA defense commonly found in the germline of animals can explain the long life of termite reproductives, implying somatic cooption of germline defense during social evolution. This presents a striking germline/soma analogy as envisioned by the superorganism concept: the reproductives and workers of a colony reflect the germline and soma of multicellular animals, respectively. Our results provide support for the disposable soma theory of aging.

All animals age. However, some have a lifespan of a few weeks only (e.g., the fruit fly Drosophila melanogaster, the nematode Caenorhabditis elegans, the water flea Daphnia longispina), whereas others live for centuries (e.g., the clam Arctica islandica) (1, 2). Social insects like honey bees, ants, and termites have lifespan differences within single colonies that can differ by two orders of magnitude (3–5). Reproducing queens (and, in termites, also kings) can live for 20 y whereas the nonreproducing workers live only a few weeks to months. This makes social insects especially promising new models, as the within-species variation in rate of senescence only exists at the range of between-species variation outside social insects. Additionally, as social insect colonies are generally composed of families, all colony members share the same genetic background and differences in longevity are caused by differences in gene expression. Studying the molecular underpinning of aging in social insects is an emerging field with groundbreaking research in the honey bee and some ants (6–14). For termites, which evolved sociality independently from social Hymenoptera (wasps, bees, and ants), we are aware of only a single study (15).

From the nematode C. elegans to humans, aging in solitary organisms has been linked to an increased activity of transposable elements (TEs) (16–21). Active TEs can “jump” within the genome, thereby disturbing the regulation and expression of other genes, for instance, by transposing into another gene or regulatory region. The germline is protected against uncontrolled TE activity by the PIWI-interacting RNA (piRNA) pathway (22–24). Proteins from the PIWI group (Aubergine, Argonaute 3, Piwi), which are part of the Argonaute superfamily, interact with piRNAs and a core of other proteins to prevent activity of TEs (24).

To gain insights into the mechanisms underlying caste-specific aging in social insects, and especially in termites, we compared gene-expression profiles by using de novo-sequenced transcriptomes of heads from old and young individuals from different castes of the fungus-growing termite Macrotermes bellicosus (Smeathman, 1781; Fig. 1 and Dataset S1). This species has one of the highest longevity skews (80–120 fold) of all social insects, with queens and kings living as long as 20 y whereas workers have a lifespan of only a few months. Whereas workers are completely sterile, Macrotermes queens are specialized egg-layers that produce ca. 20,000 eggs per day (25). Given their long lifespan, they are probably the most fecund terrestrial animals. M. bellicosus has two worker castes, majors and minors, which differ in behavior, sex, and development. Minor workers are mainly involved in “indoor” tasks, such as brood care and constructing mounds; they are females and can develop into soldiers, thereby extending their lifespan (26, 27). Major workers have a higher random extrinsic mortality already early in life, as they generally do the more dangerous tasks like foraging for food; they are males and present developmental endpoints (26, 27).

Fig. 1.

M. bellicosus and its gene-expression pattern between castes and age classes. (A) An M. bellicosus mound with royal chamber harboring the queen (large whitish individual), the king (dark brown individual), and minor and major workers (smaller individuals). (B) Results of a PCA including all transcriptomes. Castes are separated along the first component (PC1), which explains 53% of the variance. Age classes of workers are separated along the second component (PC2), which accounts for 18% of the variance. By contrast, age classes of kings and queens cluster together. (C) Heat map of the DEGs for identified TEs. The height of the rows reflects the number of different TEs (shown as numbers). The color legend reflects the log10 of respective numbers (gray in the old kings, as they showed no significantly up-regulated genes). Termite symbols with a crown indicate reproductives: ♀, queens; ♂, kings; small and large termite symbols without a crown are minor and major workers, respectively. Old individuals are shown with a walking stick, in gray, with open symbols. Young individuals are shown in green with closed symbols.

To uncover molecular underpinnings of longevity differences, we analyzed transcriptomes of old and young queens, kings, and major and minor workers sampled from field colonies of a population that we have studied for more than 20 y. This gives us the unique opportunity not only to obtain and age long-lived reproductives from the wild but also to determine life expectancy and sample young and old individuals accordingly. Thus, it was possible to compare reproductives and nonreproductives from the wild that were at the maximum range of their natural lifespan (). Although Macrotermes reproductives can live for as long as 20 y in the laboratory, in the study population, the median lifetime expectancy after reaching the epigeal mound phase is 6 y (Fig. S1). This is because of the common occurrence of army ants that kill complete colonies and function as the major cause of random extrinsic mortality after colonies are established (28). Accordingly, old reproductives in this study had an age of at least 6 y, and young reproductives of 1 y, after mound emergence. The age differences between young and old major and minor workers were in the range of months ().

Results

Age-Related Changes in Gene Expression in Different Castes.

By using a total of 24 transcriptomes, a principal component analysis (PCA) showed that gene expression clustered by caste and age ( and Fig. 1). The first component separated castes, with kings being placed close to queens. Both worker castes formed a large cluster. The second component separated young and old workers, whereas all reproductives grouped together. No axis separated age classes of reproductives, implying that there are few age-related differences in gene expression in kings and queens.

This was supported when analyzing significantly differentially expressed genes (DEGs) at the single-gene level. We focused on the differential gene expression between young and old age classes within each caste (detailed in ). Of 13,959 expressed genes analyzed (Datasets S2, S3, and S11), only 26 genes were differentially expressed between old and young queens (18 up-regulated in old and 8 up-regulated in young). Kings had only two DEGs, which were up-regulated in young kings. This contrasts with the high number of DEGs in workers: 67 in minors (17 up-regulated in old and 50 up-regulated in young) and more than 5,000 in majors (3,125 up-regulated in old and 2,491 up-regulated in young; Dataset S3). These differences of aging in terms of DEGs between reproductives and workers are especially striking given that old and young workers differed in age by a few months only, whereas the age differences in reproductives was at least 5 y. Independent of age, queens had the most distinctive gene-expression pattern, with 717 genes being queen-specifically overexpressed compared with all other castes, whereas kings and minor and major workers had 50, 88, and 179, respectively (Fig. S2).

Age-Related Changes in Gene Expression and TE Activity.

Analyzing the function of DEGs by using D. melanogaster homologs and orthologs in Database for Annotation, Visualization and Integrated Discovery (DAVID) did not reveal any significant enrichments in old or young individuals for queens, kings, and minor workers (complete functional annotation is provided in Dataset S4). However, young major workers were characterized by an overexpression of genes related to metabolism compared with old major workers (Fig. S3, Table S1, and Dataset S5).

In old major workers, a striking 356 of 2,387 annotated (with Macrotermes natalensis) up-regulated genes (14.9%) were TEs. This significantly differed compared with young major workers, for which only 13 of 2,096 annotated up-regulated genes (0.6%) were TEs (contingency analysis, χ2 = 301.87, P < 0.001). Note that the different number of genes up-regulated in old or young major workers compared with the DEGs mentioned earlier is because, here, we considered only annotated genes. Up-regulated TEs comprised DNA transposons as well as retrotransposons (Fig. 1 and Dataset S6). We identified active TEs of the types Line/L1, Mariner Mos 1 (>50%), BovB (15%), PiggyBac (10%), Pogo (8%), and others (Fig. 1). Not all of the active TEs contained functional protein domains. Hence, some are likely nonautonomous TEs (29).

High TE expression is specific for major workers, and it differed significantly between major workers compared with the other castes (contingency analysis, majors vs. all other castes combined, χ2 = 6.85, P = 0.008). We found no TEs to be differentially expressed between age classes in queens and minor workers (Fig. 1). Between old and young kings, only in the young kings were genes up-regulated, and none of them were TEs.

These differential-expression (DE) analyses were done with DESeq2. A corresponding pattern was obtained when using edgeR (, , and Dataset S7) (30, 31). Furthermore, a quantitative real-time PCR (qRT-PCR) experiment for two TEs, Mariner Mos 1 and BovB, with additional samples from several colonies, confirmed the results based on transcriptome data (Datasets S8 and S11). Both TEs were significantly up-regulated in old major workers compared with young major workers based on samples from a total of seven colonies (Mariner Mos 1, old, 0.55 ± 0.06 SEM, n = 6; young, 0.37 ± 0.03 SEM, n = 5; t test for unequal variances, t = −2.61, P = 0.028; BovB, old, 0.62 ± 0.09 SEM, n = 6; young, 0.31 ± 0.01 SEM, n = 5; t test for equal variances, t = −3.42, P = 0.018; note, we did not find old and young major workers for each colony; one sample was analyzed per colony). Hence, our results show that old major workers were characterized by a high activity of TEs. As high TE activity has been linked to senescence in D. melanogaster (19–21, 32–34) and other animals (17, 18), the short intrinsic lifespan of major workers may be a result of their high TE activity.

This leads to the following questions. Which mechanisms might prevent TE activity in reproductives and minor workers, and which mechanisms are switched on in young but off in old major workers?

Expression of TE-Silencing Pathways in Different Castes.

In D. melanogaster, TEs are silenced by piRNAs, which function in the germline (22–24, 35). To test whether similar TE-silencing mechanisms occur in M. bellicosus, and to what extent they are differentially expressed, we identified orthologs of key genes known to be involved in piRNA-mediated TE silencing (22, 24, 35), checked protein domains (Dataset S11), and analyzed their expression pattern (Fig. 2).

Fig. 2.

DEGs that are involved in TE silencing according to results for D. melanogaster. (A) Shown is a simplified version of the ping-pong amplification cycle [after Czech and Hannon (24)] of the piRNA pathway, which is involved in silencing TEs in the germline of animals. Ago3 with its sense piRNA recognizes and cleaves piRNA cluster transcripts, which are loaded into Aub1, probably assisted by a Tudor protein and Qin, which prevents Aub1/Aub1 pairings. Maturation of piRNAs involves several proteins, including Mino, Zuc, Papi, and Hen1, and results in mature antisense piRNA–Aub1 complexes. Methylation by Capsuleen (Csul) is essential for stable Ago3 and Aub proteins (36). These complexes can detect and slice transposon RNAs. The 3′ cleaved product is then loaded in Ago3 and, after modification, the cycle can restart. The expression of the following genes involved in the ping-pong cycle was down-regulated in old major workers compared with young ones (mean ± SEM): (B) aub1 (C) zuc, (D) qin, and (E) csul.

piRNA pathway.

piRNAs are a distinct class of small noncoding RNAs that form the piRNA-induced silencing complex (piRISC) in the germline of many animals (18–24, 35). The piRISC protects the integrity of the genome by silencing TEs. In contrast to siRNA or miRNA processing, it is Dicer-independent. The piRNA pathway can be separated into two parts: de novo biogenesis of piRNAs and the “ping-pong” amplification cycle. During biogenesis in somatic cells, primary piRNAs are produced from piRNA clusters (i.e., sites in the genome with high piRNA coverage) that serve as recognition and complementary binding sites for TEs in the ping-pong amplification cycle. The latter couples piRNA biogenesis to target silencing and amplifies primary piRNA in fly germ cells (Fig. 2). There is considerable overlap of genes between both parts of the piRNA pathway (Table S2). Based on the recent review by Czech and Hannon (24), we unambiguously identified nine genes from the piRNA pathway (Table S3). This includes three piRNA genes (ago3, aub1, and aub2) that are characterized by a Piwi and PAZ domain, typical for D. melanogaster PIWIs (Fig. S4; ; rationale of naming these genes is described in ).

Ping-pong amplification cycle.

Several genes involved in the ping-pong amplification cycle were significantly down-regulated in old compared with young major workers (Table S2). They include the key genes aub1 as well as capsuleen, qin, and zucchini (Fig. 2). These are genes representing major processes of the ping-pong cycle in D. melanogaster, such as piRNA processing or linking Aub and Ago3 to transfer transcripts (Fig. 2). In the edgeR analysis, aub1 was no longer significantly differentially expressed, but the remaining DEGs of the ping-pong cycle were significantly down-regulated in old major workers (Dataset S7). By contrast, none of the identified ping-pong cycle genes were differentially expressed between age classes in queens, kings, or minor workers (Table S2). This suggests that, specifically in old major workers, the ping-pong amplification cycle is less active, in line with an increased TE activity.

De novo biogenesis.

In contrast to the differential expression of genes from the ping-pong cycle, none of the genes specifically involved in the somatic de novo biogenesis of piRNAs known from D. melanogaster was differentially expressed, such as armi or putative gasz (Table S2). SoYb, which was significantly down-regulated in the old major workers, is part of the de novo biogenesis but is germline-specific (37). One shared gene, zucchini, belongs to the ping-pong cycle and the de novo biogenesis. The expression of these genes was confirmed in DESeq2 and edgeR. These results imply that the somatic, de novo biogenesis part of the piRNA pathway is not associated with TE activity in old major workers, but that the germline-specific ping-pong amplification cycle is.

Discussion

Our results show that several genes from the ping-pong cycle that, in other organisms, protect the germline against TE activity (and hence aging) are specifically down-regulated in old compared with young major workers. Their expression pattern is diametrically opposite to the activity of TEs. This supports the hypothesis that a defense mechanism against TEs known from the germline of flies and other animals functions in termites to protect the soma of reproductives and minor workers. Note that we analyzed head transcriptomes (Dataset S1); this introduces selectivity () but now allows us to exclude any germline influence. We do not rule out that other processes such as oxidative stress may play a role in aging, but we focused on the striking pattern of differentially expressed TEs in this study (the complete list is provided in Datasets S2 and S3).

Comparing Major and Minor Workers.

We found major differences between minor and major workers in age-related gene expression. Only major workers showed an increased TEs activity and a down-regulation of piRNA genes with age, whereas minor workers were more similar to the reproductives. Two factors, which may be interlinked, can account for this when applying evolutionary life history reasoning and regarding colony fitness as the relevant fitness level. The latter is reasonable because all workers in this species are sterile and cannot gain any direct fitness; colony fitness increases their indirect fitness. Hence, the workers of this species can be regarded as analogous to the soma of a multicellular organism, and we propose that aging theory (38, 39), such as the disposable soma theory (40, 41), applies equivalently. As a first factor, age-specific random extrinsic mortality differs between major and minor workers in M. bellicosus. Major workers do mainly outdoor tasks (i.e., foraging) for which random extrinsic mortality is extremely high in this species (42). By contrast, minor workers do indoor tasks when young and later do outdoor tasks (43). As extrinsic mortality is much higher outdoors than indoors, life history theory should analogously predict that selection against intrinsic aging (e.g., investment into physiological mechanisms that increase longevity such as TE defenses) should be lower in major compared with minor workers as the reproductive value of the major workers is lower. The reproductive value is measured/realized in these sterile castes via colony fitness. Second, major workers are developmental endpoints whereas minor workers can develop into soldiers. This additionally increases the reproductive value of minors compared with majors and hence selects more strongly against (intrinsic) aging in minors. Both factors give minor workers a higher reproductive value than major workers. This may explain the benefit of a higher continued investment in TE defense in minor workers. Life history models specifically tailored to social insects are required to test this explanation.

TE Activity, Aging, and the Germline.

There is substantial evidence, from C. elegans to humans, that senescence correlates with increasing TE activity (16–24, 44), and several studies suggest a causal link for D. melanogaster (32–34, 21). Uncontrolled TE activity is a threat to genome integrity as well as proper functioning of gene expression. Silencing TEs is especially important in the germline, as it “needs” more protection than disposable soma cells to sustain itself through generations. Accordingly, the piRNA pathway is mainly active in the germline, although recent evidence suggests that it can also be active in the brain and maybe some other tissues (44, 45). We can only speculate why this defense mechanism is largely restricted to the germline. Most likely, it is (energetically) costly, and a disposable soma is not worthwhile to be protected under all circumstances (40). Protection against TEs may present some of the germline-maintenance costs, which should receive more attention in aging research (41).

The Germline Analogy in Social Insects.

In highly social insects with sterile workers, the queens and kings can be regarded as equivalent to the germline of a colony, whereas the workers represent the disposable soma. Similar to organs in a multicellular organism, virtually all reproduction is channeled through long-lived queens (and, in termites, also kings), whereas the workers specialize on tasks like food provisioning or health care. This organism-like analogy is reflected in the superorganism concept. Our results provide support for this analogy and the concept of the disposable soma. Causal evidence has to be shown in further studies. The reproductives, and to some extent the “more valuable” minor workers, seem to be protected against aging by a mechanism that is mainly active in the germline of solitary animals.

Methods

Sampling, Age Determination, Laboratory Work, and Sequencing.

The M. bellicosus samples were collected in Kakpin, next to the Comoé National Park in Côte d’Ivoire (coordinates 8°39′N 3°46′W), where we have been studying this species for more than 20 y. We collected individuals of four castes (two age classes each): major and minor workers, queens and kings. The castes were identified by using the description in ref. 46. Of each of the worker castes, six individuals (three old and three young) were sampled. Sample size was constrained by the availability of colonies of known age and a low number of major workers of clearly distinct age (as detailed later). Of the sexual castes, five kings (two old, three young) and seven queens (three old, four young) were collected (Datasets S1 and S8).

We compared queens and kings from colonies that just had started to build a termite mound the previous year (i.e., young) with those that had mounds continuously inhabited for at least 6 y (i.e., old). Queens and kings are known to live as long as 20 y in the laboratory. The median longevity after emergence of an epigeal mound in our study site was estimated with Kaplan–Meier survival analysis for 36 mounds (Fig. S1) with the use SPSS 23 software (47).

For the workers, all individuals were collected in and around the queen chamber to control for a confounding effect of task on gene expression. Age was determined by checking mandible wear, a marker for age in termite workers (48, 49). Age differences between young and old reproductives were at least 5 y, whereas those for workers were only a few months.

Transcriptomes were generated from termite heads to avoid bacterial contamination from the gut. Additionally, by limiting ourselves to a single body part, we tried to avoid confounding factors/dilution of the signal. Therefore, the head was separated from all individuals. Altogether, 24 samples (one individual per sample; Dataset S1) were stored in Eppendorf tubes with RNAlater before the samples were transported to the laboratory for RNA extraction. Heads of all individuals were flushed with RNAlater (Qiagen) to ensure complete permeation of the head tissue. The samples were transported on ice to the field station and kept overnight at 4 °C to ensure full permeation. Then samples were kept frozen in a −20 °C refrigerator for no more than 2 wk. During transport, the samples were kept in an insulated container, cooled on ice, until arrival in the laboratory in Freiburg, Germany. There, samples were kept in a −80 °C freezer until further processing. Eleven additional samples were treated in the same manner for the qRT-PCR analysis; eight of these were collected in the Pendjari National Park (coordinates 11°19′ N 1°35′E; Dataset S8).

Total RNA was extracted from the head of single individuals according to a modified version of the peqGOLD TriFast protocol. Modifications were an incubation time of 10 min instead of 5 min at room temperature after adding 900 μL peqGOLD TriFast and inverting the tube. After adding chloroform and after incubation, centrifugation was modified to 10 min at 4 °C, 12,388 × g. After transfer of the aqueous phase into a new Eppendorf tube, we added another centrifugation step for 5 min at 4 °C (12,388 × g). The supernatant was then transferred into a new Eppendorf tube and 2 μL nuclease-free glycogen (5 mg/mL) was added per sample. After adding isopropanol, the mixture was incubated for 10 min at room temperature and not on ice. After washing the RNA pellet with 75% ethanol and vortexing, the sample was centrifuged for 5 min at 4 °C, 7,674 × g. This step was repeated three times.

After dissolving the pellet in nuclease-free water, samples were kept for at least 1 h at 4 °C. DNA digestion was performed by using DNase I Amplification Grade Kit (Sigma-Aldrich). The concentration and quality of the isolated RNA was checked with an Agilent Bioanalyzer (Agilent RNA 6000 Nano Kit) and sent to BGI in Hong Kong on dry ice.

Library preparation was done by BGI. Libraries for 22 samples showing slight degradation of RNA were constructed with the Ribo-Zero rRNA Removal Kit (human/mouse/rat; Illumina). Additionally, for 5 of these 22 samples, RnaseH treatment was performed. For two samples, the TruSeq RNA Kit was applied following the BGI sample preparation protocol. Details for each sample are given in Dataset S1.

Amplified libraries were sequenced on an Illumina HiSeq4000 platform with strategies of 100-bp paired-end reads, generating ∼4 Gb of raw data for each sample. After sequencing, index sequences from the machine reads were demultiplexed (sorted and removed) by a proprietary in-house tool at BGI.

Preprocessing of RNAseq Raw Reads.

Raw sequence reads from all samples provided by BGI were used for all further steps. After quality control with FastQC v. 0.11.4 (50), raw sequence reads were checked for a minimum length of 70 bp, and adapters (including two BGI in-house adapters) were trimmed using Trimmomatic v. 0.33 (51). Trimming and other subsequent steps were facilitated by GNU parallel v. 20141022 (52).

Gene Expression and Mapping Against the Genome of M. natalensis.

Because there is no sequenced M. bellicosus genome, we used the genome of M. natalensis (Haviland, 1898) v. 1.0 as a mapping backbone in our study (53). This species is closely related to our study species M. bellicosus. To test the suitability of M. natalensis as backbone (54), we downloaded mitochondrial cytochrome oxidase I (COI) sequences and sequence fragments for both Macrotermes species from GenBank (two sequence fragments of M. natalensis and six sequence fragments of M. bellicosus; Dataset S11). After aligning the sequences with the program Seaview (55) with the implemented clustal-omega (-clustalo) algorithm, we calculated with distmat (EMBOSS package v. 6.6.0.0) the distances between all COI sequences (all base positions considered, Kimura correction method) (56). The short evolutionary distances (<11) between COI sequences of M. natalensis and M. bellicosus suggest a short divergence time. This is in line with Brandl et al. (57), who determined the divergence time between the major clades of Macrotermes to be ∼8 My. Hence, M. natalensis can be used as a mapping backbone.

Trimmed and cleaned raw reads were mapped against the M. natalensis genome by using the program Tophat2 v. 2.1.0 with default settings (58). We created bowtie2-indices (59) based on the genome assembly (v. 1.0) and the corresponding GFF file as reference on which Tophat2 relies. The resulting BAM files were sorted by gene name by using SAMtools (v. 0.1.19–96b5f2294a) (60). HTSeq v. 0.6.0 was then used to count the mapped raw reads per gene (stranded = no, type = gene, mode = union, order = name) for each sample (61). A plain text table summarizing all samples was compiled (Dataset S11). Based on all valid counts identified by HTSeq, a background file with all expressed and mapped genes in M. bellicosus was built for later enrichment and pathway analyses.

Identification and Annotation of Clusters of Orthologous Groups.

We inferred a set of clusters of orthologous sequence groups (COGs) from the official gene sets, i.e., all protein-coding sequences, at amino-acid level of the fruit fly D. melanogaster (v. r6.11) and the termite M. natalensis (v. 1.2; Dataset S11). For both species, we kept only the longest isoform per gene and removed all sequences with internal stop codons and/or selenocysteines (U) by custom-made perl scripts. The latter may cause problems in downstream analyses. For inferring COGs, we used OrthoFinder v. 0.4.0, a graph-based approach that relies on an all-vs.-all BLAST and MCL clustering and corrects for gene length bias (62). We considered only 1:1 orthologs.

All transcripts of the M. natalensis genome assembly that had a positive match with the M. bellicosus raw reads after HTSeq were matched to all pathways, Gene Ontology (GO) terms, and protein domains in the Interpro database (63), including the databases Pfam, Gene3D, and PANTHER, with InterProScan (5.16–55.0, December 2015) (64) using default settings.

Gene-Expression Analysis.

DESeq2 v. 1.10.1 (30) was used in R v. 3.3.1 (65) to analyze differential gene expression between young and old individuals separately for each caste and between castes regardless of age. In addition, the obtained normalized read counts (Dataset S11) were used internally in DESeq2 to determine differential expression of all genes.

The count data were variance-transformed to perform a PCA and create a distance-clustering heat map. This allowed us to visualize the data and compare the effect of age and caste on gene expression. To test for pairwise expression differences, we used the generalized linear model with negative binominal distribution as implemented in DESeq2. We compared young and old age classes within each caste by combining caste and age factors into a group factor, as described in the DESeq2 vignette. Additionally, we also compared differential gene expression between castes (Dataset S11). False discovery rate (FDR)-adjusted P values were used to correct for multiple testing (66). Venn diagrams were produced to visualize the number of shared and unique genes between castes by using the online tool Venny 2 (67). We used the adjusted P values to determine DEGs of the piRNA pathway (Table S2 and Dataset S9). To corroborate our results, we additionally applied the algorithm implemented in edgeR to our data (31) (, , Fig. S5, and Dataset S7).

Functional Annotation and Enrichment.

A GO enrichment analysis (GO category, Biological Process) was done with the DAVID 6.8 Web tool (68). We supplemented our 1:1 ortholog set with homologs obtained by a BLASTP search between M. natalensis and D. melanogaster (BLAST v. 2.2.31+; National Center for Biotechnology Information) (69) with default settings and an e-value cutoff of 1e−5 of M. natalensis vs. D. melanogaster. We generated a tag cloud for significant results (P < 0.05 after FDR correction). Tag clouds were created in R 3.3.1 using the package tagcloud v. 0.6 (https://cran.r-project.org/package=tagcloud).

Transposable Elements.

We chose the annotation of the M. natalensis genome provided by Haofu Hu (Dataset S2) to identify TEs (51). Additionally, we checked the Pfam database (70) (Pfam A, release 30), the Pfam annotation from the InterPro (61), and the Dfam database (71) for the presence of TEs in our DEGs (Dataset S11).

We compared the number of differentially expressed TEs between young and old individuals within each caste by using contingency analyses. This was possible only for major workers, as differentially expressed TEs were missing between age classes in the other castes. Additionally, we tested whether the frequency of differentially expressed TEs differed between major workers and the other castes by using contingency analyses. All tests were two-tailed, and all analyses were performed with SPSS 23 (47).

RNAi Pathway Genes Involved in TE Silencing.

We identified genes related to the piRNA pathway and their relatives (Table S3 and Dataset S10): we obtained gene symbols and FlyBase IDs of D. melanogaster (e.g., “AGO1,” “FBgn0262739”; Dataset S10) from FlyBase (72). We used the FlyBase IDs to search in OrthoDB v. 9.1 (73) for COGs of the following species: D. melanogaster (DMEL), Apis mellifera (AMEL), Tribolium castaneum (TCAS), Zootermopsis nevadensis (ZNEV), and Blattella germanica (BGER). We downloaded amino-acid sequences for each COG of the aforementioned species. We aligned the sequences separately for each COG with MAFFT v7.294b, choosing the G-INS-i vsm algorithm to avoid overalignment of sequences [(option allowshift–unalignlevel-leavegappyregion (74)]. We used hmmbuild implemented in HMMER v3.1b2 to construct a hidden Markov model (HMM) from each multiple sequence alignment (75). The HMM was used to search with hmmsearch (HMMER) against the official gene set v. 1.2 of M. natalensis to identify candidate sequences for each COG. We searched all candidate sequences with hmmscan (HMMER) against the Pfam database (Pfam A, release 30; inclusion threshold 1e−2). The candidate sequences were searched reciprocally against the official gene set of D. melanogaster (v. r6.11) with BLASTP and an e-value cutoff of 1e−10. In cases for which the original sequence of D. melanogaster was not identical with the one after the reciprocal BLASTP search, we considered results as ambiguous.

For selected genes related to the piRNA pathway (Dataset S10), we inferred phylogenetic trees of gene groups. We included sequences of all species mentioned earlier (DMEL, AMEL, TCAS, BGER, ZNEV, Macrotermes) plus amino acid sequences from the recently sequenced genome of the termite Cryptotermes secundus (data provided by J.K.). In particular, we created sets of sequences that we previously had assigned to more than one gene with the HMM search and reciprocal BLAST search against D. melanogaster: the AGO1-3/Piwi/Aub set (AGO1, AGO2, AGO3, Piwi/Aubergine; Dataset S10). For the set of sequences, the amino acid selenocysteine (U) was replaced by “X.” We aligned each set of sequences with the G-INS-i vsm algorithm (option allowshift–unalignlevel-leavegappyregion) and subsequently checked each multiple sequence alignment for ambiguously aligned sequence regions with Aliscore v. 2 (76, 77), allowing the maximal number of pairwise comparisons, the e-option for gappy alignments, and otherwise default settings. We removed all positions suggested by Aliscore with Alicut v. 2.3 (https://www.zfmk.de/en/research/research-centres-and-groups/utilities; Dataset S11). We conducted 10 maximum-likelihood tree searches per sequence set with IQ-Tree v. 1.5.3 (78). We used completely random start trees and the automated estimation of the best model considering implemented nuclear models (with +G, +I, or +G+I, four rate categories, and the median approximation for +G site rates) plus the free rate model LG4X (79) with four rate categories. Statistical support was inferred from 1,000 nonparametric bootstrap replicates. We plotted bootstrap support on the maximum-likelihood tree with the best log-likelihood value. We visualized the unrooted trees with bootstrap support by using Seaview v. 4.5.4 (55) and graphically processed the trees with Inkscape v0.91 (www.inkscape.org/).

The sequence of Mnat_13309 in the genome of M. natalensis was short ( and Dataset S11). To obtain further information, we assembled the original raw reads from M. bellicosus into contigs with Trinity v2.0.6 (80) and confirmed the identity of the best matching contig to Mnat_13309 as ago3 using a reciprocal BLAST against D. melanogaster and a domain search against Pfam release 30. The contig was most similar to D. melanogaster Ago3 and contained Piwi and PAZ domains.