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

Significance Gene duplication is a major source of evolutionary innovation and is associated with the increases in biological complexity and adaptive radiation. Termites are model social organisms characterized by a sophisticated caste system. We analyzed the genome of the Japanese subterranean termite, an ecologically and economically important insect acting as a destructive pest. The analyses revealed the significance of gene duplication in social evolution. Gene duplication associated with caste-biased gene expression was prevalent in the termite genome. Many of the duplicated genes were related to social functions, such as chemical communication, social immunity, and defense, and they were often expressed in caste-specific organs. We propose that gene duplication facilitates social evolution through regulatory diversification leading to caste-biased expression and functional specialization.

Termites are model social organisms characterized by a polyphenic caste system. Subterranean termites (Rhinotermitidae) are ecologically and economically important species, including acting as destructive pests. Rhinotermitidae occupies an important evolutionary position within the clade representing a transitional taxon between the higher (Termitidae) and lower (other families) termites. Here, we report the genome, transcriptome, and methylome of the Japanese subterranean termite Reticulitermes speratus. Our analyses highlight the significance of gene duplication in social evolution in this termite. Gene duplication associated with caste-biased gene expression was prevalent in the R. speratus genome. The duplicated genes comprised diverse categories related to social functions, including lipocalins (chemical communication), cellulases (wood digestion and social interaction), lysozymes (social immunity), geranylgeranyl diphosphate synthase (social defense), and a novel class of termite lineage-specific genes with unknown functions. Paralogous genes were often observed in tandem in the genome, but their expression patterns were highly variable, exhibiting caste biases. Some of the assayed duplicated genes were expressed in caste-specific organs, such as the accessory glands of the queen ovary and the frontal glands of soldier heads. We propose that gene duplication facilitates social evolution through regulatory diversification, leading to caste-biased expression and subfunctionalization and/or neofunctionalization conferring caste-specialized functions.
social insect j termite j gene duplication j caste system T he evolution of eusociality (i.e., animal societies defined by the reproductive division of labor, cooperative brood care, and multiple overlapping generations) represents an important step in animal evolution that increased the level of biological complexity (1,2). Eusocial insects such as bees, wasps, ants, and termites show sophisticated systems based on the division of labor among castes, which is one of the pinnacles of eusocial evolution (3). Recent advances in molecular biological technologies and omics studies have revealed many molecular mechanisms underlying eusociality and have led to the establishment of a new field of study known as "sociogenomics" (4). The genomes of major eusocial hymenopteran lineages (i.e., ants, bees, and wasps) have been sequenced, and the differences in gene expression and DNA methylation among castes have been explored (5)(6)(7)(8)(9)(10)(11)(12)(13)(14). These sociogenomic studies in hymenopterans have revealed some of the genetic basis of social evolution, including the co-option of genetic toolkits of conserved genes, Significance Gene duplication is a major source of evolutionary innovation and is associated with the increases in biological complexity and adaptive radiation. Termites are model social organisms characterized by a sophisticated caste system. We analyzed the genome of the Japanese subterranean termite, an ecologically and economically important insect acting as a destructive pest. The analyses revealed the significance of gene duplication in social evolution. Gene duplication associated with caste-biased gene expression was prevalent in the termite genome. Many of the duplicated genes were related to social functions, such as chemical communication, social immunity, and defense, and they were often expressed in caste-specific organs. We propose that gene duplication facilitates social evolution through regulatory diversification leading to castebiased expression and functional specialization. changes in protein-coding genes, cis-regulatory evolution leading to genetic network reconstruction, epigenetic modifications, and taxonomically restricted genes (TRGs) (15,16).
Termites (Isoptera) are another representative insect lineage exhibiting highly sophisticated eusociality and a wide range of social complexities (17). Termites are hemimetabolous, diploid insects that are phylogenetically distant from hymenopterans with holometaboly and haplodiploidy. Termite societies are characterized by reproductives of both sexes, workers, and soldiers. In the termite lineage, eusociality is thought to have evolved once, although the levels of social complexity features, such as colony size, feeding habitat, symbiosis with microorganisms, and caste developmental pathways, have diverged among termite species (18). These characteristics are especially different between the two major termite sublineages: the early-branching families (referred to as "lower" termites) and the most apical family Termitidae ("higher" termites) (Fig. 1A). Based on the wholegenome sequences of a few termite species reported to date, the commonality and diversity of the genetic repertoires between Isoptera and Hymenoptera or between termite lineages and their solitary outgroup (i.e., cockroach) have been investigated (20)(21)(22)(23). Additionally, in Zootermopsis nevadensis, clear differences in gene expression levels among castes (20) and in DNA methylation between alates and workers (24) have been detected.
Evolutionary novelties are often brought about by gene duplications (25) (reviewed in ref. 26). Gene duplication allows the subsequent divergent evolution of the resultant gene copies, enabling evolutionary innovations in protein functions and/or expression patterns (27)(28)(29). The transition to eusociality in Hymenoptera has been associated with gene family expansion, including the expansion of odorant receptor (30)(31)(32)(33)(34), insulin signaling (35), and vitellogenin genes (36,37). In two termites with available genome information, Z. nevadensis and Cryptotermes secundus, copy numbers were found to be expanded for genes such as ionotropic receptor (20,21), vitellogenin (20,38), insulin receptor (39), and juvenile hormone biosynthesis genes (40). Despite the accumulating examples of duplications of some categories of genes in social insects, their genome-wide impact and the evolutionary significance remain unclear.
Among the more than 2,900 extant species of termites (Isoptera) (41), subterranean termites (Rhinotermitidae), especially those belonging to two genera, Reticulitermes and Coptotermes, occupy an important evolutionary position (Fig. 1A). Recent phylogenetic studies showed that Rhinotermitidae is paraphyletic, and a clade including Reticulitermes, Coptotermes, and Heterotermes was shown to be sister to Termitidae (42,43). In particular, Reticulitermes exhibits intermediate social complexity between higher (Termitidae) and lower (all the other families) termites (44), although it belongs to the latter group [e.g., this genus displays primitive feeding ecology and gut symbiont features, a relatively complex colony structure and a caste development mode known as the bifurcated pathway ( Fig. 1B)]. Moreover, Reticulitermes is the most common termite group in palearctic (45) and nearctic (46) regions and includes major pests causing serious damage to human-made wooden structures (47). For these reasons, members of this genus are probably among the most studied termites (17). Nevertheless, despite their evolutionary, ecological, and economic relevance, subterranean termites remain an understudied group in terms of both genetics and genomics.
In this study, we targeted the Japanese subterranean termite Reticulitermes speratus. We conducted whole-genome sequencing, caste-specific RNA sequencing (RNA-seq) analysis and whole-genome bisulfite sequencing (BS-seq) of R. speratus to understand the genomic, transcriptomic, and epigenetic bases of the social life of this termite species. We also compared the omics data of R. speratus with those of sequenced higher and lower termites (Fig. 1A). Our integrative analyses revealed that gene duplications are often associated with caste-biased gene expression and caste-specific functions, which highlights the significant role of gene duplication in eusocial evolution in the termite lineage.

Results and Discussion
Genomic Features of R. speratus. 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 proteincoding 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).
Transcriptome Differentiation among Castes. 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 (  (19). R. speratus is marked in red, and the two termites mainly used for comparison in this study are marked with bold characters. (B) Developmental pathway of R. speratus. Both sexes have the same developmental pathway. There are two larval stages before molting into a nymph (with wing buds) or worker (no wing buds). There are six imaginal stages, and the sixth-stage nymphs molt into alates, which are primary reproductives (queen and king). Secondary reproductives (neotenics or nymphoids) differentiate from third-to sixth-stage nymphs. In the apterous line, there are at least five stages of workers. Some workers in the colony molt into presoldiers and soldiers. The female neotenics used for genome sequencing and the three castes used for RNA-seq are marked with asterisks. (C) Gene repertoire of R. speratus, Z. nevadensis, and M. natalensis categorized by orthology. For each species, genes were compared to those of other 88 arthropods and grouped into three classes: orthologs shared with other arthropods (labeled "arthropod"), orthologs shared with either of the other termites but with no orthologs in other arthropods (labeled "termite"), and orphan genes unique to each species (labeled "no hit"). Genes classified as "termite" and "no hit" are designated as TRGs. (D) High conservation of synteny between termite genomes revealed by dot plots generated by comparing R. speratus with Z. nevadensis and M. natalensis. Scaffolds longer than 2.0 Mb in the R. speratus assembly are used for plotting. Forward alignments are plotted in red, and reverse alignments are plotted in blue.

Shigenobu et al.
Genomic and transcriptomic analyses of the subterranean termite Reticulitermes speratus: Gene duplication facilitates social evolution 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)(57)(58)(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.
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).
Gene Duplication and Caste-Biased Gene Expression. 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 À16 ). 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 À16 , 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.
Multigene Families Related to Caste-Specific Traits 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. Lipocalins. Lipocalins belong to a family of proteins with molecular recognition properties, such as the ability to bind a range of small hydrophobic molecules (e.g., pheromones) and specific cell surface receptors and to form complexes with soluble macromolecules (60). A previous study identified a gene of the lipocalin family, SOL1, that is exclusively expressed in the mandibular glands of mature soldiers of the rotten-wood termite Hodotermopsis sjostedti (61). SOL1 is thought to function as a signaling molecule for defensive social interactions among termite colony members (53). Moreover, RNA-seq analysis showed that a lipocalin gene, Neural Lazarillo homolog 1 (ZnNlaz1), was specifically expressed in soldier-destined larvae in an incipient colony of Z. nevadensis (62). Gene function and protein localization analyses suggested that ZnNLaz1 was a crucial regulator of soldier differentiation through the regulation of trophallactic interactions with a queen. Thus, it was of interest that the lipocalin-related motif (Pfam PF00061; lipocalin/cytosolic fatty-acid binding protein family) was significantly enriched in the list of caste-DEGs (FDR < 0.05; (SI Appendix, Table S4).
We identified 18 lipocalin family genes in the R. speratus genome (Fig. 3 A-C and SI Appendix, Table S6). The number of lipocalin genes was larger than those in other insects (Fig. 3C). Phylogenetic analysis of lipocalin family genes identified in arthropods, including three termite species, revealed a highly dynamic evolutionary history of this protein family (Fig. 3A). Two subfamilies, namely, clades A and B, were observed to have experienced extensive expansion in the termite lineage ( Fig. 3 A  and C). The most drastic expansion was found in clade A, which included H. sjostedti SOL1. In this clade, 5, 8, 3, and 4 genes were identified in R. speratus, M. natalensis, Z. nevadensis, and C. secundus, respectively, and extensive and independent gene expansions occurred in each species. Clade B was also composed of genes with termite lineage-specific duplication. In many cases, these lipocalin genes were found in tandem arrays in the R. speratus genome (Fig. 3B). The inferred phylogenetic tree indicated that duplications in each clade occurred after the divergence of termites from a common ancestor. We detected signs of positive selection in the four R. speratus genes of Clades A and B (RS008823, RS008824, RS13913, and RS008884) (SI Appendix, Fig. S4).
A comparison of the transcriptome among castes revealed that most lipocalin genes (15 of 18) showed caste-biased gene expression (Figs. 2D and 3B). The caste specificity, however, varied among genes, regardless of sequence similarity and positional proximity in the genome. In particular, the expression levels of genes in clades A and B differed greatly among castes. For example, RS008823 and RS008824 displayed solder-specific expression, the expression of RS013912 was biased toward workers, and RS013913 was down-regulated in soldiers. RS008881 and RS008884 were exclusively expressed in the queen body (thorax + abdomen), while RS008882, a gene located next to the two aforementioned genes, showed quite different expression patterns, with high expression levels in heads, especially those of workers. These results indicate that termite lipocalin genes have undergone dynamic expansion in terms of the gene repertoire and regulatory To address the caste-specific function of the lipocalin paralogs, the expression patterns of several selected caste-biased lipocalin genes were examined by in situ hybridization (Fig. 3 D-H and  SI Appendix, Fig. S5). RS008881, a queen-biased lipocalin gene, was found to be expressed exclusively in the accessory glands of the ovary (Fig. 3D). The next gene on the same scaffold, RS008882, was shown to be specifically expressed in worker antennae and maxillary/labial palps (Fig. 3 F-H). RS008823, a soldier-biased gene, was expressed exclusively in the frontal gland cells of the soldier heads (Fig. 3E). Note that ovaries and frontal glands develop during postembryogenesis in a caste-specific manner (i.e., ovaries in queens and frontal glands in soldiers) in the pathway of caste differentiation in R. speratus. Antennae and maxillary/labial palps are not caste specific but are crucial sensory organs, especially for blind termite immatures, such as workers. Given that animal lipocalins generally work as carrier proteins (63), there is a possibility that focal termite lipocalins bind and convey some molecules to targets from caste-specific organs [e.g., egg-recognition pheromone and soldier defensive and/or inhibitory substances (64-66)], or participate in sensory reception, as observed for odorant-binding proteins (67). Cellulases. Lignocellulose degradation in termites is achieved by a diverse array of CAZymes produced by the host and its intestinal symbionts. The repertoire of CAZyme families in the genome of R. speratus did not show considerable differences from those of other nonxylophagous insects, such as honeybees and fruit flies (SI Appendix, Fig. S6). However, we found gene family expansion and expressional diversification of glycoside hydrolase family (GH) 1 and GH9 members ( Fig. 4A and SI Appendix, Table S7). The majority of GH1 and GH9 members are β-glucosidases (BGs; EC 3.2.1.21) and endo-β-1,4-glucanases (EGs; EC 3.2.1.4), respectively, which are essential for cellulose digestion in termites (68). We identified 16 GH1 paralogs (SI Appendix, Table S7). Such gene expansion of GH1 is also observed in the genome of other termites, but the reason for this expansion remains elusive (69). Although the phylogenetic tree divided these GH1 paralogs into seven distinct groups (clades A to G in Fig. 4A), most of them were tandemly located in the genome of R. speratus (Fig. 4 B and  C). The predominantly expressed BG gene was RS004136, and the expression of this gene was clearly biased toward the body (thorax + abdomen) in reproductives and workers (Fig. 4B). This gene formed a rigid clade with bona fide BGs reported from the salivary glands or midgut of termites (clade A in Fig. 4A) (70), suggesting that this gene is involved in cellulose digestion in R. speratus. Indeed, in situ hybridization analysis showed that RS004136 was specifically expressed in the salivary glands of workers ( Fig. 4 D and E and SI Appendix, Fig. S7A). Other GH1 members showed a wide variety of expression patterns across castes and body parts (Fig. 4B). Some of these genes might have diversified their functions, other than wood digestion that are related to termite sociality, as observed for egg-recognition pheromones (71). A typical example of such diversification was provided by RS004624, which was expressed specifically in the abdomens of queens (Fig. 4B). The peptide sequence of this gene showed a monophyletic relationship with that of Neofem2 of C. secundus (clade G in Fig. 4A), which is a queen-recognition pheromone that probably functions in the suppression of reproductive emergence (72). In situ hybridization showed that RS004624 was specifically expressed in the accessory glands of queen ovaries (Fig. 4 F and G and SI Appendix, Fig. S7B), suggesting that RS004624 is involved in enzymatic activities in queen-specific glands. Together with the results for a queen-biased lipocalin (RS008881), this finding indicates that the queen accessory glands may produce some queen-specific pheromones. Like lipocalins, GH1 paralogs are also typical examples of multigene family members participating in caste-specific tasks that may have been acquired by gene duplication resulting in neo-or subfunctionalization.
We found four paralogs of GH9 in R. speratus (SI Appendix, Fig. S8 and Table S7). Although several insect GH9 EGs have acquired the ability to hydrolyze hemicellulose (73), the neo-or subfunctionalization of termite EGs has yet to be clarified. Intriguingly, we found that the GH9 member RS006396 was weakly but uniformly expressed across all termite body parts and castes. This result suggests that some GH9 members also perform a function other than that of cellulase, as is the case for GH1. Lysozymes. The immune system of termites is of particular interest because the group living of termites with nonsclerotized and nonpigmented epidermis and microbe-rich habitats puts them at high risk for pathogenic infections (74). Thus, defense against pathogenic microbes is important for termites. In the R. speratus genome, we identified 251 immune-related genes (SI Appendix, Table S8). The repertoire and number of immune-related genes of R. speratus showed no large differences compared to those of other insect species with the notable exception of lysozymes (SI Appendix, Fig. S9). Lysozymes are involved in bacteriolysis through the hydrolysis of β-1,4-linkages in the peptidoglycans present in bacterial cell walls, and three distinct types of lysozymes, chicken-or conventional-type (c-type), goose-type (g-type), and invertebratetype (i-type) lysozymes, have been found in animals (75). We identified 13 and 3 genes encoding c-type and i-type lysozymes, respectively, and the number of lysozyme genes was larger than those in other insects (SI Appendix, Fig. S9 and Table S9). Phylogenetic analysis revealed that c-type lysozymes underwent extensive gene duplications in the sublineage leading to R. speratus (SI Appendix, Fig. S10A). Six c-type lysozymes formed a tandem array on scaffold_859 (SI Appendix, Fig. S10B), probably generated by repeated tandem gene duplication events. Interestingly, most of the c-type lysozyme genes showed caste-biased expression. Three genes (RS014698, RS100022, and RS100023) exhibited high expression levels compared to those of other lysozyme genes and were expressed in a soldier-specific manner, while RS100026 was expressed in a worker-specific manner, and RS100024 and RS100025 were highly expressed in both workers and soldiers (SI Appendix, Fig. S10B and Table S9). The differential expression patterns of the lysozyme genes of R. speratus may correspond to the division of labor among castes in terms of colony-level immunity.
It is also possible that duplicated lysozymes may have functions in processes other than immunity. A previous study indicated that the salivary glands of R. speratus secrete c-type lysozymes to digest bacteria ingested by termites through social feeding behavior (76). The same lysozyme genes are also expressed in the queen ovaries and eggs and play a role in egg recognition as proteinaceous pheromones in R. speratus (66,71). We could not find identical sequences of these lysozyme genes in our gene models, but these sequences were most closely related to RS002400 with 88% nucleotide identity, which occupied the basal position of the lineage-specific gene expansion (SI Appendix, Fig. S10A). GGPP synthase. Whole-genome comparison of R. speratus with Z. nevadensis and M. natalensis revealed a 270-kb R. speratusspecific fragment in scaffold_31, while the rest of this scaffold showed very high syntenic conservation among the three termites (Fig. 5A). We found that the R. speratus-specific region was encompassed by a tandemly duplicated gene cluster composed of 13 genes encoding GGPP synthase ( Fig. 5B and SI Appendix, Table S10). GGPP synthase catalyzes the consecutive condensation of an allylic diphosphate with three molecules of isopentenyl diphosphate to produce GGPP, an essential precursor for the biosynthesis of diterpenes, carotenoids, and retinoids (77)(78)(79). The extensive duplication of GGPP synthase paralogs observed in R. speratus is unusual because the genomes of other insects surveyed have only a single copy of GGPP synthase gene. The phylogenetic analysis of GGPP synthase homologs revealed two clusters, corresponding to a possibly ancestral group (including RS007484) and an apical group (including the other paralogs identified) (Fig. 5C). The latter cluster also contained some GGPP synthase paralogs obtained from the termitid Nasutitermes takasagoensis (56).
Transcriptome data indicated that all of the GGPP synthase genes except for RS007484, which was a member of the ancestral group in the phylogenetic tree, showed caste-biased expression, and caste specificity varied across the paralogs (Fig. 5B). Specifically, RS100010, RS007480, RS100012, RS100015, RS100016, RS100017, and RS007483 showed soldier-specific expression, while RS007481, RS007482, and RS100013 showed reproductivespecific expression (Fig. 5B). Several GGPP synthase genes have been identified in some termite species and are known to function in a caste-specific manner [e.g., the soldiers of N. takasagoensis synthesize defensive polycyclic diterpenes for use in chemical defense via a process mediated by high expression of the GGPP synthase gene in the frontal gland (80)]. It has been reported that the soldiers of Reticulitermes have a frontal gland in which diterpenes are synthesized, although their biological role is not fully understood (81-83). Consequently, it is possible that the soldier-specific GGPP synthases identified to date are involved in chemical defense. Indeed, in situ hybridization revealed that the soldier-specific GGPP synthase RS100016 was expressed exclusively in the soldier frontal gland, as shown in a previous study (84) (Fig. 5D and SI Appendix, Fig. S7C). It is also possible that reproductive-specific GGPP synthases are involved in the metabolism of other diterpenes, such as the synthesis of pheromones; in particular, RS007481 shows strong queenspecific expression in the thorax and abdomen and may play a role in the synthesis of queen substances.
Under the branch-site model of codon substitutions (85), significant positive selection was detected in five branches of the R. speratus GGPP synthase family tree (SI Appendix, Fig. S11): ancestral branches 1 and 2, and the branches leading to RS100017 (branch 3), RS100012 (branch 4), and RS007483 (branch 5). These results suggest that all GGPP synthase paralogs of R. speratus except the ancestral type RS007484 have experienced positive selection and finally acquired novel roles in the production of defensive and/or pheromonal substances.
The TY Family, a Novel Gene Family Restricted to Termites. 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.
Facilitation of Caste Specification by Gene Duplication. 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][31][32][33][34][35][36][37][38][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 castebiased 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 castebiased 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 castespecific 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.

Materials and Methods
Insects. 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.
Genome Sequencing and Assembly. 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 Tru-Seq 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.
Gene Prediction. 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 Inference and Gene Duplication Analysis. 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 Genomic and transcriptomic analyses of the subterranean termite Reticulitermes speratus: Gene duplication facilitates social evolution the OrthoDB gene orthology database version 8 (87 arthropod species) (99). We grouped R. speratus genes with the OrthoDB ortholog group using a twostep 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.
RNA-seq. 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.