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