4. Discussion
This study reports one of the first analyses comparing genome-wide
natural selection across multiple lineages of beetle species, and the
first analysis identifying host species differences in the selection
dynamics of putative immune genes against parasitoid wasps. Our study
system included three closely-related leaf beetle species
(Galerucella spp. ) with previously detected phenotypic
differences in their capacity to encapsulate eggs from the shared
parasitoid wasp (Asecodes parviclava ), as well as differences in
gene expression following wasp attack (Fors et al., 2016; Fors et al.,
2014; Yang et al., 2020). Using a modified MK-test based approach and
whole genome re-sequencing data from 44 individuals of the three
species, our analysis identified, as hypothesized, a higher number of
immune genes under positive selection in the species with the strongest
capacity to encapsulate wasp eggs (G. pusilla , N = 8).
Species lacking the capacity to encapsulate wasp eggs (), or having an
intermediate encapsulation capacity (G. tenella, N = 4), had a
lower number of immune genes under positive selection (G.
calmariensis, N = 1). Even though these genes have yet to be
functionally characterized in the three beetle species, previous
comparative gene expression analyses (Yang et al., 2020) as well as the
differential selection between beetles with high and low immune
phenotype in this paper suggest they are important in mediating defence
against parasitoid wasps. Additionally, a genome-wide scan based on
orthologous loci among the three species characterized different
selection patterns for a range of other genes, suggesting divergent
directions of selective pressure between species.
Even though we could analyze only a subset of immune genes in theGalerucella genome (72 out of 166), the detected genes indicate
the type of processes most likely under selection. Most importantly, the
analysis suggests that the observed phenotypic differences in the
capacity to encapsulate wasp eggs (Fors et al., 2014) are likely
explained by genes involved in wasp egg detection and in the recruitment
of hemocytes building the capsule (Yang et al., 2020), whereas genes
involved in melanisation or wound healing processes seem to be more
conserved. In the species with the strongest encapsulation capacity
(G. pusilla ), we identified two recognition genes
(santa-maria and Corin ), two signalling genes (Tak1and grass ), one protease coding gene (CG32483 ) and three
genes involved in haematopoiesis (Raf , cher andzfh1 ) that were all under positive selection. This set of eight
genes can be contrasted against the set of four immune genes that were
identified as under positive selection in G. tenella , the species
with intermediate encapsulation capacity, or the single positively
selected immune gene in G. calmariensis , the species with the
weakest encapsulation capacity. It is however notable that two genes
(Corin and cher ) were under positive selection in bothG. tenella and G. pusilla , whereas two other genes
(recognition gene PGRP-LE and haematopoiesis gene Cyt-b5 )
were only found as being under positive selection in G. tenella ,
indicating both parallel and diverging selection patterns. The
importance of these genes for the immune system in this and other
species is also evident as the two genes Corin and cherwere similarly differentiated among Drosophila species with
different encapsulation ability (Salazar-Jaramillo et al., 2014),
whereas three other genes (Grass , zfh1 and Cyt-b5 )
were differentially expressed following wasp infection involvingA. parviclava and G. pusilla (Yang et al., 2020).
The specific functions for these genes in Galerucella have not
been identified, but some studies exist from Drosophila . First,
both santa-maria and Corin have been suggested to encode
proteins with scavenger receptor activities, whereas PGRP-LEencodes an intracellular protein that binds to diaminopimelic acid-type
peptidoglycans to activate the IMD/Relish pathway (Bosco-Drayon et al.,
2012). Second, both Tak1 and Grass have been predicted to
be involved in the Toll signalling pathway, which may be the most
important signalling pathway in the immune response against wasp attack
in Drosophila (Carton et al., 2008; Lemaitre & Hoffmann, 2007;
Wertheim et al., 2005). Third, the protease coding gene CG32483is predicted to enable serine-type carboxypeptidase activity but there
is no known defense function against wasp attack. Finally, the four
haematopoiesis genes are widely known to be involved in the regulation
of lamellocyte differentiation and proliferation, which is the hemocyte
type that makes up most of the capsule around the wasp egg in bothGalerucella (Fors et al., 2014) and Drosophila (Kim-Jo et
al., 2019). Among these genes, Raf is suggested to play a crucial
role in the transition from pro-hemocytes to lamellocytes (Luo, Rose,
Roberts, & Dearolf, 2002), cher is predicted to be involved in
the negative regulation of lamellocyte differentiation (Rus et al.,
2006), zfh1 is one element in a transcription factor cascade that
plays a role as the switch between plasmatocyte and lamellocyte fate,
and Cyt-b5 encodes a conserved hemoprotein that is required for
hemocyte regulation (Kleinhesselink, Conway, Sholer, Huang, & Kimbrell,
2011). The specific immune genes selected for likely vary between
species (Wertheim, 2022), but our results indicates that phenotypic
differences in the immune response in these species is related to
multiple functional genes.
When examining other gene functions identified as experiencing positive
selection, we found several interesting candidates, where genes coding
for imaginal disc pattern formation and wing disc pattern formation were
under positive selection in all species. Imaginal discs are epithelial
sacs found in insect larva that later develop into cuticular structures
(e.g., head, wing, limbs, thorax) of adult insects, and the wing disc is
among the largest imaginal discs in insects (Blair, 2009). These
findings are not surprising as the beetle species are morphologically
differentiated, both in size and colour. Other genes varied between
species, but those under positive selection in G. calmariensisspecifically involved genes coding for metabolic functions such as the
metabolism of carbohydrate derivatives, oligosaccharides, amino sugars,
sulphur compounds and catechol-containing compounds. These sets of
positively selected gene functions indicate the importance of energy
allocation and dietary transitions during evolution in this species.
Moreover, pathways related to the nervous systems were also found to be
positively selected in G. calmariensis . Positive selection on
nervous system-related genes were previously documented in social
insects such as bees and ants (Roux et al., 2014; Woodard et al., 2011)
but have rarely been reported in beetles. These patterns may indicate
changes either in the capacity to detect host plants or mates and may
thus be involved in the species differentiation. In G. tenella ,
several pathways related to lipid metabolic processes were enriched in
positively selected genes. Some of these processes, such as sterol
metabolism, may be linked to needs to handle differences in plant
chemistry from the different host plants (Rosaceae vs. L.
salicaria ). Other lipids, such as sphingolipids, have been suggested to
be involved in cell defences and could be interesting to study in
relation to wasp attack. Finally, a unique pathway under positive
selection in G. tenella was involved in the pigmentation
metabolic process, which may potentially explain the different spotting
and pigmentation patterns in the species.
Some methodological concerns are relevant in relation to our study.
First, we were only able to include 4,154 orthologous genes in the
analysis, which may have reduced the power to detect important GO terms.
However, the removals should not have biased our conclusions as genes
were likely removed randomly across the genome. Moreover, by discarding
genes of low coverage and removing bias due to recent gene birth death
dynamics, this smaller set of genes provides a robust gene set for
inferring selection dynamics and is sufficient for estimating genome
wide selection dynamics. Second, a general problem in MK tests is their
sensitivity to the demographic history (Eyre-Walker, 2002; McDonald &
Kreitman, 1991). In our analysis, we approached this problem by using
HDMKPRF, which attempts to correct for demography effects during the
detection of selection patterns by accounting for the changes in
effective population size and other population genetic parameters. One
potential concern, however, is that the proportion of substitutions
driven by positive selection were highest in G. calmariensis even
though this species had the smallest Ne . The
correction for demography in this species may have been insufficient as
nearly neutral nonsynonymous substitutions have fixed during the recent
period of low Ne and thereby incorrectly indicate
evidence for positive selection. Third, an additional issue when
examining selection patterns is slightly deleterious mutations. If such
mutations are segregating in the population, it becomes difficult to
detect positive selection and the degree of positive selection will be
underestimated. Removing deleterious mutations during selection analysis
is therefore widely accepted, by either setter a threshold for removing
minor alleles (Fay, Wyckoff, & Wu, 2001; Fay et al., 2002) or by
removing singleton polymorphism (Bierne & Eyre-Walker, 2004; Proschel,
Zhang, & Parsch, 2006). We used the latter strategy, which resulted in
a skew towards positively selected sites while keeping the relative
ranking of α among species. Thus, while absolute α-estimates in our
analysis are certainly underestimated, the relative ranking among
species is likely correct.
To summarize, our study species diverged recently (Hambäck et al., 2013)
but have during a limited time evolved key phenotypic differences in the
defence against parasitoid wasps. These phenotypic differences allowed
us to pinpoint important links between natural selection on immune genes
and host-parasitoid interactions where the species with the highest
immunocompetence had the highest number of positively selected immune
genes. When examining candidate immune-related genes, our results are
consistent with the variable selection pattern of immune genes in other
studies (Heger & Ponting, 2007; Waterhouse et al., 2007), but also show
that selection acts on multiple genes in the immune pathways and not
only on a few major genes (Kraaijeveld et al., 1998). While the
evolutionary process underlying the speciation and immune system
differences is unknown, comparisons with studies involving experimental
evolution in microbial systems provide potential hypotheses. These
studies show that natural enemies generally are not expected to cause
host populations to differentiate sympatrically, and that the role of
natural enemies instead seems to be that the arms-race between host and
their natural enemies increases the rate of divergence between
allopatric populations (Buckling & Hodgson, 2007; Buckling & Rainey,
2002). If these processes are similar for the evolution of
host-parasitoid interactions in this beetle system, it is likely that
the phenotypic differences in encapsulation capacity initially evolved
during periods when the beetles where geographically isolated and
thereby coincided with an allopatric speciation process. In either case,
it seems likely that the observed phenotypic and genotypic differences
and the different relationships among species with the parasitoid wasps
have been pivotal during the speciation process.