No Evidence for Single-Copy Immune-Gene Specific Signals of Selection in Termites

Selection pressures from pathogens appear to play an important role in shaping social evolution. Social behavior, in particular brood care, is associated with pathogen pressure in wood-dwelling “lower” termites. Yet, generally pathogen pressure is predicted to be low in wood-dwelling termite species that never leave the nest except for the mating flight. In comparison, pathogen pressure is predicted to be higher in species that leave the nest to forage, and thus constantly encounter a diversity of microbes from their environment. We hypothesized that such differences in predicted pathogen pressure are also reflected by differences in the intensity of natural selection on immune genes. We tested this hypothesis in a phylogenetic framework, analyzing rates of non-synonymous and synonymous substitutions on single-copy immune genes. Therefore, we leveraged recent genomic and transcriptomic data from eight termite species, representing wood-dwelling and foraging species as well as 14 additional species spanning the winged insects (Pterygota). Our results provide no evidence for a role of pathogen pressure in selection intensity on single-copy immune genes. Instead, we found evidence for a genome-wide pattern of relaxed selection in termites.


INTRODUCTION
Like in other organisms, pathogens seem to be important drivers of evolution in social insects. On the one hand, social insects are a "desirable" target for pathogens as an insect colony represents a large source of many potential hosts, all with a similar genetic background. Thus, if pathogens manage to enter a colony, they can exploit many individuals. However, social insects are also well protected as they evolved "social immunity" (Traniello et al., 2002;Cremer et al., 2007), a repertoire of defensive mechanisms that work at the colony level. For instance, behavioral task division limits contact to potentially infected individuals and can lead to their eviction. Also molecular mechanisms of social immunity exist, for example indirect immunization of colony members (e.g., Traniello et al., 2002;Cremer et al., 2007;Masri and Cremer, 2014) or impregnation of the nest walls with fungicidal compounds (Bulmer et al., 2009;Rosengaus et al., 2011). Thus, social immunity can be considered a selected emergent property of insect colonies where the whole is more than the sum of the individual parts (Rosengaus, personal communication). The evolution of social immunity aligns with the complexity of social organization, suggesting that selection pressure by pathogens could be a driver of complex social organization.
In termites, there is evidence for an association between selection pressure from pathogens, ecology, and social complexity. Although all termites are eusocial, the degree of worker altruism differs between termite lineages and aligns with ecology (Korb et al., 2012). Species of most early branching lineages of the termite phylogeny share a similar ecology called "one piece nesting" or "wood-dwelling" life type (Abe, 1987;Korb, 2007). Wood-dwelling termites nest in a single piece of wood that serves as food and shelter. "Workers" of wooddwelling termites are developmentally totipotent immatures that become reproductives, and thus are sometimes considered "false" workers. In particular, species of the family Kalotermitidae display little brood care by false workers and form less complex societies, where the interactions are generally not altruistic but cooperative (Korb, 2007). Yet, there are differences in the degree of brood care between wood-dwelling species. Brood care appears to be malleable depending on nesting ecology (Korb et al., 2012). Brood care is negligible or absent in several dry wood termites such as Cryptotermes secundus (Kalotermitidae) that nest in sound dry wood. More intensive brood care (especially allogrooming) occurs in the dampwood termite (Archotermopsidae), Zootermopsis nevadensis, which nests in rotten, decaying wood (Korb et al., 2012). Brood care is also present in Zootermopsis angusticollis (Rosengaus and Traniello, 1993). Looking beyond termites having a wood-dwelling life type, brood care is more intensive in termite lineages that have altruistic workers with reduced developmental options (true workers) and complex social organization (Roisin and Korb, 2011) like Macrotermes. These species are central-place foragers (multiple-pieces nesting, see Abe, 1987;Korb, 2007) with nests separated from the foraging ground.
The increase in social complexity and changes in ecology seem to align with the degree of immune challenge in two ways. First, wood-dwelling termites are confined within their nests, and hence it seems reasonable to assume that they are not exposed to microbial challenges as frequently as foraging termites. Second, within the group of wood-dwelling species, the nests of termites that live in sound drywood can be assumed to be microbe-and pathogen-poor when compared to nests of species that live in decaying dampwood. Rosengaus et al. (2003) provide evidence for the latter if we assume that pathogen pressure can be extrapolated from the cultivable fungal and bacterial loads that these authors measured. In fact, the dampwood termite, Z. angusticollis that was investigated in Rosengaus et al. (2003) seems to have huge constitutive investments in immune defense at different levels. This includes the individual as well as the colony level, potentially with socially acquired immunity (Rosengaus et al., 1999(Rosengaus et al., , 2011Traniello et al., 2002;Cole and Rosengaus, 2019). The investment in immunity in dampwood termites of the genus Zootermopsis is also visible at the genome level. In the genome of Z. nevadensis, six copies of Gram-negative binding proteins (GNBPs) were found, more than in many other insect species (Terrapon et al., 2014). GNBPs can serve as microbial detectors and effectors alike. Four of these presumably termite-specific genes were also found in the genome of the fungus-growing termite Macrotermes natalensis (Termitidae, see Poulsen et al., 2014;Korb et al., 2015). M. natalensis is a foraging species with intensive brood care and a complex social organization. For other central-place foraging termites, such as several Australian Nasutitermes species (Termitidae), which also have a complex social organization, as well as for Reticulitermes, GNBPs had previously been shown to be under positive selection (Bulmer and Crozier, 2006). Additionally, selection on three immune genes (GNBP1, GNBP2, and Relish) differed between Australian Nasutitermitinae. The rate of adaptive evolution of GNBP2 and Relish were increased during the transition from feeding on dry grass stored in epigeal nests to feeding on decaying wood (Bulmer and Crozier, 2006;Rosengaus et al., 2011), providing additional evidence for immune challenges that are specific to species that live in decaying wood.
Based on these results, we hypothesized that selective pressure on immune defense genes (IGs) differs across termites depending on their life type (Korb et al., 2015). We tested the hypothesis that wood-dwelling termites, which do not leave the nest to forage outside show relaxed selection on IGs and fewer signs of positive selection compared to soil-foraging species. Additionally, we tested whether within the wood-dwelling life type, the dampwood termite Z. nevadensis had stronger signs of selection than other wood-dwellers that nest in sound wood.
In order to test for differences in selective forces acting on IGs between wood-dwelling and foraging species, we analyzed a set of 81 previously identified single-copy IGs (see section "Materials and Methods") in eight termite species (published data see Misof et al., 2014;Poulsen et al., 2014;Terrapon et al., 2014;Harrison et al., 2018;Evangelista et al., 2019; data sources are provided in Supplementary Table S1). Four of the species are wood-dwelling species: C. secundus, Incisitermes marginipennis, Prorhinotermes simplex, and Z. nevadensis. The remaining four species are foraging species: Mastotermes darwiniensis, Reticulitermes santonensis (i.e., Reticulitermes flavipes), Coptotermes sp. and M. natalensis. These were analyzed in a phylogenetic framework of 22 species (one mayfly, 12 polyneopteran insects including above listed termites and their closest relatives Cryptocercidae, two paraneopteran, and six holometabolous insects) spanning winged insects to increase the statistical power of lineage specific tests for selection.

Phylogenetic Relationships
The species tree was inferred from a supermatrix including 1,178 single-copy orthologs (SCOs) and spanning an alignment length of 555,906 amino acid positions (partition coverage 100%, sitecompleteness score C a = 74.61%, see Supplementary Material and Supplementary Figure S1). Termites were monophyletic with Cryptocercus as sister group, consistent with earlier work (e.g., Lo et al., 2000;Klass and Meier, 2006;Inward et al., 2007;Legendre et al., 2008). Phylogenetic relationships within termites are largely consistent with Evangelista et al. (2019) and are statistically maximally supported (Figure 1). Consistent with earlier work (e.g., Legendre et al., 2008), neither wood-dwellers nor foragers constitute monophyletic groups, confirming that several independent switches in life type were included in our analyses. More details on phylogenetic analyses are provided in the Supplementary Material.

Patterns of Selection on Termite Immune Genes
Between 13 and 78 SCOs of IGs per species were included in the analyses (Table 1). We found no evidence for positive selection on the IGs ( Table 2 and Supplementary Table S2).
Next, we tested the hypothesis that selection on the IGs of wood-dwelling species is relaxed compared to foragers. We found 47 cases of significantly relaxed selection across all termite species analyzed (Table 1 and Supplementary Table S3, P < 0.05, FDR < 0.2). There was no evidence for a difference in the number of IGs under relaxed selection between the life types (generalized linear mixed effects model with binomial error distribution: df = 7, z = 0.096, P = 0.92). There was also no evidence for differences between species (generalized linear model assuming binomial error distribution: df = 7, z = −1.75-0, P = 0.08−1, ranges of z and P are for the different species). Because changes in the selection intensity on IGs could be obscured by genome-wide differences in selective constraint, it is important to test these hypotheses against the genomic background. To this end, we generated sets of background genes (BGs) that consisted of genes matching the GC-content and sequence length of each IG for each species (see section "Materials and Methods"). The number of IGs under relaxed selection did not differ significantly from that expected from the analysis of the BGs for any of the species [see 95% confidence intervals (CIs) for BG sets in Table 1]. In order to take genome-wide effects of selective constraint into account, when comparing selection intensity between wood-dwellers and foragers, (i) we calculated the ratio of genes under significantly relaxed selection between wood-dwellers and foragers for IGs and (ii) compared this ratio to that expected from BGs (see section "Materials and Methods"). If selection is relaxed specifically in the IGs of wood-dwellers relative to their genomic background, we expect the ratio of the number of significant genes under relaxed selection between wood-dwellers and foragers to be larger for the IGs than for the BGs. However, the ratio of the number of IGs under relaxed selection between wood-dwellers and foragers did not differ significantly from the expectation derived from BGs (Figure 2A), supporting the view that patterns of relaxed selection on IGs follow genome-wide trends in selective constraint.
Because we did not find any evidence for an increase in the number of IGs under significantly relaxed selection in wooddwellers, we reasoned that a putative signal of relaxation of selective constraint might be more diffuse and only become visible as a general trend over all IGs investigated. In order to capture such more general trends, we assessed potential differences in k, a measure for the intensity of selection, for all IGs between life types. k did not differ significantly between species (Kruskal-Wallis test: df = 7, X 2 = 3.95, n = 322, P = 0.79, for n per species see Table 1) nor was it lower for wooddwellers (Mann-Whitney U-test, one-sided: U = 11,690, n = 322, P = 0.41). This indicated similar selection intensity on IGs for wood-dwellers and foragers. For the comparison of k between life types it is, as above, important to take the selective constraint on the genomic background into account. k for the IGs did not differ significantly from k for sets of BGs for any of the species investigated (Table 1). Following the same rationale as above, we used the ratio of medians of k between wood-dwellers and foragers as a test statistic. This ratio can be interpreted as the relative intensity of selection between wood-dwellers and foragers. We found that the relative intensity of selection between wood-dwellers and foragers on IGs matched that of the BG sets ( Figure 2B), again suggesting that the IGs follow genome-wide trends of selective constraint.
Finally, we hypothesized that our results might be affected by the particular selection pressures that act on Z. nevadensis, which is a wood-dwelling species, but lives in dampwood nests. Nests in dampwood have a high microbial loads, as has been shown for Z. angusticollis (Rosengaus et al., 2003). Hence, selection would be expected to be stronger on Z. nevadensis IGs than on IGs in the other wood-dwellers, resulting in a smaller fraction of genes under significantly relaxed selection in Z. nevadensis. We found no evidence for this hypothesis (generalized linear model assuming binomial error distribution: df = 3, z = 1.73, P = 0.084). The overall intensity of selection on IGs (k) also did not differ significantly between Z. nevadensis and the other wooddwellers (Mann-Whitney U-test, U = 5,471, n = 215, P = 0.77). Similarly, it could be argued that the assumption of relaxed selection only holds for the dry wood-dwellers C. secundus and I. marginipennis, assuming that only dry wood is a truly pathogen poor substrate. We could not find a difference in the number of genes under relaxed selection between the dry wood-dwellers and the other species (generalized linear model with binomial error distribution: df = 7, z = 0.048, P = 0.96) nor for the intensity of selection over all IGs as measured by k (Mann-Whitney U-test: U = 10,538, n = 322, P = 0.37). FIGURE 1 | Phylogenetic relationships. Best ML species tree (out of 50 trees) inferred from a super alignment of 1,178 single-copy orthologs (SCOs) with 555,906 amino acid positions (schematized for Holometabola and Paraneoptera, all ML trees had the same topology). Statistical bootstrap support was derived from 200 replicates. The tree was rooted with the mayfly Ephemera danica (the full tree is provided with the Supplementary Files on DRYAD). Color code: brown indicates wood-dwelling species (wd), green indicates foraging species (f). TR, derived from transcriptome assemblies (Misof et al., 2014;Evangelista et al., 2019); OGS, derived from official gene sets (Supplementary Table S1). The symbol " * " indicates reference species whose OGS was used to create the ortholog set. Pictograms were kindly provided by H. Pohl, Jena. Genes that are under significantly relaxed selection were counted (k < 1, P < 0.05, FDR < 0.2) for columns containing #. BG, background single-copy ortholog. The symbol "*" indicates species with annotated genomes.
To our surprise, we observed that median k for IGs was smaller than one for seven of the eight investigated termite species (Table 1), indicating an overall relaxation of selection (Mann-Whitney U-test: U = 20,958, n = 322, P < 0.01). This signal was not IG specific: k for the sets of BGs was also significantly smaller than one for all species (P = 4 × 10e−18-4.8 × 10e−5), suggesting genome-wide relaxation of selection on the termite lineages compared to the background branches of the phylogeny.

DISCUSSION
In this study we combined recent genomic and transcriptomic resources to test whether termite ecology, in particular exposure to pathogens, might affect the evolution of immune genes. Surprisingly, and in contrast to studies from Drosophila Sackton et al., 2007;Hill et al., 2019), we could not find evidence for positive selection on the IGs for eight termite species. We extended our analyses, employing recently developed tests to explicitly assess relaxed selection , revealing 47 cases of significantly relaxed selection in IGs. We expected that the intensity of selection would differ between termites of the wood-dwelling and foraging life types, as foraging termites are assumed to experience higher selection pressure from pathogens due to higher exposure. Contrary to our expectation, we did not detect an effect of life type on signs of selection for immune genes. Neither did wooddwelling species differ from the soil foraging species, nor did the dampwood termite Z. nevadensis differ from the other wood-dwelling termites. A possible explanation may be that IGs that occur in multiple copies in at least one of the species that we analyzed contribute to adaptation to selection pressure from pathogens. Such multi-copy genes were excluded from our analyses. Furthermore, recent selective sweeps, as described for termicin in Reticulitermes (Bulmer et al., 2010), are difficult to detect with our methodology. For the detection of recent selective events, comprehensive polymorphism data would be required. Also, social mechanisms such as the exclusion of infested individuals and the impregnation of the nest walls with fungicidal compounds may protect efficiently against pathogens entering a colony or infecting individuals, thus buffering selection pressure on IGs (Traniello et al., 2002;Cremer et al., 2007;Bulmer et al., 2009;Rosengaus et al., 2011;Masri and Cremer, 2014). Another possible explanation for similar selection pressure on IGs between wood-dwelling and foraging species is that both harbor complex gut microbial communities that are essential for termite survival (Waidele et al., 2017), and perform, in principle, similar functions in lignocellulose digestion and nitrogen acquisition across several of the species that we analyzed (Waidele et al., 2019). It seems reasonable to assume that the immune system is likely to play a role in modulating these communities in wood-dwellers and foragers alike resulting in constant selection pressure. Our results support the results of a study in ants (Roux et al., 2014) that spanned a similar evolutionary time frame and number of focal species. These authors also found no evidence for a relaxation of selective constraint specific to IGs that could be related to social immunity. However, the authors found several instances of positively selected genes showing that it is in principle possible to detect positive selection with dn/ds based methods in a similar setup. Harrison et al. (2018) were able to pinpoint several selected codons using genome sequences from Blattella germanica and C. secundus, again spanning similar divergence times. In general, BUSTED is more powerful than other methods to detect positive selection because it has decent power to detect not only pervasive, but also episodic selection, while being fairly independent from evolutionary divergence times .
We applied state-of-the-art methods (Pond et al., 2005;Murrell et al., 2015;Wertheim et al., 2015) and carefully controlled our study by analyzing BGs that matched IGs in GC-content and sequence length (see section "Materials and Methods"). We also combined, to our knowledge, the largest data set of termite genomic and transcriptomic resources for our study that has been used in that context so far. Nonetheless, our study might lack statistical power: First, we have relatively few species for life type comparisons (4 vs. 4). Second, for five of these species (M. darwiniensis, I. marginipennis, R. santonensis, P. simulans, and Coptotermes sp.) only transcriptome data are available. Naturally, the number of genes under significantly relaxed selection is lower for the species for which only transcriptome data are available, as fewer genes could be annotated and investigated in the transcriptome data (Tables 1, 3  and Supplementary Table S4). However, potential biases in data availability were taken into account in the procedure to sample sets of matching BGs. Note that we did not apply the groupbased approach of the BUSTED and RELAX tests that might have increased power to detect selection or its relaxation. This was due to fragmented data coverage of the species representing both life types. This fragmented coverage for many genes would have led to only a limited gain in power. Furthermore, the careful selection of BGs would have become infeasible because too few BGs matched the species composition in the group-based tests on IGs. We are hopeful that the group-based approach will become more powerful in the future, when more termite genome data become available. Nonetheless, we were able to detect 47 instances of relaxed selection, showing that the RELAX approach can be powerful in a setup like ours. Because RELAX specifically tests for relaxed selection, while other studies often infer potentially relaxed selection indirectly by faster evolution in the absence of positive selection (Roux et al., 2014;Partha et al., 2017), we think it should be the preferred method. When taking the genomic background and the sampling procedure into account (Figure 2), our data provide no evidence that selection intensity on SCO IGs differs between life types. Do our results imply that selection pressure on immune genes is the same between termite life types? No, we cannot conclude this as our study excluded IGs that are not SCOs. For example, several GNBPs were excluded from our study because they seemed to be present in multiple copies in at least one termite species that we analyzed (Znev_03257, Znev_03259). Different copies can be caste-specifically expressed in termites as shown for Z. nevadensis (Terrapon et al., 2014), and for Reticulitermes speratus (Mitaka et al., 2017), and they may be under positive selection as indicated by a study on GNBPs for several foraging Nasutitermes species (Bulmer and Crozier, 2006). Note that the original hypothesis for a difference in immune defense between wood-dwelling and foraging termites considered specifically GNBPs and AMPs (Korb et al., 2015). Other GNBPs were excluded because they were restricted to only some of the lineages that we based our ortholog set on (Znev_03260, Znev_02878, and Znev_00933). Thus, more studies are warranted that test selection on genes that are lineage-specific or might occur in multiple copies. However, orthology is difficult to infer for multi-copy genes making the restriction to SCOs a standard procedure (e.g., Dowling et al., 2016;Pauli et al., 2016;Mitterboeck et al., 2017;Ran et al., 2018;Brandt et al., 2019;Hill et al., 2019). Thus, we restricted our analysis to SCOs because orthology is a basic assumption of the current methods that identify selection in a powerful phylogenetic framework [e.g., codeML implemented in PAML, Yang, 2007; HyPhy applying BUSTED, see Murrell et al. (2015), and RELAX, see Wertheim et al. (2015)].
Applying state-of-the-art methods to a comprehensive termite data set, for which genome or transcriptome data are currently available, we found no evidence for differences in selection on immune genes that correlate with termite life type. Our results suggest that the putative evolutionary response to differences in expected pathogen exposure can not be found in single-copy immune genes. Interestingly, we detected a signal of genomewide relaxation of selective constraint in termites. We speculate that this could be related to their social organization that might lead to smaller effective population size (Romiguier et al., 2014;Rubenstein et al., 2019) because only the kings and queens reproduce, and hence contribute to the effective population size. In smaller populations, natural selection becomes less effective at purging deleterious mutations as well as at driving advantageous mutations to fixation (Ohta, 1973). This is equivalent to a relaxation of selection in smaller populations. Thus, small effective population sizes compared to the other insects in our study could have manifested as the genome-wide signal of relaxed selection that we observed in termites.

MATERIALS AND METHODS
A comprehensive diagram summarized major analysis steps in Supplementary Figure S1.

Identifying Orthologous Sequence Groups of Protein-Coding Single-Copy Genes
As basis to identify SCOs, we designed an ortholog set from official gene sets (OGS) from available full (draft) genomes of four species: C. secundus, M. natalensis, Z. nevadensis, and B. germanica (Supplementary Table S1). The set of SCOs was inferred with the software OrthoFinder v.1.1.4 (Emms and Kelly, 2015) using default settings. As input, the OGS of respective species were downloaded from public databases as amino acid and nucleotide sequences. The OGS of C. secundus was kindly provided by the C. secundus consortium (via J. Korb) before it was published (Harrison et al., 2018). We only kept the longest isoform per orthologous group (OG). All OGs that included the amino acid Selenocysteine (U) were removed to avoid difficulties in downstream analyses as many software packages are not able to handle Selenocysteine. This was done using the package BioBundle [script isoformCleaner with boost 1.61.0 environment (Kemena, 2017, available from github 1 )]. SCOs inferred with OrthoFinder, were summarized with custom-made Python scripts (kindly provided by A. Faddeeva and L. Wissler, available upon request). This resulted in a set of 5,382 SCOs across the four selected species.

Taxon Sampling
We included the four reference species that were used to create the ortholog set as well as genome and transcriptome data of 18 additional species in our analyses. Eight of the included species are termites: Coptotermes sp., I. marginipennis, M. darwiniensis, P. simplex, and R. santonensis (with published transcriptome data); C. secundus, Z. nevadensis, and M. natalensis (with published OGS), see Evangelista et al. (2019). Other included species were a representative of Cryptocercidae as it is supposed to be the sister group of termites (e.g., Lo et al., 2000;Inward et al., 2007), two other non-social cockroach species, and representatives from other polyneopteran, paraneopteran and holometabolous insects, and a mayfly as outgroup (Adams et al., 2000;Xia et al., 2004;Sinkins, 2007;Tribolium Genome Sequencing Consortium et al., 2008;Bonasio et al., 2010;International Aphid Genomics Consortium, 2010;Elsik et al., 2014;Misof et al., 2014;Poulsen et al., 2014;Terrapon et al., 2014;Wang et al., 2014;Mesquita et al., 2015;Pauli et al., 2016;Harrison et al., 2018;Evangelista et al., 2019; the full species list is provided in Supplementary Table S1). Access to transcriptome data (see Figure 1) was kindly granted by 1KITE before they were published, access to the OGS of the locust and the mayfly was granted by the i5K community.

Assignment of Putative Orthologous Transcripts to the SCOs
The ortholog set was used as input for the assignment of putative SCOs (provided as Supplementary Files on DRYAD, doi: 10.5061/dryad.j6q573n98). Inference and assignment of putative orthologs from genome and transcriptome data of the 18 species that were not included for generating the ortholog set was performed with OrthoGraph v.0.6.1 (Petersen et al., 2017). OrthoGraph is recommended to infer orthologs from transcriptome data for which no OGS are available (see Petersen et al., 2017). OrthoGraph analyses resulted in 5,366 SCOs that were identified in at least one species that was not used as reference species to create the ortholog set.

Multiple Sequence Alignments, Species Tree Inference and Testing for Selection
Individual SCOs were aligned at the amino acid level with MAFFT v7.310 using the L-INS-i algorithm (Katoh and Standley, 2013).

Inferring Natural Selection Alignment processing and clean-up
Methods to identify selection are sensitive to misalignments (Markova-Raina and Petrov, 2011;Privman et al., 2012). Therefore we performed extensive alignment clean-up. First, we identified and deleted badly aligned or gap-rich sequences on amino acid level with MaxAlign v1.1 (Gouveia-Oliveira et al., 2007). This procedure resulted in five SCOs with only one sequence which were excluded from further analyses. We subsequently compiled corresponding nucleotide (i.e., codon) MSAs with PAL2NAL (Suyama et al., 2006, v14.1, see Misof et al., 2014 using the 5,361 amino acid MSAs as blue-print. The nucleotide MSAs were then used for all following analyses. Second, we deleted all SCOs with less than four sequences (223 SCOs) leaving 5,138 SCOs. Third, we identified ambiguously aligned sections on amino acid level with Aliscore v2.2 (Misof and Misof, 2009;Kück et al., 2010) with the same settings as described for the species tree inference. Suggested sections were removed from the amino acid and correspondingly from the nucleotide MSAs with AliCUT v2.3 (Kück, 2011). Subsequent analyses were performed on the masked nucleotide MSAs. First, we classified 5,138 SCOs into 86 immune single-copy genes (IGs) and into the remaining 5,052 SCOs based on Supplementary Table S25 from Terrapon et al. (2014, for Z. nevadensis) and Korb et al. (2015, for Z. nevadensis and M. natalensis), see Supplementary Table S5. The 5,052 nonimmune SCOs were used to generate gene sets from the genomic background, i.e., BGs that had similar GC-content and sequence length (see below) as the examined IGs. Note that from the 86 IGs (Supplementary Table S5) five IGs were excluded because there was no SCO fulfilling the criteria to serve as BG and these were not listed by Terrapon et al. (2014) or not reported by Korb et al. (2015). This left 81 IGs for analyses (detailed information are provided in the Supplementary Material).
To further reduce potential false positives that may originate from misalignments, we trimmed trailing ends of each MSA, i.e., each MSA started and ended with unambiguous nucleotides for all species. Because visual inspection of the trimmed MSAs still revealed putative misaligned nucleotides, we applied the GUIDe tree based AligNment ConfidencE approach (GUIDANCE) Guidance2 (Landan and Graur, 2008;Sela et al., 2015) version 2.02 using MAFFT as implemented alignment method on the trimmed MSAs (options: codon as sequence type, sequence cutoff = 0 and the default column cutoff = 0.93).

Inferring positive selection and selection intensity
To test for evidence of positive selection we used BUSTED  as implemented in the software package HyPhy (Pond et al., 2005). BUSTED uses a branch-site test for positive selection on entire genes in a foreground branch relative to the background branches in a phylogeny. A significant Pvalue means that at least one codon in the foreground branch has experienced at least an episode of positive selection. The high sensitivity of the method compared to tests from alternative packages (see e.g., Enard et al., 2016;Ebel et al., 2017;Venkat et al., 2018;Hill et al., 2019) and the option to define the foreground branches according to our research question made it perfectly suited for our study.
For inferring potential relaxation of selection, we used RELAX  as implemented in the software package HyPhy (Pond et al., 2005). RELAX has been designed to identify changes in the intensity of selection on a given protein-coding gene in a codon-based phylogenetic framework (see Wertheim et al., 2015). The basic expectation of RELAX is that under relaxed selection, the ω of sites under purifying and positive selection will move closer to neutrality. The change of ω for the selected sites relative to the background branches is quantified with the selection intensity parameter k, where If parameter k is significantly larger than one, selection has been intensified along the test branches. If k is significantly smaller than one, selection has been relaxed.
We used BUSTED and RELAX as implemented in the software package HyPhy, version 2.4.0-alpha.2 (access: April, 2019). We performed BUSTED and RELAX for each gene in each termite species separately (focal species: foreground, remaining species in the alignment: background) and calculated false discovery rates (FDR) to correct for multiple testing. We then determined k for each of the termite species relative to all other species in the tree. We chose this species-wise analysis setup because we wanted to take species level differences in potential pathogen exposure into consideration. For example, Z. nevadensis that resides in dampwood might differ in microbial exposure from Cryptotermes and Incisitermes that reside in dry wood, which in turn might affect selection pressure on IGs. Furthermore, an effective selection of BGs was only feasible in the species-wise framework (see section "Discussion"). Results from the specieswise analyses were summarized by life type after the BUSTED and RELAX analyses.

Comparison of IG selection parameter to the genomic background
To test whether or not signals of selection were specific to IGs or the consequence of genome-wide trends, we generated sets of BGs. To this end, we searched the 5,052 non-immune SCOs for genes that closely matched the GC-content (±5%) and the sequence length (±5%). Following this procedure, we generated lists of matching BGs for each IG and each species. From these lists, we randomly sampled 100 gene sets such that there was a matching BG for each IG that was analyzed in the respective species (e.g., for Z. nevadensis, 78 IGs were analyzed, thus each of the 100 sets of BGs contained therefore 78 BGs, see also Table 3. Lists of analyses BG that are similar in GC-content and sequence length of the IGs are provided for each species as Supplementary Files on DRYAD). We performed BUSTED and RELAX analyses for each termite species on all IGs and RELAX on the species' respective BG set with default settings. The same cutoffs as for the IGs (P < 0.05, FDR < 0.2, k < 1) were applied to the BGs.
All analyses were performed on Linux Desktop PCs at the University of Freiburg, Germany and on the Linux HPC CSIRO Cluster Pearcey, Australia. Analyses results of all IGs are summarized in Supplementary Tables S2, S3; results of BGs are provided species-wise on DRYAD.

Statistical Analyses
In order to assess potential differences in selection intensity on the IGs between species and life types, we summarized the RELAX results by counting the number of genes under significantly relaxed selection: genes with k < 1, P < 0.05, and FDR < 0.2 were considered. With our FDR cut-off, we followed the recommendation from Efron (2007) for genomewide analyses. Potential differences were tested for statistical significance with generalized linear models with binomial error distribution using the functions glm and glmer from the lme4 R package [Bates et al., 2015, version 1.1-21 with R Core Team (2018), version 3.4.4]. The number of significant genes divided by the total number of genes analyzed was used as response variable. Species or life type were used as potential predictors. Varying sampling depths between species, as represented by the number of IGs analyzed per species, were taken into account as weights in the model. When comparing life types, species were treated as a random effect. See Supplementary File (RanalysisscriptfortermiteIGs.R on DRYAD, doi: 10.5061/ dryad.j6q573n98) for a detailed R analysis script with all models, commands and functions used.
We also analyzed parameter k to search for more diffuse trends in selection intensity that are distributed over the IGs so that individual IGs do not reach significance. According to its definition, k should map linearly on a logarithmic scale. However, we found six strong outliers on the logarithmic scale that were more than three standard deviations away from the mean [log (k) < −9, see Supplementary Table S3] that could make the analysis in a linear framework error prone. Visual inspection of the alignments underlying these extremely small values of k did not reveal any obvious misalignments that would justify their exclusion. Therefore, potential differences in k were assessed with non-parametric tests (Mann-Whitney U-test, Kruskal-Wallis test) that are robust to outliers.
Genome-wide trends in selection intensity can potentially obscure IG specific patterns or generate false positives. For example, changes in population size can affect the efficiency of both purifying and positive selection (Ohta, 1973) on a genome-wide scale. Population sizes might differ between species and life types in our study depending on reproductive rates and degrees of sociality. Therefore, it is essential to put the results for IGs into the context of the genomic background. To this end, we generated expected values for the number of significant genes and for k based on 100 sets of BGs (see above) per termite species, representing the genomic background. The median of k and 95% CIs from the BGbased distributions for each species were calculated with R (version 3.4.4), using the median and quantile functions with standard settings. In order to compare differences between life types while taking the genomic background into account, we calculated (i) the ratio of the number of genes that were significantly relaxed between wood-dwellers and foragers R relax = # genes wood-dwellers | P < 0.05 ∩ FDR < 0.2 ∩ k < 1 # genes foragers | P < 0.05 ∩ FDR < 0.2 ∩ k < 1 and (ii) the ratio of median parameter k (Rk, relative selection intensity): Rk =k wood-dwellers k foragers . These ratios were calculated for the IGs and the BGs. Then the ratio of the IGs was compared to their expectation from the BGs. Significant shifts in selection intensity that are specific to IGs should lead to shifts of R relax and R ki only for IGs. Thus, if there were IG specific patterns of relaxed selection, the ratios R relax and R ki for the IGs should represent extremes of the distribution of sets of BGs. Therefore, we only considered a signal as significantly specific for IGs if our test-statistics of the IGs were outside of the 95% CI calculated for the BGs. However, this was not the case in our study.

DATA AVAILABILITY STATEMENT
Supplementary Data (see Supplementary Material) used in this study can be found at the DRYAD digital repository, doi: 10.5061/ dryad.j6q573n98.

AUTHOR CONTRIBUTIONS
JK and FS conceived the study. KM, MS, and FS performed all analyses. FS, KM, and JK wrote the manuscript. All authors edited and approved the manuscript.

FUNDING
This project was supported by the German Science Foundation (DFG; KO1895/20-1, KO1895/20-2, STA 1154/2-1 Projektnummer 270882710). The article processing charge was funded by the German Research Foundation (DFG) and the University of Freiburg in the funding programme Open Access Publishing. Funders had no role in research design or decision to submit.

ACKNOWLEDGMENTS
We acknowledge Bernhard Misof, Panagiotis Provataris, Coby Schal, Xavier Belles, Stephen Richards and the i5K community the usage of the official gene set of Ephemera danica, and Blattella germanica, and Xianhui Wang and Le Kang for access and usage of the official gene set of Locusta migratoria. We thank the 1KITE Dictyoptera group who granted us access to transcriptome assemblies before they were published. We thank David Enard (University of Arizona, United States), Dario Valenzano (MPI, Cologne, Germany), Ondrej Hlinka (CSIRO, Australia), Ryan Velazquez, and Sergej Pond for their useful input and help with analyses of the molecular sequence data using Guidance2 and HyPhy. KM thanks Thomas Pauli (University of Freiburg, Germany) for fruitful discussions and Hans Pohl (University of Jena, Germany), who kindly allowed the usage of several pictograms in Figure 1. We also thank the two reviewers for helpful comments on the manuscript.