Introduction
Parasitoids are major natural enemies of many insects, leading to high mortality rates of host species (Asgari & Rivers, 2011; Hawkins, Cornell, & Hochberg, 1997). One particularly important group are koinobiont endoparasitoid wasps, which lay eggs and complete their entire larval development inside their hosts, leading to host death. To defend themselves, insect hosts respond immunologically within hours by encapsulating and killing the wasp eggs (Carton, Poirié, & Nappi, 2008). Due to the ever-changing virulence of wasps, host species have been shown to rapidly evolve novel defence mechanisms (Kraaijeveld, Van Alphen, & Godfray, 1998), and this host-parasitoid coevolutionary dynamic likely explains the large amount of variation in defence mechanisms between host populations and species (Carton et al., 2008; Wertheim, 2022). Therefore, traits associated with host defence against parasitoids are widely regarded to be under strong selection pressures with a potentially crucial role in species divergence (Carton et al., 2008). However, while many studies have investigated the physiological and behavioural responses upon parasitoid attack in single host species (De Roode & Lefèvre, 2012; Fors, Markus, Theopold, Ericson, & Hambäck, 2016; Godfray, 1994), and evolutionary responses of hosts to various pathogens (Kim-Jo, Gatti, & Poirié, 2019), few studies have compared divergence in immune genes among closely related host species that are attacked by a common parasitoid species except for wasps attackingDrosophila (Wertheim, 2022).
We have previously investigated host-parasitoid interactions in a system with three closely-related beetle species (Galerucella calmariensis , G. pusilla and G. tenella , Coleoptera: Chrysomelidae) that differ in their hemocyte composition, in their capacity to induce potent hemocytes upon attack and in their encapsulation success against the parasitic eulophid wasp Asecodes parviclava (Fors et al., 2016; Fors, Markus, Theopold, & Hambäck, 2014). At the cellular level, G. pusilla exerts a much more potent immune response, with strong encapsulation and melanisation responses against wasp eggs, compared to G. calmariensis where encapsulation events of wasp eggs in host larvae are rarely observed. The immune response of G. tenella , in turn, is intermediate between G. pusilla and G. calmariensis . Recently, these differences in immune responses among species were also confirmed by gene expression analysis, especially in showing differential expression of genes involved in signalling, haematopoiesis, and melanisation inG. pusilla (Yang et al., 2020).
Insects rely on their innate immune responses to defend themselves against parasitoid and pathogen attack, as they lack an adaptive immune response. This innate immune system consists of two parts, humoral and cellular immunity. The humoral defence system is mainly active against pathogens and is mediated by genes regulating the expression of antimicrobial peptides (AMPs), wound healing, coagulation, and melanisation pathways (Bulet, Hetru, Dimarcq, & Hoffmann, 1999; Gillespie, Kanost, & Trenczek, 1997). Immune responses against parasitoid wasps, on the other hand, are usually regulated by cellular immunity, which can be classified into seven broad functional categories: recognition, signalling, effector, proteases, haematopoiesis, melanisation, and wound healing (see also Yang et al., 2020). Immune genes, particularly those involved in the recognition phase, have been found to evolve faster than other genes (Nielsen et al., 2005; Sackton et al., 2007; Schlenke & Begun, 2003; Waterhouse et al., 2007), but selection dynamics are complex and depend upon specific functions and taxa studied (Keehnen, Hill, Nylin, & Wheat, 2018; Keehnen, Rolff, Theopold, & Wheat, 2017). Recent studies based on amino-acid substitutions in Drosophila species suggest that some classes of immune-related genes show more evidence of positive selection than other genes (Heger & Ponting, 2007; Waterhouse et al., 2007) while other functional classes such as AMPs appear to experience more balancing or purifying selection (Unckless, Howick, & Lazzaro, 2016; Unckless & Lazzaro, 2016). These findings are however restricted to the scope of anti-microbial immune genes, whereas the selection dynamics of genes involved in the host-parasitoid specific immune response have received less attention. The observed phenotypic differences between species in our study system suggest differences in historic and current rates of evolution in key immune-related genes involved in the defence against parasitoids, and makes the system highly suitable for studying host selection dynamics in response to parasitic wasps. In particular, we expect stronger positive selection on parasitoid related immune genes in the species with the strongest immunocompetence against parasitoid wasps (G. pusilla ) compared with the species with the least immunocompetence (G. calmariensis ).
Apart from immunity related genes, other gene groups are also expected to differ between the three species, though the link between such genes and parasitoid attack are less clear or absent. For instance, G. calmariensis are slightly larger in size and develop faster than the other two species, indicating potential differences in pathways concerning development and metabolic processes (Carreira, Mensch, & Fanara, 2009). Moreover, Galerucella larvae differ in colour; bright yellow in G. calmariensis , pale white-yellow in G. pusilla and yellow with many black/brown spots in G. tenella(Hambäck, 2004). Such colour differences and spotting patterns may involve gene pathways related to pigmentation and melanisation (Ito et al., 2010; Linnen, O’Quin, Shackleford, Sears, & Lindstedt, 2018), or master regulatory genes (Deshmukh, Baral, Gandhimathi, Kuwalekar, & Kunte, 2018). Finally, the species may also differ in mate finding traits, such as the occurrence and detection of pheromones and cuticular hydrocarbons, and host-plant finding traits, as G. tenella uses a different host plant than the two other species.
A widely-used approach for investigating positive selection based on DNA sequencing data is the McDonald–Kreitman test (MK test), which infers the direction of natural selection by comparing the ratio of non-synonymous and synonymous polymorphism (Pn/Ps ) within species to the ratio of synonymous and nonsynonymous divergence (Dn/Ds ) between species (McDonald & Kreitman, 1991). When positive selection favours the phenotypic impact of novel amino acid changes, and the corresponding advantageous alleles go to fixation, positive selection is expected to yieldDn/Ds >Pn/Ps , under the assumption that synonymous mutations are evolving neutrally and there is no change in constraint over time. In contrast, if weak purifying selection is prevalent, deleterious alleles can segregate in the population for extended periods, yet rarely fix and therefore contribute little to divergence. Under this scenarioDn/Ds <Pn/Ps . The contribution of positive selection to amino acid divergence can be estimated usingα (= 1-DsPn/DnPs ) (Smith & Eyre-Walker, 2002), which represents the proportion of substitutions driven by positive selection.
The classic MK test was designed in the pre-genomics era, analyzing each gene locus separately and usually limited to comparisons of pairs of species (McDonald & Kreitman, 1991). It was not designed to infer the directionality or relative strength of selection across 1000s of genes across several study species. To overcome the limitations, here we used the high-dimension McDonald-Kreitman Poisson random field method (hereafter HDMKPRF, Zhao et al., 2019). This method is an extension of the MKPRF method developed by Sawyer and Hartl (1992), applying a Bayesian model across multiple gene loci to simultaneously estimate population genetic parameters of multiple target species, including lineage specific mutation rates and relative effective population size (Ne ), which improves inference of directionality and relative strength of selection along the lineages unique to each species.
To examine selection dynamics in our study system, we generated whole genome re-sequencing data from 15 individuals from each of the three beetle species (G. calmariensis , G. pusilla and G. tenella ). We then investigated the selection dynamics acting on coding regions from single-copy orthologs using HDMKPRF as well as a set of candidate immune genes. Based on previous studies, our primary hypothesis was that immune genes involved in wasp attack would more frequently have experienced positive selection in the species with the strongest immune response (i.e. G. pusilla >G. tenella > G. calmariensis ).