No obvious transcriptome‐wide signature of indirect selection in termites

Abstract

The evolution of sterile helper castes in social insects implies selection on genes that underlie variation in this nonreproductive phenotype. These focal genes confer no direct fitness and are presumed to evolve through indirect fitness effects on the helper's reproducing relatives. This separation of a gene's phenotypic effect on one caste and its fitness effect on another suggests that genes for this and other forms of reproductive altruism are buffered from selection and will thus evolve closer to the neutral rate than genes directly selected for selfish reproduction. We test this hypothesis by comparing the strength of selection at loci associated in their expression with reproductive versus sterile castes in termites. Specifically, we gather caste‐biased gene expression data from four termite transcriptomes and measure the global dN/dS ratio across gene sets and phylogenetic lineages. We find that the majority of examined orthologous gene groups show patterns of nucleotide substitution that are consistent with strong purifying selection and display little evidence for distinct signatures of direct versus indirect selection in reproductive and sterile castes. For one particular species (Reticulitermes flavipes), the strength of purifying selection is relaxed in a reproductive nymph‐biased gene set, which opposes the nearly neutral idea. In other species, the synonymous rate (dS) alone was often found to be the highest in the sterile worker caste, suggesting a more subtle signature of indirect selection or an altogether different relationship between caste‐biased expression and rates of molecular evolution.

Abstract

Termites are eusocial insects with reproductive and non‐reproductive castes, the latter of which can only evolve via indirect selection. We looked at whether genes associated with sterile workers and soldiers have similar patterns of molecular evolution to those of reproductive castes. One species, the Eastern subterranean termite (shown), yields dN/dS ratios indicative of relaxed purifying selection in queen‐biased (QB) genes relative to caste‐unbiased (UB) genes. This pattern was not universal across species but does suggest caste can affect strength of selection. Photo: Brock Fenton (Western University, Canada).

Open in new tabDownload slide

INTRODUCTION

Eusociality is a type of breeding system in which individuality is subsumed into an interdependent society with reproductive division of labour (Beshers & Fewell, 2001; Hölldobler & Wilson, 2009). For most eusocial groups of bees, wasps, ants (Hymenoptera) and termites (Isoptera), the division of labour is so pronounced that some individuals are physically specialized for reproduction, whereas others perform nonreproductive tasks or defend the colony (Queller & Strassmann, 1998; Thompson & Chernyshova, 2021). How these nonreproductive and otherwise self‐sacrificial qualities can evolve in insects is best understood from Hamilton's (1964) inclusive fitness theory, which describes how so‐called 'genes for altruism' can evolve indirectly via their effect on reproducing relatives who carry and transmit, but do not generally express, these genes (Bourke, 2011; Charlesworth, 1978; Parker, 1989). The notion of indirect selection on genes is therefore fundamental to our understanding of social evolution, yet we know little about how real genes respond to this type of selection at the sequence level.

This shortcoming from the field of insect sociobiology may stem in part from the rather short list of genes that have been unequivocally linked to reproductively altruistic behaviour or reproductively altruistic castes in social insects (Schwander et al., 2010; Thompson et al., 2013). This lack of precise knowledge from molecular sciences hampers the widespread study of social gene evolution. At the same time, however, transcriptomic and genomic data are rapidly emerging from the social Hymenoptera (Kapheim et al., 2015) and Isoptera (Harrison et al., 2018). Such data are amenable to testing ideas from inclusive fitness theory, provided that connections to the underlying theory are made (Akçay et al., 2015; Linksvayer & Wade, 2016; Mullen & Thompson, 2015).

One opportunity to bridge the gap between sociobiological theory and gene expression data is to assume that genes uniquely associated in their expression with altruistically sterile castes evolve primarily under indirect selection. If genes associated with inherently altruistic worker and soldier phenotypes have indirect effects (Crozier & Pamilo, 1996), then we can use lists of caste‐associated genes, of which there are many, as a proxy to test how direct versus indirect selection affects them at the molecular level. In so doing, we can potentially reveal any molecular signatures of indirect (i.e. 'kin') selection.

Even this proposed analysis, however, is constrained by uncertainty from population genomic models on how genes ought to evolve under caste‐mediated selection (Helanterä & Uller, 2014). Linksvayer and Wade (2009) predict that, all else equal, genes indirectly selected for nonreproductive helper traits (i.e. genes associated in their expression with workers or soldiers) will show signatures of having evolved under relatively weak selection—that is, weak positive or weak purifying selection—and correspondingly harbour more deleterious polymorphisms, relative to genes associated with the direct selection of reproductive traits. If so, we expect worker‐ or soldier‐associated genes from any eusocial insect to more closely approach the neutral rate, as might generally be expected for genes under weak selection (Kimura, 1983; Nielsen, 2005). This argument is similar to that proposed by Hall and Goodisman (2012) who predict weaker selection on worker‐associated genes, in particular for species with low nestmate relatedness.

In our study, we propose to test for molecular signatures of kin selection by comparing the rates of adaptive molecular evolution across caste‐biased and unbiased genes from social insects in the order Blattodea. Termites—still colloquially referred to as 'Isoptera'—are a relatively small (~2,500 spp.) epifamily (Termitoidae; Eggleton et al., 2007) of social roaches that live in caste‐based societies (Roisin, 2000). Sexual kings and queens are specialized for reproduction, whereas more‐or‐less sterile workers and soldiers perform nonreproductive roles that are associated with colony growth, colony maintenance and defence (Shellman‐Reeve, 1997). Workers and soldiers can therefore be considered reproductively altruistic because they help produce large numbers of nondescendent kin (i.e. siblings, half‐siblings.) at the expense of their own direct fitness (West et al., 2007).

If the nearly neutral idea is applied to the termite data, we expect that queen‐associated genes will evolve under strong selection and thus show a pattern of nonsynonymous (dN) to synonymous (dS) substitution rate that deviates strongly from a neutral expectation of dN/dS = 1, either in the positive dN /dS > 1 or purifying dN /dS < 1 direction. Conversely, we expect soldier‐ and worker‐associated genes to show evidence of having evolved under weak or relaxed selection and thus not deviate as strongly from neutrality in either positive or purifying direction.

Other scenarios outside of this rationale are also possible. Harpur et al. (2014), for example, suggest that it is the nonreproductive castes (specifically, honey bee workers) that may well be more exposed to and affected by selection than are queens that remain in‐hive and are shielded from the outside environment by the workers themselves. If so, then worker‐associated genes may be under relatively strong selection. This adapted worker hypothesis differs from the nearly neutral hypothesis and generally predicts strong, not weak, selection on worker‐biased genes. Importantly however, neither idea has been tested outside of the social Hymenoptera. Yet, another possibility is for genes biased in their expression towards any caste, sterile or otherwise, to evolve rapidly under drift or selection (Gadagkar, 1997; Hunt et al., 2010), once freed from the constraints of multi‐caste pleiotropy, relative to caste‐unbiased genes. Each of these three theoretical possibilities has intuitive appeal but a general signature of indirect or caste‐mediated selection, if there is one, remains unknown.

METHODS

Defining caste‐biased genes

To identify sets of caste‐associated genes, we made use of four recently published, caste‐specific RNA‐seq data sets, corresponding to Zootermopsis nevandensis (Archotermopsidae; Terrapon et al., 2014), Cryptotermes secundus (Kalotermitidae; Harrison et al., 2018), Reticulitermes flavipes (Rhinotermitidae; Wu et al., 2018) and Macrotermes natalensis (Termitidae; Harrison et al., 2018). Despite the differences in sampling and analysis among these source studies that we used to infer caste‐associated genes, our meta‐genomic data set is ultimately derived from whole‐body tissue extracts of a small number (three‐four individuals) of each caste, typically sampled across colonies and populations (Table S1). From each data set, we aimed to identify genes with expression biased towards sterile castes (workers and soldiers) or towards reproductive castes (queens and nymphs). We also identified genes from each species without any caste‐biased expression. We defined caste‐biased genes as those with significantly higher log fold expression (adjusted p < .05; Wald test in DESeq2, v1.28; Love et al., 2014) in one caste relative to all others in each data set. By contrast, we considered genes to be unbiased if their pairwise comparisons had no significant differences in gene expression across all castes. The differentially expressed gene sets for C. secundus, Z. nevadensis and M. natalensis were obtained from Harrison et al. (2018) and for Reticulitermes flavipes from Wu et al. (2018).

Identifying orthologs

For evolutionary analyses, it was necessary to identify orthologous genes between the four termite species, together with those from the German cockroach, Blattella germanica, and the invasive walking stick, Medauroidea extradendata. To initiate this search, we first obtained gene annotation files for each termite as well as from the cockroach (Harrison et al., 2018) and walking stick (Brand et al., 2018) and then extracted the proteomes and coding sequences from each genome as presented in the original study or otherwise using GffRead (Pertea & Pertea, 2020). In all cases, orthologous groups were identified using OrthoFinder (Emms & Kelly, 2015) at default settings. We retained only those groups containing at least one gene copy per species. Where more than one was present, we retained only the longest copy, to maximize the power of our comparative sequence analysis and to ensure this analysis was consistently based on an alignment of exactly six orthologous genes.

Evolutionary analyses

For each ortholog set, amino acid sequences were aligned with mafft (v7.397; Katoh & Standley, 2013) using the E‐INS‐I strategy, and coding sequences were aligned using pal2nal (Suyama et al., 2006). Evolutionary rates (dN /dS) were estimated for each orthologous gene set using the 'free ratio' model of codeml within the paml suite (Yang, 1997), which allows dN/dS rates to vary throughout the species tree. To minimize saturation effects, alignments with values for dS > 3 or dN/dS > 10 were excluded from further analyses. For each species, dN/dS ratios on the corresponding terminal branches were compared between gene groups. Differences in dN/dS ratio between caste‐biased gene groups were tested with the Wilcoxon–Mann–Whitney test adjusting for false discovery rate (FDR) using the Benjamini–Hochberg procedure in r (v3.5.1; R Core Team, 2015).

To test for positive selection along specific branches of the phylogenetic tree, we used the adaptive branch‐site random effects likelihood method (Smith et al., 2015), as implemented in the aBSREL program within the HyPhy suite (v2.5.5; Kosakovsky Pond et al., 2020). Here, we allowed one of three rate classes to assume a dN /dS value of ω > 1 and applied likelihood ratio tests (LRTs) to gauge significance of the adaptive model's fit to our data relative to a null model in which ω does not exceed 1. Further, to test for evidence of relaxed or intensified purifying selection along specific branches of the phylogenetic tree, we applied a basic species‐level tree ((((((M. natalensis, R. flavipes), C. secundus), Z. nevadensis), B. germanica), M. extradentata); after Bourguignon et al., 2014) and the selection‐testing method of Wertheim et al. (2015), as implemented in the RELAX module of HyPhy. Here, we fit a three rate‐class codon substitution model (ω₁ ≤ ω₂ ≤ 1 ≤ ω₃) to each termite‐specific test branch (T) and to the remaining subset of reference (R) branches on the phylogenetic tree. For each species‐specific caste‐biased gene set, we compared the fit of a null model in which a selection intensity parameter k is constrained to a value of 1, against an alternative model in which this parameter is free to vary, such that $ωRk$ = ω_T. Caste‐biased genes for which k differs significantly from 1 indicate that the strength of selection has been intensified (k > 1) or relaxed (k < 1). We repeated this test for each caste‐biased gene set over the corresponding species‐specific branches of the tree. Finally, we tested for differences in the proportion of genes detected as evolving under positive, relaxed or intensified purifying selection using sets of chi‐squared tests (chisq.test), FDR‐corrected in r.

Determining protein functions

We inferred functions of the ortholog groups using evidence combined from four different methods: (1) Orthology to Drosophila melanogaster was determined by performing a blast search using the M. natalensis gene from each of our ortholog gene sets against the D. melanogaster proteome (v6.12; flybase.org) and inferring each gene's function from its 'reciprocal best hit' (Ward & Moreno‐Hagelsieb, 2014). (2) Where no clear one‐to‐one match (i.e. similarity score of E‐value < 10^–5) could be determined, protein sequences were blast‐ed against the NCBI nr database (latest version December 2019), and functional descriptions contained on uniprot.org were consulted. (3) All protein sequences were annotated with pfam domain collection of protein families database v32 using pfamscan (Finn et al., 2016). (4) Gene ontology (GO) terms were mapped to identify pfam domains based on the Pfam2GO table, as accessed from geneontology.org. For small gene sets (n < 10), individual genes were described based on orthology to known genes (1 and 2 above). For larger gene sets, functional enrichment analyses were carried out with the TopGO package (Alexa & Rahnenfuhrer, 2020) in r using the classic algorithm with a p‐value cut‐off at .05 and retention of GO terms with at least two genes mapped.

RESULTS

Caste‐associated genes

To enable comparative sequence analyses, we retrieved a set of 5,304 ortholog groups, containing at least one gene copy within each of four termite and two outgroup species (Table S2). Within this set of orthologs, we identified a total of 2,349 caste‐biased genes (44%; Table S3). Of these, a majority (1981, 84%) were up‐regulated in the reproductive caste with the remaining minority up‐regulated in nonreproductive workers (301) or soldiers (171), depending on the species (Figure 1). Each of the four species tended to express unique caste‐biased genes, as evidenced by the low proportion of overlap. In the queen‐biased set, for example, only 64 genes (~3%) were consistently queen‐biased across three or more species. Similarly, in the worker‐biased set only two genes (<1%; amylase‐distal and scabrous; Table 1) were shared between three termite species and no genes were worker‐biased across all four species. An additional two genes were determined as worker‐biased in the two Neoisoptera species (R. flavipes, M. natalensis)—namely the odorant binding protein Obp56d and the Gram‐negative binding protein GNBP1.

FIGURE 1

Venn diagrams depicting numbers of caste‐biased genes across four termite species. The majority of caste‐biased genes are uniquely expressed in queens. There are few identified core caste‐biased genes, most of which demonstrate species‐specific expression

Open in new tabDownload slide

The total queen‐biased gene set was enriched (p‐values < .05 with a minimum of n = 2 genes per GO category) for 89 GO terms related to a broad range of global functions including cellular macromolecule metabolism, gene expression, cell cycle and signalling (Table 2). The only five core genes that are queen‐biased in all four termites were histone H2B, MrgBP of the TIP60 complex (chromatin remodelling), CycB3 (cell cycle), Skp2 (cell cycle) and yolkless (yl, uptake of vitellogenin; Table 1). By contrast, worker‐biased genes were enriched for 25 GO terms with functions mainly related to carbohydrate metabolism, oxidation–reduction processes and transmembrane transport. We identified 26 genes that were soldier‐biased in both M. natalensis and R. flavipes, indicating their possible general importance for derived termite soldiers (Table 1). The total soldier‐biased set was enriched for protein homo‐oligomerization and cellular homeostasis.

Evolutionary analyses

The average evolutionary rate was well below a neutral value of ω = 1 for all gene sets and species, suggesting that these genes are generally conserved with respect to amino acid altering substitutions (Figure 2). Further, we found no widespread differences in this rate between sets of caste‐biased genes in any of the four species, with the exception of R. flavipes in which the median dN/dS ratio was higher in reproductive nymph‐biased versus unbiased genes (median dN/dS = 0.178 versus 0.133; U = 161,260; p = .016). We did, however, find other caste‐biased patterns of nucleotide substitution. The synonymous substitution rate (dS) was consistently highest in worker‐biased genes, significantly so in the dampwood (Z. nevadensis) and drywood (C. secundus) species, relative to queen‐biased and unbiased genes (Figure 2). Similarly, a higher nonsynonymous substitution rate (dN) was detected in worker‐biased genes of Z. nevadensis, relative to queen‐biased and unbiased sets. A similar pattern for dN was also observed in R. flavipes.

FIGURE 2

Violin plots showing variation in synonymous dS and nonsynonymous dN substitution rates across genes sets with queen‐biased (QB), worker‐biased (WB), soldier‐biased (SB) or unbiased (UB) expression for each of the four termite species. Letters above the plots denote differences between tested gene groups, if any (Wilcoxon–Mann–Witney U tests). 'QB' in R. flavipes refers to the sampled reproductive nymph caste, rather than true queens, which is present in these sometimes‐invasive species (Scaduto et al., 2012)

Open in new tabDownload slide

Within an overall context of purifying selection, we examined how gene‐specific estimates of dN/dS might vary over lineages of the termite tree, as might be expected, for example, if specific genes experienced changes in the strength or direction of selection at different points in termite social evolution (Bulmer & Crozier, 2006). We found evidence for a mix of positive selection, relaxed or intensified purifying selection on single genes over most of the six prespecified branches of the termite phylogeny (Figure 3). For example, a relatively high proportion of genes were under positive selection on the terminal branch leading to the socially complex and fungus‐growing termite M. natalensis (12.1%–16.1% of genes compared to only 0%–7% of genes on other terminal branches, depending on caste). Other genes examined are under positive selection, including two queen‐biased genes (OG0007671: Galphas and OG0003652: putative ortholog of NRIF1; Table 1) at the base of our working termite tree (‘Lower termite root’ in Figure 3) and one soldier‐biased gene (OG000361: aspartyl/asparaginyl beta‐hydroxylase) on the branch leading to the Neoisoptera (‘Higher termite root’ in Figure 3).

FIGURE 3

Number of genes under positive selection or under relaxed or intensified purifying selection within each caste‐biased set and relative to an unbiased gene set, as inferred from aBSREL and RELAX analyses over specific branches of termite phylogeny. Each bar graph shows the number of genes against the proportion tested in that set displayed on the y‐axis. Results of chi‐squared tests are shown above each bar graph. Total gene numbers are indicated below the header of each analysed branch

Open in new tabDownload slide

Genes under intensified purifying selection are evident at the root of the termite tree (the worker‐biased gene scabrous and 19 queen‐biased genes) and at the base of Neoisoptera (the queen‐biased gene yolkless, the worker‐biased GNBP1 and the soldier‐biased upheld). By comparison, caste‐biased genes under relaxed purifying selection were limited to the four terminal branches, with a majority (115 of 123) being queen‐biased (Figure 3). The functions of all single genes under selection are listed in Supplementary Files.

In most cases, the proportions of genes that are evolving under different selective regimes do not vary between caste‐biased groups. This homogenous pattern suggests that sets of genes associated in their expression with one caste or another are subject to similar selective forces, with one exception: the proportion of genes under intensified purifying selection differed among queen‐, worker‐ and unbiased genes at the root of the termite phylogeny (chi‐square = 12.6; P₂ = .002; Figure 3). Overall, however, the strength and direction of gene‐level selection does not typically vary by caste in termites and, by extension, a clear signature for direct or indirect selection remains obscure.

DISCUSSION

We applied a hypothesis testing framework to reconcile patterns of nucleotide substitution in termite transcriptomes with evolutionary scenarios that describe how selection drives the evolution of castes. Specifically, we test for signatures of direct versus indirect selection on genes associated in their expression with reproductive and nonreproductive castes in termites. From a relatively small but evolutionary diverse sample of four termite species, we found no clear pattern of gene‐level selection associated with reproductive (queen, nymph) or nonreproductive (worker, soldier) castes. Instead, we found that all caste‐biased genes showed similar rates of molecular evolution, with the predominant theme of strong purifying selection and lack of widespread evidence for caste‐associated differences in the strength or direction of selection. We did, however, determine that nymph‐biased genes in one species (R. flavipes) were evolving closer to the neutral rate than caste‐unbiased genes, in a manner suggesting relaxed purifying selection in that quasi‐reproductive caste. Moreover, we report that synonymous dS and nonsynonymous dN rates were often elevated in workers, but not in soldiers, which suggest a more complex association between sterile caste, gene expression and substitution rate than is currently understood. Finally, a gene‐level test for sites positively selected or, alternatively, sites under relaxed or intensified selection, identified dozens of genes with particularly complex signatures of selection. Overall, however, we did not find a consistent difference in the pattern of nucleotide substitution between sets of caste‐biased genes across termite species. At this point, we suggest that any differences in relative strength of gene‐level selection between termite castes are likely species‐specific and potentially related to ecology or life type rather than universally related to direct versus indirect selection. We note that the synonymous rate dS is elevated for worker‐biased genes of the two species (C. secundus, Z. nevadensis) with a 'one‐piece' ecological life type (Abe, 1991; Shellman‐Reeve, 1997) whereby colonies consume solely the wood in which they nest. Future studies could explore the potential for this and other (e.g. Korb et al., 2015) broad ecological correlates to termite gene diversity.

Few core caste‐biased genes

Our initial summary of genes showing caste‐biased expression in four termite species revealed few core genes with conserved caste function. Just five genes (of 1981) were queen‐biased in all four species, whereas for worker‐biased genes this number was zero (Figure 1). Such low overlap and strong species‐specific pattern in gene expression is consistent with the idea that different social animals use socially acting genes in novel ways (Behl et al., 2018; Sumner, 2014) and that caste‐biased expression is not necessarily conserved even in closely related species of termites (Weil et al., 2009). The particularly low overlap for worker‐biased genes suggests that the worker caste is a relatively generic phenotype with few caste‐defining genes expressed across termite species. Unlike the defensively specialized and monophyletic soldier caste, this lack of consistent gene expression signature in workers (Behl et al., 2018; Wu et al., 2018) may, in part, be due to their uncertain evolutionary history (Bignell & Roisin, 2011; Thompson et al., 2004). Nonetheless, we discovered two genes with relatively consistent worker bias—amylase‐distal is a digestive enzyme involved in the hydrolysis of starch and is worker‐biased in all but one species (C. secundus); and scabrous with neurogenic function, which likewise is worker‐biased in all but one species (R. flavipes). These two genes may therefore be important for the molecular biology of the worker caste. We also identified two genes that were consistently worker‐biased only in two socially complex species, R. flavipes and M. natalensis. These genes code for an odorant binding protein (Obp56d) and a bacteria binding protein (GNBP1), suggesting the importance of chemical communication and innate immunity, at least within more complex termite colonies.

No obvious transcriptome‐wide signature of indirect selection

Contrary to our expectation, the strength of gene‐level selection in termites on sets of orthologous genes does not typically vary between genes unbiased in their expression versus those biased towards a specific reproductive or nonreproductive caste. If indirect selection is generally weaker than direct selection, as if diluted by kinship between actor and receiver (Hamilton, 1964), then all else equal (Hall & Goodisman, 2012; Linksvayer & Wade, 2009) we expect worker‐ and soldier‐associated genes to more closely approach the neutral rate than queen‐biased or unbiased caste genes. Instead, and unlike one study on the social Hymenoptera (Warner et al., 2017), we found little evidence that selection on sterile castes is less efficient than on reproductive castes. Our results are more similar to Imrit et al. (2020) who found no consistent evidence from four social (hymenopteran) insect species that the strength of purifying selection was stronger on queen‐biased genes relative to worker‐biased genes. This convergent result with Imrit et al. (2020) suggests generality in our findings across several origins of eusociality.

One possibility for termite species is that the strength of selection is similar on both directly and indirectly selected loci. For example, many termite species tend to inbreed, either through close‐relative reproduction within their colonies or upon consanguineous mating at the time of colony foundation (e.g. Thompson et al., 2007). These in‐breeding patterns might facilitate the purging of recessive lethal mutations and contribute to the consistently strong pattern of purifying selection that we observe. Importantly, however, in one species 'queen'‐biased genes (technically, nymph‐biased in R. flavipes) did show some evidence of relaxed purifying selection relative to unbiased genes. That is, the reproductive caste's nymph‐biased genes in R. flavipes edged closer to the neutral rate than did unbiased genes (Figure 2).

Although this pattern of relaxed purifying selection on genes associated with the reproductive rather than the helper caste is inconsistent with the nearly neutral hypothesis, it also does not provide support for the adapted worker hypothesis. In some social Hymenoptera, including honey bees (Harpur et al., 2014) and fire ants (Roux et al., 2014), signatures of positive selection are obvious on worker‐biased genes. This suggests that in some taxa, workers are not buffered from the full‐strength of selection, rather they appear to be the focal targets of adaptive diversification. This is not a pattern we observed in our current termite data set. First, worker‐biased genes are generally under purifying rather than positive selection. Second, as a matter of degree, worker‐biased genes are not actually different in their selection profile from sets of caste‐unbiased genes for any termite species tested here. A similar pattern in which the evolutionary rate does not vary between worker‐biased and caste‐unbiased genes was determined in a primitively eusocial Polistes wasp (Ferreira et al., 2013). Our analyses contained two one‐piece species (C. secundus, Z. nevadensis), where workers retain some reproductive potential as pseudergates, and two separate‐type species (R. flavipes, M. natalensis) in which workers are well differentiated from reproductive castes and are effectively sterile (Abe, 1991; Shellman‐Reeve, 1997). Despite these and other natural history details that may influence the relative strength of selection on castes (e.g. colony size, queen longevity, mode of colony founding), dN/dS did not vary. Moreover, termite soldier‐biased genes in our study, with conditional expression and presumably strong indirect effects of their own, show no evidence of relaxed selection – a pattern, which again opposes the nearly neutral idea. Instead, the relaxed purifying selection that we observed in nymph‐biased genes of R. flavipes is nominally consistent with the genetic release hypothesis of Gadagkar (1997), with the exception that information on pleiotropy is absent and here the 'release' from purifying selection is limited to just one reproductive caste (i.e. nymph), rather than all castes.

Estimates of dN/dS rates across entire gene sets can reveal average differences in selection but a large proportion of genes under, say, strong purifying selection (dN/dS << 1) can mask individual genes that deviate from this pattern. Our gene‐specific test for selection revealed many loci with such deviant patterns of selection—for example, worker‐associated genes (Apoltp and chitin deacetylase‐like 5) under relaxed purifying selection on the branch leading to Z. nevadensis or on the branch leading to M. natalensis (E3 ubiquitin ligase complex and DPR interacting protein zeta). Other examples include genes under relaxed purifying selection for queens, from a low of six genes in M. natalensis to a high of 68 genes in Z. nevadensis, and likewise some relaxed genes in soldier caste (Cyp4g15 in M. natalensis and three uncharacterized genes in R. flavipes). These and other interesting deviations associated with intense purifying or positive selection (Figure 3) are nonetheless consistent with our global estimates of dN/dS that show no consistent difference in strength of selection across genes associated in their expression with reproductive versus nonreproductive castes. Rather, the differences in selection that we do detect are likely associated with specific lineages and species.

The synonymous rate dS on its own, however, was determined to be significantly higher in workers among the two most basal termite species, relative to both queen‐biased and unbiased genes (Figure 2). Under the assumption that synonymous substitutions are neutral, a higher dS value may be caused by higher rates of synonymous mutations or, alternatively, by similar mutation rates but with less efficient DNA repair mechanisms for worker‐associated genes. Similarly, the nonsynonymous rate dN was often highest in workers versus other castes. This tendency for higher dN and dS in worker‐associated genes, whereas dN/dS ratio for this caste remains low, suggests that worker‐associated genes may harbour more polymorphisms despite the presence of strong purifying selection for this caste. We suggest that elevated dN and dS rates under strong purifying selection in workers may be explained by a reduced DNA repair rate for genes expressed in this caste. With low expression in queens, worker genes may be less accessible to the DNA repair machinery, and thus, queens may pass unrepaired copies to future generations in a manner analogous to the error‐prone transmission of sex‐biased genes in humans (Supek & Lehner, 2015). If so, this type of molecular signature may be linked to caste‐biased expression but is not a precise diagnostic of direct and indirect selection and, instead, may accrue as a by‐product of neutral processes (Helanterä & Uller, 2014).

CONCLUSIONS

If direct selection is generally stronger acting than indirect, then we expect genes associated in their expression with nonreproductive altruistic helper castes to evolve as if selection were, among other considerations, diluted through relatedness to those they are helping. As such, we expect worker‐ and soldier‐associated genes to more closely approach the neutral rate relative to queen‐biased genes, from either positive or purifying direction (Chernyshova et al., 2018). By comparing worker‐, soldier‐, queen‐associated and caste‐unbiased genes across four termite species, we present a preliminary test for this sociogenomic hypothesis. We did not reveal any obvious signatures of indirect selection. Our transcriptome‐wide tests measuring strength of selection were either not significant (Figure 2a, b and d) or showed that queen‐biased genes, not worker‐ or soldier‐biased genes, are under relaxed purifying selection (Figure 2c). This latter pattern is very similar to the one observed between adult castes in the fire ant Solenopsis invicta (Hunt et al., 2011). Future tests might therefore increase power by focusing on single‐gene patterns, potentially by use of candidates identified in our study (e.g. Galphas, NRIF1, scabrous, aspartyl/asparaginyl beta‐hydroxylase, yolkless, GNBP1, upheld). This candidate gene approach seems reasonable given that the number of specific genes under selection to maintain caste phenotypes or to regulate caste‐biased expression might be few (Figure 3; Thompson & Chernyshova, 2021). To ensure that future comparative sequence analyses are comparing 'like with like', candidate genes would ideally include paralogues with common‐origin copies that are functionally differentiated by caste within a single species—if such genes exist in sufficient number (Chau & Goodisman, 2017). Alternatively, genome‐wide analyses could focus on single species and make use of population level data, which conciliate the all else equal criterion and help isolate the directness effect of selection on protein divergence.

To render our tests operational with interspecific data in this study, we made several simplifying assumptions. First, that our gene sets are essentially comparable with respect to the mean fitness distribution, effective population size, evolutionary history and network connectivity. Second, we assumed that each caste‐biased gene set is exclusively expressed in one caste or another and thus evolves under direct or indirect selection and not a mix of the two. Therefore, in our current analysis, the indirect effect would indeed need to be large to be clearly detected. Given that our gene sets are ultimately derived from whole‐body tissue extracts, some important tissue‐specific and thus caste‐specific genes may inadvertently be absent from our analysis. Our choice of developmental stage (adult or late instar, depending on caste) may also have obscured particular gene expression patterns, and we acknowledge the difficulty of identifying truly caste‐biased genes. Notwithstanding these and other potential complexities associated with tests for indirect selection (Linksvayer & Wade, 2016), our interspecific analyses do reveal interesting and informative patterns—namely that substitution rates do vary by caste, albeit mostly in a species‐specific manner. Our study thus highlights an emerging challenge in insect sociobiology—to test and explain the precise relationship between caste‐biased gene expression and corresponding rates of molecular evolution.

ACKNOWLEDGMENTS

This work was supported by a Natural Sciences and Engineering Research Council (NSERC) of Canada Discovery Grant to GJT and a DFG grant (BO2544/11‐1) to MCH. The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

All three authors conceived the project at the World Congress of the International Union for the Study of Social Insects meeting (IUSSI2018) in Guarujá, Brazil. MCH and AMC assembled and compiled the data. MCH analysed the data. GJT and MCH drafted the paper. All three authors edited and finalized the paper.

Peer Review

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13749.

DATA AVAILABILITY STATEMENT

Our data are publicly available at GenBank https://www.ncbi.nlm.nih.gov/ under the following accessions: C. secundus: PRJNA382129; M. natalensis: PRJNA382034 (published in Harrison et al NatEcoEvo 2018) Z. nevadensis: PRJNA203244 (published in Terrapon et al., 2014). R. flavipes: PRJNA379073 (published in Wu et al., 2018). HyPhy scripts used in our analysis are available and can be accessed with the following link: https://github.com/MCH74/TermiteSelection/

REFERENCES