Hemimetabolous genomes reveal molecular basis of termite eusociality

Abstract

Around 150 million years ago, eusocial termites evolved from within the cockroaches, 50 million years before eusocial Hymenoptera, such as bees and ants, appeared. Here, we report the 2-Gb genome of the German cockroach, Blattella germanica, and the 1.3-Gb genome of the drywood termite Cryptotermes secundus. We show evolutionary signatures of termite eusociality by comparing the genomes and transcriptomes of three termites and the cockroach against the background of 16 other eusocial and non-eusocial insects. Dramatic adaptive changes in genes underlying the production and perception of pheromones confirm the importance of chemical communication in the termites. These are accompanied by major changes in gene regulation and the molecular evolution of caste determination. Many of these results parallel molecular mechanisms of eusocial evolution in Hymenoptera. However, the specific solutions are remarkably different, thus revealing a striking case of convergence in one of the major evolutionary transitions in biological complexity.

Fig. 1. Phylogenetic, genomic and proteomic comparisons of 20 insect species.

From left to right: a phylogenetic tree of 20 insect species with Strigamia maritima (centipede) as the outgroup (species of newly sequenced genomes presented in this study are underlined); level of eusociality (one red insect: simple eusociality; two red insects: advanced eusociality; black fly: non-eusocial); fractions of repetitive content (yellow) within genomes of selected species (for sources, see Supplementary Information); proportions of species-specific gene family expansions (green), contractions (red) and stable gene families (black), the size of the pies represents the relative size of the gene family change (based on total numbers); a bar chart showing protein orthology across taxonomic groups within each genome. Ma, million years ago. Smar, Strigamia maritima; Edan, Ephemera danica; Rpro, Rhodnius prolixus; Nvit, Nasonia vitripennis; Amel, Apis mellifera; Pcan, Polistes canadensis; Hsal, Harpegnathos saltator; Lhum, Linepithema humile; Cflo, Camponotus floridanus; Pbar, Pogonomyrmex barbatus; Sinv, Solenopsis invicta; Aech, Acromyrmex echinatior; Acep, Atta cephalotes; Tcas, Tribolium castaneum; Aaeg, Aedes aegypti; Dmel, Drosophila melanogaster; Lmig, Locusta migratoria; Bger, Blattella germanica; Znev, Zootermopsis nevadensis; Csec, Cryptotermes secundus; Mnat, Macrotermes natalensis.

Fig. 2. Comparison of developmental pathways between B. germanica, the lower termites Z. nevadensis and C. secundus, and the higher termite M. natalensis.

Shown from left to right are: a simple phylogeny describing important novelties along the evolutionary trajectory to termites (numbers in brackets are genome sizes); life cycles; differential expression (log₂(fold change) > 1 and P < 0.05; DESeq2; sample sizes are shown in the last column) between workers and queens (between nymphs and adult females in B. germanica) of the selected gene families (Desat, desaturases; Elong, elongases; H’ween, Halloween genes) and total numbers within all genes; the numbers denote total numbers of genes in each gene family.

Fig. 3. Expansions, contractions and positive selection within iRs and ORs in termites.

a,b, IR (a) and OR (b) gene trees of 13 insect species. In each tree, only well-supported clades (support values > 85) that include B. germanica or termite genes are highlighted within the gene trees. The lengths of the coloured bars represent the number of genes per species within each of these clades. The red asterisk in a denotes the putative root of intronless IRs. c, The upper schematic diagram depicts the 2D structure of an IR, containing ligand-binding lobes (S1 and S2), transmembrane regions (TM1–3) and the pore domain (P). Below, the sequence of the domains along the peptide is represented, showing that the sites, which are under significant positive selection (red bars; codeml site models 7 and 8) within Blattodea IRs for M. natalensis (P < 1.7 × 10⁻¹⁰; likelihood-ratio test, 5 sequences, 413 aligned codons), are all situated within the ligand-binding lobes and on or around the putative ligand-binding sites (asterisks). d, The same representation for ORs, which include eight transmembrane regions. Positive selection was found for M. natalensis (P = 1.1 × 10⁻¹⁰; 5 sequences, 1,001 aligned codons) and C. secundus (P = 5.6 × 10⁻¹⁶; likelihood ratio test, 26 sequences, 1,913 aligned codons) of the orange clade, each at two codon positions within the second transmembrane region and at a third position within the carboxy-terminal extracellular region for M. natalensis.

Fig. 4. CpG_o/e of seven hemimetabolous insects.

a, Principal component analysis (PCA) of predicted DNA methylation patterns among 2,664 1-to-1 orthologues, estimated via CpG_o/e. The spheres represent the positions of the species within the 3D PCA, with the distance between the spheres representing the similarity of CpG_o/e between species at each orthologue; the curves are the distribution of CpG_o/e, with the dotted line showing CpG_o/e = 1. b, Tag clouds of enriched (P < 0.05; Fisher test, weight algorithm, topGO) GO terms (biological processes) among the lower (left) and the higher quartile (right) of CpG_o/e within termites (top) and B. germanica (bottom). For termites, genes were merged from all three species for analysing GO term enrichment. Number of enriched genes and total number of genes in topGO enrichment tests (low CpGo/e/high CpG_o/e/gene universe): B. germanica (3,291/1,842/11,409); termites (6,754/4,600/25,910). High CpG_o/e indicates a low level of DNA methylation and vice versa.