Genomic and transcriptomic analyses of the subterranean termite Reticulitermes speratus: Gene duplication facilitates social evolution

Genome sequencing of R. speratus was performed with genomic DNA isolated from female secondary reproductives (nymphoids) (Fig. 1B). R. speratus nymphoids are almost exclusively produced parthenogenetically by automixis with terminal fusion in primary queens, such that the genome should be homozygous at most loci (48), which provides an advantage in genome assembly. We generated a total of 86 Gb of Illumina HiSeq sequence data and performed de novo assembly into 5,817 scaffolds with an N50 of 1.97 Mb and total size of 881 Mb (SI Appendix, Table S1), covering 88% of the genome based on the genome size (1.0 Gb) estimated by flow cytometry (49). The assembled R. speratus genome shows high coverage of coding regions, capturing 99.2% (98.5% complete; 0.7% fragmented) of 1,367 Insecta benchmarking universal single-copy orthologs (50) (SI Appendix, Table S1). The R. speratus genome is rich in repetitive elements, which make up 40.4% of the genome. The high accumulation of repetitive elements may be relevant to the evolution of eusociality, as proposed for termites and snapping shrimp (22, 51). A total of 15,591 protein-coding genes were predicted by combining the reference-guided assembly of RNA-seq reads (36 libraries derived from different castes, sexes, and body parts; see Transcriptome Differentiation among Castes for details) and homology-based gene prediction followed by manual curation of gene families of interest (Fig. 1C). Whole-genome BS-seq revealed extensive gene body methylation of the R. speratus genome, accounting for 8.8% of methylated cytosines in the CG context (SI Appendix, Fig. S1). The genome-wide DNA methylation landscape was similar to that of the dampwood termite Z. nevadensis (12%) (24). These omics data and a genome browser are available at http://www.termite.nibb.info/retsp/.

We compared the R. speratus gene repertoire with those of 88 other arthropods, including the two termites Z. nevadensis and Macrotermes natalensis (20, 52). Ortholog analysis showed that 12,032 (82.9%) of the 15,591 genes in R. speratus genes were shared with other arthropods, while 1,773 were taxonomically restricted (TRGs) to Isoptera, among which 430 were shared with the other two termites and 1,343 were unique to R. speratus (Fig. 1C). Whole-genome comparisons with two sequenced termites, M. natalensis and Z. nevadensis, showed a high degree of synteny conservation (Fig. 1D). We identified 2,799 syntenic blocks (N50: 858.4 kb) shared with M. natalensis, covering 95.1% of the R. speratus genome, in which 560.4 Mb of nucleotides were aligned, while 3,650 syntenic blocks (N50: 591.1 kb) were shared with Z. nevadensis, covering 72.1% of the R. speratus genome, 116.7 Mb of which were aligned. Only a few cases of large genomic rearrangements were found between the termite genomes, at least at the contiguity level of the current assemblies, suggesting overall conservation of genome architecture in the lineage of termites over 135 million y (Fig. 1A). Interestingly, despite such high conservation of macrosynteny, interruptions or breaks in local synteny were observed and were often associated with tandem gene duplications. For example, when we examined regions containing large tandem gene duplications (>5-gene tandem duplications), the synteny between the R. speratus and M. natalensis genomes was found to be interrupted in 10 of 21 regions (examples shown in SI Appendix, Fig. S2).

Distinct castes arise from the same genome via a phenomenon called caste polyphenism, which is a distinctive hallmark of social insects (53, 54). To elucidate caste-biased gene expression to understand the mechanism underlying caste-specific phenotypes, we compared the transcriptomes of three castes (primary reproductives, workers, and soldiers) of R. speratus (Fig. 1B and SI Appendix, Table S2). We sequenced 36 RNA-seq libraries, representing three biological replicates of both sexes and two body parts (“head” and “thorax + abdomen”) for each of the three castes.

The results clearly showed that the termite castes were distinctively differentiated at the gene expression level. The multidimensional scaling plot depicted the three castes as clearly distinct transcriptomic clusters for both the head and thorax + abdomen transcriptomes (Fig. 2A). However, little sexual difference was detected within each caste, although reproductives showed substantial transcriptomic differences in thorax + abdomen samples between queens and kings, probably due to the differences in the reproductive organs (Fig. 2A). Using a generalized linear model (GLM) with caste and sex as explanatory variables, we identified 1,579 and 2,076 genes that were differentially expressed among castes in the head and thorax + abdomen samples, respectively, according to the criterion of a false discovery rate (FDR)-corrected P < 0.01, while we identified only 6 and 79 genes that were differentially expressed between the sexes in the head and thorax + abdomen samples, respectively, with the same criteria. We focused on the genes that were differentially expressed among castes (caste-DEGs) and further classified them into three categories of caste-biased genes (i.e., reproductive-, worker-, and soldier-biased genes), according to the criterion of >twofold higher expression than that in the other two castes (Fig. 2B and Dataset S1). These caste-biased genes should account for the specialized functions of each caste. Soldier samples exhibited the highest number of caste-specific genes, suggesting the highly specialized functions of the soldier caste (Fig. 2B). This was consistent with the finding of a previous RNA-seq analysis of the Eastern subterranean termite Reticulitermes flavipes indicating that a majority of the identified DEGs were soldier specific (55); 73 of the 93 identified DEGs were up- or down-regulated specifically in the soldier caste. In addition to these soldier-specific R. flavipes genes (e.g., troponin C and fatty acyl-CoA reductase), the caste-biased genes identified in our transcriptome analysis of R. speratus included genes previously reported to be caste-biased genes in other termites (56–59), such as vitellogenin (reproductives), geranylgeranyl pyrophosphate synthase (soldiers), and beta-glucosidase genes (probably associated with cellulase; workers). This consistency between the transcriptome analyses of different termite species indicates that the RNA-seq analysis conducted in this study is reliable and that the regulation and perhaps the functions of these caste-biased genes are conserved across the termite lineage.

The caste-DEGs were enriched for gene ontology terms related to a wide array of functions (SI Appendix, Table S3), such as terpenoid metabolism, hormone metabolism, chitin metabolism, hydrolase activity, oxidoreductase activity, lipid metabolism, signaling, eye pigmentation, serine-type peptidase, thermotaxis, heme binding, and lysozyme activity. Protein motifs enriched in the caste-DEGs were also identified, including cytochrome P450, lipocalin, lysozyme, glycosyl hydrolase family, TGF-beta, and trypsin motifs (SI Appendix, Table S4). The genes of these categories were further investigated and discussed in the context of termite biology related to 14 categories (lipocalins, GGPP synthases, cellulases, lysozyems, sex determination, epigenetics, chemosensory system, biogenic amines and neuropeptides, juvenile hormone, ecdysone, insulin signaling, wing formation, immunity, insecticide target, and detoxification) in Gene Duplication and Caste-Biased Gene Expression and SI Appendix.

Among the 1,773 TRGs that were restricted to Isoptera, the termite-shared TRGs showed strong enrichment for caste-DEG (Fisher's exact test, P < 1.0 × 10⁻⁷ for head samples and P < 1.0 × 10⁻¹⁰ for thorax + abdomen samples), while the TRGs found only in R. speratus (orphan genes) did not (P = 0.99 and P = 0.97, respectively).

To investigate the relationship between caste-biased gene expression and DNA methylation, we analyzed differential methylation levels among three castes (reproductives, workers, and soldiers). The BS-seq data showed that the global CpG methylation patterns were very similar among the castes (SI Appendix, Fig. S3 A and B), in contrast to the methylation pattern of Z. nevadensis, in which DNA methylation differs strongly between castes (winged adults versus final-instar larvae) and is strongly linked to caste-specific splicing (24). Instead, gene body DNA methylation in R. speratus seems to be important for the expression of housekeeping genes, as reported in the drywood termite C. secundus (21). Housekeeping genes exhibited a high degree of gene body methylation in all castes of R. speratus, while caste-biased genes showed a significantly lower level of DNA methylation (SI Appendix, Fig. S3 C and D).

Evolutionary novelties are often brought about by gene duplications (25, 26). Our ortholog analysis comparing the R. speratus gene repertoire with those of 88 other arthropods identified 1,396 multigene families that are duplicated in the R. speratus genome (SI Appendix, Table S5). Interestingly, compared to the genome as a whole, the set of caste-DEGs described in Transcriptome Differentiation among Castes was significantly enriched for genes in multigene families (X-squared = 218.62, df = 1, P < 2.2 × 10⁻¹⁶). We also calculated the tau score as a proxy of the caste specificity of gene expression for all genes and found that duplicated genes were significantly more caste specific than single-copy genes (P < 2.2 × 10⁻¹⁶, Wilcoxon rank sum test) in both transcriptome datasets (head and thorax + abdomen) (Fig. 2C). Additionally, gene set tests showed that sets of duplicated genes were differentially expressed in all pairwise comparisons between castes (Fig. 2D). These data highlight the important roles of duplicated genes for the genetic regulation of caste differentiation in R. speratus.

Caste-biased multigene families were associated with diverse functional categories, some of which were strongly related to caste-specific behaviors and tasks. Here, we highlight five families, namely, lipocalins (protein transporters for social communication and physiological signaling), cellulases (carbohydrate-active enzymes [CAZymes] or worker wood digestion), lysozymes (immune-related genes for social immunity), geranylgeranyl diphosphate (GGPP) synthases (metabolic enzymes for the production of soldier defensive chemicals), and a termite-specific gene family with unknown functions, as examples of multigene families relevant to termite sociality. Molecular evolution studies have shown that the redundancy caused by gene duplication may allow one paralog to acquire a new function (neofunctionalization) or divide the ancestral function among paralogs (subfunctionalization) (25, 26). We are particularly interested in the evolutionary impact of gene duplication on caste specialization through neo/subfunctionalization.

Numerous studies have shown that novel genes (e.g., TRGs) play important roles in the evolution of novel social phenotypes in hymenopteran social insects (9, 86, 87). We found that termite-shared TRGs showed strong enrichment for caste-DEGs. A striking example of caste-biased TRGs is a tandem array of three novel genes (Fig. 6A and SI Appendix, Table S11), RS001196, RS001197, and RS001198, for which no significant homologs were found in any organisms outside of termite clades in our sequence similarity search, although the German cockroach, Blattella germanica, seems to have a single copy of a homolog of these genes. These three genes were expressed at extremely high levels (up to 250,000 reads per kilobase per million [RPKM]) in R. speratus, which constituted ∼30% of the worker head transcriptome, and their expression was strongly biased across the three castes (Fig. 6A). Each gene encodes a short peptide of ∼60 amino acids in length that contained a secretion signal peptide in the N-terminal region followed by a middle part rich in charged amino acid residues and a C-terminal region rich in polar amino acids, with an unusually high number of tyrosine residues (Fig. 6B). Here, we refer to this novel class of peptides as the tyrosine-rich peptide family (TY family). The three TY genes shared modest sequence similarity with each other, suggesting that they are paralogs derived by tandem duplication. Duplicated TY family orthologs were also found in the genomes of Z. nevadensis, C. secundus, and M. natalensis (SI Appendix, Fig. S12). We estimated pairwise evolutionary rates (the ratio of nonsynonymous to synonymous substitutions, K_a/K_s) between R. speratus and Z. nevadensis for RS001196, RS001197, and RS001198. The K_a/K_s of each gene was 0.15, 0.16, and 0.03, respectively, indicating that they evolved under strong purifying selection and suggesting a conserved function in the termite lineage.

Recent advances in sociogenomics in different social insects are improving our understanding of the genetic bases of social evolution, which include the co-option of genetic toolkits of conserved genes, changes in protein-coding genes, cis-regulatory evolution leading to genetic network reconstruction, epigenetic modifications, accumulation of transposable elements, and TRGs (16, 22, 51, 88). In addition to these components, our genomic and transcriptomic analyses of R. speratus highlighted the significance of gene duplication for caste specialization. Gene duplication is, in general, a key source of genetic innovation that plays a role in the evolution of phenotypic complexity through sub- or neofunctionalization (26). Regarding eusocial evolution in insects, Gadagkar (89) first pointed out the importance of gene duplication, indicating that “genetic release followed by diversifying evolution” made possible the appearance of multiple caste phenotypes in social insects. According to Gadagkar's hypothesis, duplicated genes can be released from the constraints of original selection, leading to new directional evolution for caste-specific functions (e.g., queen- or worker-trait genes). Many decades later, genomic analyses revealed many instances of gene family expansions related to chemical communication, insulin, and vitellogenin signaling pathways in hymenopteran insects and termites (20, 21, 30–39), leading to the appreciation that gene duplications might have been commonly co-opted during social evolution in many social insects (23).

This study revealed that gene duplication associated with caste-biased gene expression is prevalent in the R. speratus genome. The list of duplicated genes encompasses a wide array of functional categories related to the social behaviors in termites as exemplified by transporters such as lipocalins (communication and physiological signaling; compare refs. 53 and 62), digestive enzymes such as CAZymes, immune-related genes such as lysozymes (social immunity), and metabolic enzymes such as GGPP synthase (social defense). This study also demonstrated that caste-specific expression patterns differed among in-paralogs. Although such paralogous genes were often observed in tandem in the genome, their expression patterns were often independent of one another, showing differential caste biases in many cases. Additionally, discordant caste biases in transcriptional expression were observed among closely related paralogs with similar coding sequences, as represented by the low correlation between phylogenetic position and caste specificity (Fig. 3A). Although the regulatory and evolutionary mechanisms underlying caste-biased expression patterns are elusive, these examples strongly suggest that gene duplications have facilitated caste specialization, leading to social evolution in termites.

After the gain of caste-biased gene regulation, subfunctionalization and/or neofunctionalization seems to have occurred, leading to caste-specific expression and caste-specialized functions. For example, in the case of the lipocalin family, lipocalin paralogs were generated by lineage-specific functional expansion in caste-specific organs or tissues: a queen-specific lipocalin (RS008881) was expressed specifically in the ovarian accessory glands, while a soldier-biased lipocalin (RS008823) was expressed exclusively in the frontal glands of soldier heads (Fig. 3 D and E). Taken together, these results lead us to hypothesize that, in termites, caste specification through gene duplication proceeded via following three steps: 1) gene family expansion by tandem gene duplication, 2) regulatory diversification leading to an expression pattern restricted to a certain caste, and 3) subfunctionalization and/or neofunctionalization of the gene products conferring caste-specific functions. As an exaptation of these steps, the case in which one (or several) of the multiple functions of pleiotropic genes is allocated and specialized to a duplicated gene copy might have led to caste-specific subfunctionalization (25, 26).

Recently, it was suggested that the evolution of phenotypic differences among castes of honeybees was associated with the gene duplication, based on the finding that duplicated genes presented higher levels of caste-biased expression than singleton genes (90). It was also shown that the level of gene duplication was correlated with social complexity in bees (superfamily Apoidea) (90). On the other hand, it was noted that a large number of genomics studies of hymenopteran species have not revealed such patterns (91, 92). Thus, it is still an open question whether gene duplication is a shared mechanism facilitating the evolution of caste systems in social insects with independent origins. Future studies targeting a broad range of taxa with different social systems will answer this question.

All mature colonies of R. speratus used for genome, RNA, and BS-seq were collected in Furudo, Toyama Prefecture, Japan (SI Appendix, Table S12). Detailed information on the samples is provided in SI Appendix, Supplementary Methodology.

We used female secondary reproductives (nymphoids I and II) for genome sequencing. Genomic DNA was isolated from each individual using a Genomic-tip 20/G (Qiagen). We examined five microsatellite loci to confirm whether the individuals were homozygous at these loci and shared the same genotype. We generated two paired-end libraries using a TruSeq DNA Sample Preparation Kit (Illumina) with insert sizes of ∼250 and ∼800 base pairs (bp) (SI Appendix, Table S13). Four mate-pair libraries with peaks at ∼3, ∼5, ∼8, and ∼10 kb were also generated using a Nextera Mate Pair Sample Preparation Kit (Illumina) (SI Appendix, Table S13). These libraries were sequenced using an Illumina HiSeq system with the 2 × 151-bp paired-end sequencing protocol. The reads of the paired-end and mate-pair libraries were assembled using ALLPATHS-LG (build No. 47878) (93) with the default parameters.

A protein-coding gene reference set of R. speratus was generated based on two main sources of evidence: aligned R. speratus transcripts and aligned homologous proteins of other insects and a set of ab initio gene predictions. The R. speratus RNA-seq reads were assembled de novo using Trinity (r2013_0814) (94) and then mapped to the genome using Exonerate (version 2.2.0) (95). We processed homology evidence at the protein level using the reference protein sequences of seven sequenced insects including Z. nevadensis and Blattodea protein sequences predicted from RNA-seq of Periplaneta americana and N. takasagoensis (detailed dataset information is provided in SI Appendix, Supplementary Methodology). These proteins were split-mapped to the R. speratus genome with Exonerate. These models were merged using the EvidenceModeler (r2012-06-25) (96), which yielded 15,584 gene models. A total of 74 genes were manually inspected and corrected. In particular, tandemly duplicated genes were prone to incorrect gene prediction with erroneous exon–exon connections across homologs. The final set of 15,591 genes was designated Rspe OGS1.0. The quality of OGS1.0 was evaluated by assessing two types of evidence: homology and expression data. Among the 15,591 genes, 12,996 (83.3%) showed hits in the National Center for Biotechnology Information nonredundant protein sequences database, 10,440 (70.0%) included known protein motifs defined in the Pfam database, and 14,302 (91.7%) showed evidence of expression with a threshold of RPKM = 1.0 in at least one sample of caste-specific RNA-seq data. In total, 15,577 genes (99.9%) showed evidence of the existence of homologs and/or expression.

Orthology determination among three termites: Orthologous genes among three termite species, R. speratus, Z. nevadensis, and M. natalensis (protein sequences of gene models RspeOGS1.0, ZnevOGSv2.229, and MnatOGS3, respectively), were determined by pairwise comparisons with InParanoid version 4.1, followed by a three-species comparison with MultiParanoid (97, 98). The M. natalensis gene set (MnatOGS3) was built in this study using a pipeline similar to that used for R. speratus.

Orthology analysis with arthropods: The orthology relationships of R. speratus genes (OGS1.0) with other arthropod genes were analyzed by referring to the OrthoDB gene orthology database version 8 (87 arthropod species) (99). We grouped R. speratus genes with the OrthoDB ortholog group using a two-step clustering procedure. For each R. speratus protein, BLASTP was used to find similar proteins among the arthropod proteins, and the ortholog group of the top hit was provisionally assigned to the query R. speratus gene. Then, ortholog grouping was evaluated by comparing the similarity levels (BLAST bit scores) among members within the focal ortholog group. We retained the grouping if the BLAST bit score between the query R. speratus gene and the top arthropod gene was higher than the minimal score within the original cluster members. Among the 15,591 R. speratus OGS1.0 genes, 12,434 were clustered into 9,033 OrthoDB Arthropod ortholog groups. Gene duplication was assessed based on this clustering. If two or more members of one species were included in a single ortholog group, it was regarded as a multigene family.

Old workers (sixth and seventh instars) and soldiers were collected from each colony. To collect primary reproductives, dealated adults were randomly chosen from each colony in accordance with the method described in the literature (100), and female–male pairs were mated (SI Appendix, Table S2). Kings and queens were sampled after 4 mo. Each individual was divided into the head and body (thorax + abdomen). We prepared RNA-seq libraries for 12 categories based on castes (reproductives, workers, and soldiers), sexes (males and females), and body parts (head and thorax + abdomen). Three biological replications of the 12 categories were performed from three different field colonies totaling 36 RNA-seq libraries (SI Appendix, Table S2). All Illumina libraries prepared using a TruSeq Stranded mRNA Library Prep kit were subjected to single-end sequencing of 101-bp fragments on HiSeq 2500. The cleaned reads were mapped onto the genome with TopHat version 2.1.0 guided by the OGS1.0 gene models. Transcript abundances were estimated using featureCounts and normalized with the trimmed mean of M-values algorithm in edgeR. DEGs among castes and between sexes were detected in each body part (head/thorax and abdomen) according to a GLM with two factors, namely, caste and sex using edgeR with the conditions set to a FDR < 0.01 and a log2-fold change of the expression level > 1.

To quantify caste specificity from the RNA-seq data, we used the tau score, which was originally developed to measure the tissue specificity of an expression profile (101). The index tau is defined as follows:where N is the number of categories (i.e., castes in this study), and $x_i$ is the expression level normalized by the maximal component value. We used RPKM+1 values as representations of expression levels.

See SI Appendix, Supplementary Methodology for whole-genome alignment, annotation of gene models/repeat sequences, methylome analysis, and RNA in situ hybridization.