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 ).