Introduction
Parasitoid attack always leads to the death of either host or parasitoid, and genes affecting the host capacity to defend themselves against the intruder should therefore evolve rapidly (Kraaijeveld, Van Alphen, & Godfray, 1998). In insects being attacked by endoparasitoid wasps, which inject the egg inside the host body, host survival depends on an immunological reaction that ends in the encapsulation of the enemy egg (Carton, Poirié, & Nappi, 2008). Hence, a major task in evolutionary studies of host-parasitoid systems is to identify the genes underlying phenotypic differences in immunocompetence (Wertheim, 2015). Work on model organisms, such as Drosophila melanogaster , and a range of non-model organisms have identified key pathways involved in the immunological process (Carton et al., 2008). We also know that host species or genotypes may differ greatly in their defence capacity (Carton et al., 2008; Fors, Markus, Theopold, & Hambäck, 2014; Wertheim, 2022), but we lack information outside of theDrosophila system on the importance of canonical immune genes as determinants of defence phenotypes.
In our work, we have observed large phenotypic differences between three closely related leaf beetle species (Chrysomelidae: Galerucellaspp.) in defense capacity when attacked by the same parasitoid wasp species (Asecodes parviclava , Eulophidae). WhereasGalerucella pusilla shows a high capacity to encapsulate wasp eggs, this capacity is much lower in G. tenella and almost absent in G. calmariensis (Fors et al., 2014). These defense phenotypes have also been connected to differences in the induction of those immune cells involved in the encapsulation of wasp eggs (Fors et al., 2014) and in the expression of immune-related genes following parasitoid attack (Yang et al., 2020), suggesting a genetic basis for the observed differences in immune performance. Importantly, all three of these closely related species are attacked by the same parasitoid, suggesting that these large differences in immunocapacity evolved recently. Thus, this system provides a promising opportunity outside ofDrosophila to connect defense phenotypes and genotypic differences by identifying genes under positive selection. Within this top-down framework, here we ask whether species-level differences in phenotypes with clear connections to Darwinian fitness in response to a common, strong selective pressure are mirrored in molecular tests of selection, with the expectation that parasitoid-relevant immune genes in the most defended species should be enriched for positive selection.
The insect immune system consists of two parts, where the humoral system is mainly active against pathogens (Bulet, Hetru, Dimarcq, & Hoffmann, 1999; Gillespie, Kanost, & Trenczek, 1997), whereas the cellular immune system seems more important when encapsulating parasitoid wasp eggs (Carton et al., 2008). The cellular immunity can be separated into seven broad functional categories: recognition, signalling, effector, proteases, haematopoiesis, melanisation, and wound healing (see also Yang et al., 2020) and particularly genes 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). However, selection dynamics are complex and depend on specific functions and taxa studied (Keehnen, Hill, Nylin, & Wheat, 2018; Keehnen, Rolff, Theopold, & Wheat, 2017), but some classes show more evidence of positive selection while other classes such as AMPs appear to experience more balancing or purifying selection (Heger & Ponting, 2007; Unckless, Howick, & Lazzaro, 2016; Unckless & Lazzaro, 2016; Waterhouse et al., 2007). So far, however, studies focused upon the selection dynamics of parasitoid specific immune genes have received less attention. Here we make the prediction that species exhibiting the strongest immunocompetence against wasp attack (G. pusilla ) will also exhibit stronger positive selection on relevant immune genes compared with the species with lower immunocompetence, as we expect immunocompetence to have evolved via historical selection upon relevant genes.
Apart from defence phenotypes, there are also other phenotypic differences between the three species of unknown or absent relation to parasitoid defence capacity. First, G. calmariensis is slightly larger and develops faster than the other two species, indicating potential differences in development and metabolic pathways. Second, 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), which may be traced to differences in gene pathways related to pigmentation and melanisation (Ito et al., 2010; Linnen, O’Quin, Shackleford, Sears, & Lindstedt, 2018). Finally, the species likely 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. Thus, each species is expected to exhibit unique signatures of selection, highlighting the importance of our a priori hypothesis framework with regards to parasitoid-relevant immune genes.
To detect selection, we employ two molecular tests that are reliant upon population level genetic variation but primarily differ in their genomic focus (coding regions vs. genome wide), and the time scale over which they have the most statistical power. For coding regions, 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 nonsynonymous and synonymous polymorphism (Pn/Ps ) within species to the ratio of nonsynonymous and synonymous fixed differences (Dn/Ds ) between species (McDonald & Kreitman, 1991). When positive selection favours the phenotypic impact of a novel amino acid change, the corresponding advantageous allele goes to fixation, and when this happens repeatedly, 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 at genes can be estimated using α (= 1-DsPn/DnPs ) (Smith & Eyre-Walker, 2002), which represents the proportion of non-synonymous substitutions driven by positive selection during divergence between focal species and an outgroup, allowing for a powerful test of selection dynamics.
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 these limitations, we here 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.
For regulatory regions, we employed a genome wide scan of positive selection dynamics that leverages information in the frequency of genetic variation found within populations of a given species to detect both soft and hard selective sweeps (aka positive selection). Due to the reliance on the site frequency spectrum (SFS), this test covers a more recent time history of selection compared to the aforementioned MK test, yet allows for the detection of selection anywhere in the genome, such as in the regulatory region flanking genes. We investigated genome signatures of positive selection using LASSI Plus (Harris & DeGiorgio, 2020). This method estimates haplotype frequency spectrum statistics within sliding windows of population genomic data to detect soft and hard sweeps, with stronger power for the latter category due to the more dramatic signature left in the SFS (Harris, Garud, & DeGiorgio, 2018).
To test our primary hypothesis 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 ), we quantified selection dynamics in our study system using the two molecular tests of selection (HDMKPRF and LASSI plus) using whole genome re-sequencing data from 15 individuals from each of the three beetle species (G. calmariensis , G. pusilla and G. tenella ).