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.