Complex Genomic Landscape of Inversion Polymorphism in Europe’s Most Destructive Forest Pest Anastasiia Mykhailenko 1,2, Piotr Zieliński 1, Aleksandra Bednarz 1, Fredrik Schlyter 3,4, Martin N. Andersson 5, Bernardo Antunes 1, Zbigniew Borowski 6, Paal Krokene 7, Markus Melin 8, Julia Morales-García 1,2, Jörg Müller 9,10, Zuzanna Nowak 1, Martin Schebeck 11, Christian Stauffer 11, Heli Viiri 12, Julia Zaborowska 1, Wiesław Babik 1, Krystyna Nadachowska-Brzyska 1,* 1Institute of Environmental Sciences, Faculty of Biology, Jagiellonian University, 30-387 Kraków, Poland 2Doctoral School of Exact and Natural Sciences, Jagiellonian University, 30-348 Kraków, Poland 3Chemical Ecology, Department of Plant Protection Biology, Swedish University of Agricultural Sciences Alnarp, 234 22 Lomma, Sweden 4ETM, Faculty of Forestry and Wood Sciences, Czech University of Life Sciences Prague, 165 00 Praha, Czechia 5Department of Biology, Lund University, 223 62 Lund, Sweden 6Departament of Forest Ecology, Forest Research Institute, 05-090 Raszyn, Poland 7Division of Biotechnology and Plant Health, Norwegian Institute of Bioeconomy Research, 1433 Ås, Norway 8Forest Health and Bidiversity Group, Natural Resources Institute Finland, 80100 Joensuu, Finland 9Field Station Fabrikschleichach, Animal Ecology and Tropical Biology, Biocenter, University of Würzburg, 96181 Rauhenebrach, Germany 10Bavarian Forest National Park, 94481 Grafenau, Germany 11Institute of Forest Entomology, Forest Pathology and Forest Protection, Department of Forest and Soil Sciences, University of Natural Resources and Life Sciences Vienna (BOKU), 1190 Vienna, Austria 12UPM Forest, UPM-Kymmene, 33100 Tampere, Finland *Corresponding author: E-mail: krystyna.nadachowska@uj.edu.pl. Accepted: December 02, 2024 Abstract In many species, polymorphic genomic inversions underlie complex phenotypic polymorphisms and facilitate local adaptation in the face of gene flow. Multiple polymorphic inversions can co-occur in a genome, but the prevalence, evolutionary significance, and limits to complexity of genomic inversion landscapes remain poorly understood. Here, we examine genome-wide genetic variation in one of Europe’s most destructive forest pests, the spruce bark beetle Ips typographus, scan for polymorphic inver- sions, and test whether inversions are associated with key traits in this species. We analyzed 240 individuals from 18 populations across the species’ European range and, using a whole-genome resequencing approach, identified 27 polymorphic inversions covering ∼28% of the genome. The inversions vary in size and in levels of intra-inversion recombination, are highly polymorphic across the species range, and often overlap, forming a complex genomic architecture. We found no support for mechanisms such as directional selection, overdominance, and associative overdominance that are often invoked to explain the presence of large inversion polymorphisms in the genome. This suggests that inversions are either neutral or maintained by the combined action of multiple evolutionary forces. We also found that inversions are enriched in odorant receptor genes encoding elements of recognition pathways for host plants, mates, and symbiotic fungi. Our results indicate that the genome of this major forest pest of growing social, political, and economic importance harbors one of the most complex inversion landscapes described to date and raise questions about the limits of intraspecific genomic architecture complexity. Key words: polymorphic inversions, genome complexity, spruce bark beetle, Ips typographus, forest pest. © The Author(s) 2024. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 1 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Introduction Large structural variants (SVs), such as chromosomal inver- sions, translocations, insertions, and duplications, were among the earliest mutations to be described in natural po- pulations (McClung 1905; Sturtevant 1921; Dobzhansky and Sturtevant 1938; Dobzhansky 1970), but their system- atic identification and in-depth analysis have only become possible with recent advances in sequencing technology (Wellenreuther et al. 2019). The increasing accumulation of genomic data has revealed not only the presence of structural variation between and within many species but also its important role in species adaptation and speciation (Lamichhaney et al. 2015; Van’t Hof et al. 2016; Cheng et al. 2018; Wellenreuther and Bernatchez 2018; Faria, Chaube, et al. 2019; Todesco et al. 2020; Mérot et al. 2021; Harringmeyer and Hoekstra 2022; Saitou et al. 2022). One particular form of SVs, polymorphic chromo- somal inversions, has been at the center of recent debate about their potential to influence the evolutionary process (Berdan et al. 2023). Polymorphic inversions are chromosomal segments that occur in two orientations within populations: collinear and inverted haplotypes/arrangements. Inversions have been shown to be involved in speciation, local adaptation, and maintenance of complex phenotypes (Lamichhaney et al. 2015; Lohse et al. 2015; Wellenreuther and Bernatchez 2018; Fuller et al. 2019). This is due to a key property of inversions: they suppress recombination within heterozy- gotes and thereby prevent the separation of alleles present in each inversion haplotype. Because of their role as recom- bination modifiers, inversions can act as supergenes—large elements of genomic architecture containing multiple linked functional elements (Thompson and Jiggins 2014). Supergenes keep alleles together in the face of gene flow and suppress the formation of recombinant genotypes. Classic examples of supergenes include inversions asso- ciated with different mating strategies in ruffs (Calidris pugnax; Lamichhaney et al. 2015), mimicry phenotypes in Heliconius butterflies (Joron et al. 2011), and social organ- ization in fire ants (Solenopsis spp.; Purcell et al. 2014; Gutiérrez-Valencia et al. 2021). In many other species, poly- morphic inversions define locally adapted ecotypes (Lowry and Willis 2010; Kirubakaran et al. 2016; Koch et al. 2021; Matschiner et al. 2022; Reeve et al. 2023) or exhibit spatial haplotype frequency differences, e.g. by forming geographic and climatic gradients in inversion distributions (Ayala et al. 2014, 2017; Kapun et al. 2016; Mérot et al. 2021). While most of the described cases are organisms with one or a few inversions, several recent studies have re- ported species with many polymorphic inversions (Littorina snails, Faria, Chaube, et al. 2019; Reeve et al. 2023; sun- flowers, Todesco et al. 2020; deer mice, Harringmeyer and Hoekstra 2022; humans, Porubsky et al. 2022). These recent findings raise questions about the prevalence and evolutionary significance of polymorphic inversions in nat- ural populations. Are inversion-rich genomes the exception or the rule? How much of the genome can be situated with- in polymorphic inversions and, consequently, how large can the fraction of the genome that undergoes reduced recom- bination be? The latter question is particularly important because, in addition to suppressing recombination and keeping coadapted alleles together, inversion heterozy- gotes also suppress the formation of new allelic combina- tions and thus reduce the efficacy of natural selection (Roesti et al. 2022). Long-term persistence of two inversion haplotypes can be facilitated by both divergent and balancing selection (Faria, Johannesson, et al. 2019). Divergent selection can favor different inversion genotypes in different environ- ments and, when coupled with reduced migration between divergent populations, can lead to speciation. Even when intraspecific gene flow is high, divergent selection can lead to divergent ecotypes associated with locally advanta- geous inversion genotypes. Alternatively, balancing selec- tion may maintain balanced inversion polymorphisms over time via several, not mutually exclusive, mechanisms, such as overdominance, negative frequency dependence, antag- onistic pleiotropy, and spatially or temporally varying selec- tion (Connallon and Clark 2014; Faria, Johannesson, et al. 2019). Importantly, regardless of the selection mechanism, Significance The spruce bark beetle (Ips typographus) is one of the most destructive forest pests, causing mass mortality of spruce stands under favorable weather conditions. Here, we performed a large-scale population genomic analysis and provided information on the genetic structuring and variation along the species genome. We found that the spruce bark beetle harbors one of the most complex polymorphic inversion landscapes described to date. We also found that inversions are enriched in odorant receptor genes, which are important for recognizing hosts, mates, and symbiotic fungi. Our study raises questions about the limits of the complexity of intraspecific genome rearrangements and the causes and conse- quences of inversion-rich genomes. Our findings not only shed light on the genomic variation of the major forest pest but also contribute to the understanding of the abundance of inversion polymorphisms that are increasingly found across the tree of life. Mykhailenko et al. GBE 2 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 inversions will accumulate mutations independently since recombination is suppressed in inversion heterozygotes. This will lead to differentiated allelic content and increased differentiation between inversion haplotypes over time (Faria, Johannesson, et al. 2019). Each inversion haplotype can thus be treated as a separate “population” with a size that corresponds to the frequency of that arrangement within the studied population or species. Rare inversion ar- rangements will experience a high mutational load due to reduced recombination and limited purging because there are few homozygotes. However, deleterious mutations can also accumulate on more frequent inversion haplotypes (Berdan et al. 2021) and lead to associative overdominance. This helps to maintain inversion polymorphisms since independent accumulation of mutations continues over time, but recessive deleterious alleles private to an inver- sion arrangement are invisible to selection in inversion heterozygotes. The Eurasian spruce bark beetle (Ips typographus [L.]: Curculionidae: Scolytinae; hereafter the spruce bark bee- tle), plays a key role in Eurasian forest ecosystems. Under endemic conditions, it attacks weakened Norway spruce (Picea abies [L.] H. Karst.) trees. However, if spruce resist- ance is compromised by certain abiotic disturbances (e.g. snowbreaks, windfalls, high temperatures, and drought), an increased availability of stressed trees can trigger mass propagation of beetles and lead to rapid population in- crease and devastating outbreaks. The extent of recent beetle outbreaks in Europe is unprecedented, and impacts will likely increase in the coming decades in response to cli- mate change (Biedermann et al. 2019; Bentz et al. 2021; Müller et al. 2022). For example, during the first decade of this century, the spruce bark beetle killed an estimated 14.5 million m3 of timber per year on average. The Czech Republic provides a particularly striking example of the bee- tles’ destructive potential. During the peak outbreak years of 2017 to 2019, the beetles killed 3.1% to 5.4% of the country’s entire growing stock of Norway spruce each year, translating to 23 million m3 spruce killed in 2019. Historically, Central Europe has been most heavily affected by bark beetle outbreaks, while outbreaks have been less frequent and destructive in northern Europe (Hlásny et al. 2019). However, this may change as climate warming is ex- pected to make the boreal forests of Northern Europe more vulnerable to bark beetle outbreaks (Lange et al. 2006; Jönsson et al. 2007; Bentz et al. 2021; Müller et al. 2022). As an example, heatwaves and severe summer drought in Sweden in 2018 initiated a bark beetle outbreak killing over 30 million m3 Norway spruce in the next years (Öhrn et al. 2021). Increasing bark beetle attacks and other forest distur- bances have already triggered social and political conflicts in parts of Europe and have highlighted the urgent need for improved management strategies (Vega and Hofstetter 2014; Hlásny et al. 2021). The increasing import- ance of the spruce bark beetle is reflected in a rapidly grow- ing body of research focuses on the species’ ecology and the causes and consequences of outbreaks (reviewed in Vega and Hofstetter 2014; Biedermann et al. 2019; Hlásny et al. 2021). Despite this enormous interest, one as- pect of the species’ biology remains largely unexplored: we know almost nothing about the species’ genome-wide vari- ation and the evolutionary mechanisms that shape this vari- ation. The lack of population genomics studies that analyze whole-genome variation restricts our understanding of the genomic basis of adaptation and adaptive potential in the spruce bark beetle. This is particularly important because, as revealed in this study, the spruce bark beetle genome har- bors a complex inversion polymorphism landscape that may play a critical role in many evolutionary processes, including adaptation (discussed above). Here, we investigated genome-wide variation across spruce bark beetle populations with a special focus on chromosomal inversion polymorphisms. We found one of the most complex polymorphic inversion architectures de- scribed to date and investigated several evolutionary me- chanisms, including directional selection, overdominance, and associative overdominance, that can maintain inversion polymorphisms in the genome. We also tested associations between inversion polymorphisms and two key fitness- related traits in the spruce bark beetle—diapause and olfac- tion. Our results suggest that inversions may play a role in the evolution of key traits in this species and raise questions about the prevalence and role of inversion polymorphisms in species characterized by large effective population sizes, little geographic subdivision, and no clear phenotypic vari- ation across their geographical range. Results and Discussion We resequenced whole genomes of the spruce bark beetle from 18 European populations (Fig. 1; supplementary table S1, Supplementary Material online). After initial raw data processing from 253 individuals (including data quality check, trimming and Genome Analysis Toolkit (GATK)-guided variant calling, genotyping, recalibration, and filtering; for details, see Supplementary Material), we used 240 individuals (141 females and 99 males) in down- stream analyses. The mean per individual sequencing depth was 23.2× (range: 5 to 53×). Sequencing coverage on IpsContig9 was consistently lower in individuals sexed as males (on average 0.57 individual coverage). Thus, we con- sidered IpsContig 9 to be a sex (X) chromosome (females carry XX and males XYp chromosomes). After quality filter- ing, we retained 5.245 million single nucleotide poly- morphisms (SNPs) covering the entire length of the genome assembly (236.8 Mb) but analyzed a subset of 5.067 million SNPs located on the 36 longest contigs and Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 3 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 representing 78% of the assembly. Alignment of the I. typographus and Ips nitidus genome assemblies indicated no major misassembly problems (I. typographus contigs align to single I. nitidus chromosomes) and demonstrated high synteny between karyotypes of both species (supplementary fig. S1, Supplementary Material online). Given these results on contig sizes, genome assembly qual- ity scores, and species karyotypes, we find it likely that many large spruce bark beetle contigs represent entire chromosome arms. Such a high level of genome quality is sufficient for all downstream analyses, in particular SNP-based polymorphic inversion identification (Mérot et al. 2021; Reeve et al. 2023). Complex Genomic Inversion Landscape Twenty-nine candidate inversions (“inversions” hence- forth) were identified following the criteria described in the Materials and Methods section (Table 1; Fig. 2; supplementary figs. S2 and S3, Supplementary Material on- line). Briefly, we considered a genomic region to harbor a polymorphic inversion if local principal component analysis (PCA) identified the region as an outlier and/or the region exhibited high linkage disequilibrium (LD) and PCA per- formed on SNPs from this region separated individuals not by geography but into three distinct groups (with one group classified as putative heterozygous individuals char- acterized by higher heterozygosity compared to the other two groups). Two putative inversions (Inv16.1 and Inv23.1) were most likely part of the same inversion: they are in strong LD (supplementary fig. S4, Supplementary Material online) and the same beetles were always geno- typed as homo- and heterozygotes in both, which is unlike- ly to occur for two unlinked inversions. The same situation was found for Inv16.2 and Inv23.2. No other inversions were in strong LD with each other, which would suggest cosegregation (supplementary fig. S4, Supplementary Material online). Thus, overall, we found 27 inversions in 17 contigs, including one located on the X chromosome (Inv9). Approximate inversion sizes varied from 0.1 to 10.8 Mb (Table 1) and inversions constituted ∼28% of the analyzed part of the genome (48 Mb of 170 Mb) and at least 18.6% of the entire assembly. Estimated inversion ages ranged from 0.5 to 2.6 My (Table 1, assuming a mu- tation rate of 2.9 × 10−9; see supplementary table S1, Supplementary Material online for results for different mu- tation rates), and younger inversions had higher major in- version haplotype frequency (supplementary fig. S5, Supplementary Material online). Inversion regions exhibited a reduced population recombination rate in heterozygous individuals (supplementary fig. S2, Supplementary Material online) and moderate to high genetic differenti- ation between inversion arrangements (FST between homo- zygous individuals was 0.15 to 0.64; Fig. 1; supplementary fig. S2, Supplementary Material online). While the variation patterns observed on contigs with single inversions (12 contigs) were clear and left no doubt about the presence of polymorphic inversions, we also Fig. 1. Genomic structure and differentiation in I. typographus. a) Whole-genome PCA, where colors correspond to what population an individual beetle belongs to. Colors indicate geographical grouping of populations: S, south; P, Poland; N, north. b, c) Geographical distribution and genetic differentiation of the 18 beetle populations analyzed (see supplementary table S11, Supplementary Material online for details about the sampling locations). The two main genetic groups are shown in red (southern) and blue (northern). d) Genome-wide genetic differentiation (FST) calculated between northern and southern groups. Vertical lines separate different contigs (shown in gray and black). All panels were prepared using data that included inversions. Mykhailenko et al. GBE 4 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 observed more complex patterns such as overlapping inver- sions and inversions with double-crossover events. These patterns were more difficult to interpret and deserve add- itional explanation and caution (Fig. 2; supplementary figs. S2 and S3, Supplementary Material online). Putative overlapping inversions were identified on the basis of overlapping clusters of high LD, PCA clustering into three distinct groups along PC2 or PC3 (in addition to the three distinct groups along PC1), and moderate to high FST between individuals classified as alternative homo- zygotes (Fig. 2; supplementary fig. S3, Supplementary Material online). In summary, IpsContig14 and IpsContig22 contained complexes of multiple adjacent and sometimes overlapping inversions (Fig. 2j to l; supplementary figs. S2 and S3, Supplementary Material online), and three other contigs contained two overlapping inversions each (IpsContigs: 7, 16, and 23; supplementary figs. S2 and S3, Supplementary Material online). In the case of putative double-crossover events, the observed LD clusters were separated by regions of lower LD and intermediate groups of individuals were visible between the main clusters along PC1 (supplementary fig. S3, Supplementary Material online). The FST between indivi- duals classified as alternative homozygotes was significant- ly reduced in the region of the putative double crossover (Fig. 2; supplementary figs. S2 and S3, Supplementary Material online). In total, four putative double-crossover events were identified within four inversions (Inv5, Inv18, Inv22.3, and Inv22.4; Fig. 1; supplementary figs. S2 and S3, Supplementary Material online). The most difficult patterns to disentangle were asso- ciated with Inv2, Inv5, Inv7.1, and Inv7.2. Putative Inv2 showed a signature of overlapping LD clusters (overlapping inversions) but ambiguous PCA clustering (supplementary fig. S3, Supplementary Material online). The patterns were clear and consistent with inversion polymorphisms when only populations from the southern region were ana- lyzed (supplementary fig. S6, Supplementary Material online; thus, Inv2 was genotyped only in this part of the species range). This suggests the existence of other Table 1 Identified chromosomal inversions in the I. typographus genome ID Contig Size Start End Age Odorant receptors Inv2 IpsContig2 4.04 12.67 16.71 0.5 … Inv3 IpsContig3 0.14 1.11 1.25 1.1 … Inv5 IpsContig5 10.84 0.00 10.84 1.3 ItypOR33, ItypOR41, ItypOR40, ItypOR10, ItypOR47, ItypOR50, ItypOR29, ItypOR43JOI, ItypOR34, ItypOR52NTE, ItypOR4, ItypOR3, ItypOR53, ItypOR2, ItypOR19 Inv6 IpsContig6 0.34 8.93 9.27 1.0 … Inv7.1 IpsContig7 0.67 0.00 0.67 1.7 … Inv7.2 IpsContig7 6.92 0.00 6.92 1.1 ItypOR1, ItypOR17 Inv9 IpsContig9 3.30 1.71 5.01 0.6 … Inv10 IpsContig10 0.08 6.05 6.13 1.8 … Inv12 IpsContig12 0.07 3.63 3.70 1.5 … Inv13 IpsContig13 4.50 0.00 4.50 1.7 ItypOR28, ItypOR23, ItypOR49, ItypOR27 Inv14.1 IpsContig14 2.08 0.00 2.08 1.9 ItypOR36, ItypOR44, ItypOR18JF, ItypOR20NTE Inv14.2 IpsContig14 0.67 2.08 2.75 1.0 … Inv14.3 IpsContig14 0.76 2.78 3.54 2.1 … Inv14.4 IpsContig14 0.11 3.73 3.84 2.1 … Inv14.5 IpsContig14 0.57 4.23 4.80 2.3 … Inv14.6 IpsContig14 2.48 0.00 2.48 0.6 ItypOR36, ItypOR44, ItypOR18JF, ItypOR20NTE Inv15 IpsContig15 1.92 0.86 2.78 1.9 … Inv16.1 IpsContig16 4.83 0.00 4.83 1.6 ItypOR58, ItypOR9, ItypOR11, ItypOR31, ItypOR30, ItypOR16, ItypOR35JF Inv16.2 IpsContig16 4.83 0.00 4.83 0.7 ItypOR58, ItypOR9, ItypOR11, ItypOR31, ItypOR30, ItypOR16, ItypOR35JF Inv17 IpsContig17 0.21 2.48 2.69 2.2 … Inv18 IpsContig18 2.32 0.00 2.32 2.1 … Inv22.1 IpsContig22 0.32 0.00 0.32 2.1 … Inv22.2 IpsContig22 0.23 0.42 0.65 2.1 … Inv22.3 IpsContig22 1.23 0.68 1.91 2.0 ItypOR22CTE Inv22.4 IpsContig22 0.20 1.92 2.12 2.1 … Inv22.5 IpsContig22 2.24 0.00 2.24 0.6 ItypOR22CTE Inv23.1 IpsContig23 2.10 0.00 2.10 1.4 ItypOR58, ItypOR9 Inv23.2 IpsContig23 1.84 0.26 2.10 0.6 ItypOR58, ItypOR9 Inv26 IpsContig26 0.10 0.12 0.22 2.6 … Note that Inv16.1 and Inv23.1 are parts of the same inversion and Inv16.2 and Inv23.2 are a part of another single inversion. ID, inversion name; Contig, contig name; Size, size of the inversions (Mb); Start and End, coordinates of the inversion (Mb); Age, approximate age of the inversion in Myr; Odorant receptors, odorant receptors present within inversion (in the order they appear along the contig sequence). Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 5 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Fig. 2. Identification of chromosomal inversions in I. typographus. Each row of figure panels shows the results of per-contig PCA a, d, g, j), per-contig LD analysis b, r, h, k), and genetic differentiation (FST) analysis for a selected contig c, f, i, l). Dots in the PCA panels represent individual beetles and are colored according to the heterozygosity of the individual (darker colors represent lower heterozygosity). a to c) Results for a contig with no inversions (IpsContig8) and little differentiation between southern and northern populations (c, FST between southern and northern populations). d to l) Different contigs with increasingly complex inversion patterns: a single inversion (d to f, Inv15); a single inversion with a double-crossover signal (g to i, Inv18); and multiple adjacent inversions with one inversion overlapping with the first two inversions on the contig (j to l; Inv14.1, Inv14.2, Inv14.3, Inv14.4, Inv14.5, and Inv.14.6). m) Contigs with inversions indicated in red (contigs for which inversion genotyping was not possible are not shown; see Materials and Methods section for details). Colored dots indicate the inversions that are presented in detail in a) to l). Colored bars e, h, and k) indicate the position of each inversion. The same colors are used in f), i), and l), where FST is calculated between major (MJ) and minor (MN) inversion haplotypes. Two additional FST lines (blue and green) in i) show FST calculated between inversion homozygotes and a recombinant haplotype (R) formed after a putative double-crossover event. Both axes in b), e), h), and k) represent positions along the contig in megabases; low levels of linkage are shown in blue and higher levels in yellow to dark red. Sex chromosome is indicated with a star in m). Mykhailenko et al. GBE 6 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 structural differences between the southern and northern regions or unidentified genome assembly problems. PCA of Inv5 identified three clusters that most likely correspond to three inversion genotypes, but it also identified several smaller clusters. The smaller clusters may indicate multiple independent recombination events (double-crossover events that occurred at different locations within Inv5; supplementary figs. S2 and S3, Supplementary Material online) or assembly problems within Inv5 (unlikely, since IpsContig5 aligns perfectly with I. nitidus chromosome 9; supplementary fig. S1, Supplementary Material online). Thus, Inv5 was only genotyped for three major genotype groups for which multiple lines of evidence (LD, heterozygos- ity, and FST) suggested an inversion polymorphism. Inv7.1 and Inv7.2 are characterized by overlapping LD clusters, and PCA and FST patterns consistent with overlapping inver- sions. However, Inv7.2 showed low LD in the middle of the inversion and unexpected clustering of individuals in one of the clusters along PC1. Here, the cluster was split into two distinct groups where individuals clustered according to their geographic location within each group (supplementary figs. S2 and S3, Supplementary Material online). This suggests possible misassembly, and Inv7.2 is probably shorter than predicted from the size of the LD regions. Genotyping of Inv16.2 and Inv23.2 was only possible for part of the species range, as PCA clustering into three genotype groups was only visible when analyzing the northern group (supplementary fig. S6, Supplementary Material online). Even though we cannot rule out that some of the inver- sion patterns we identified were formed or distorted by misassembly or additional structural rearrangements (e.g. duplications; Kim et al. 2022), there is strong evidence for most of the polymorphic inversions in the spruce bark beetle. Extensive collinearity with I. nitidus suggested no ma- jor misassembly problems and supports our identification of polymorphic inversions in the spruce bark beetle genome (supplementary fig. S1, Supplementary Material online; Wang et al. 2023). Nevertheless, long-read sequencing and Hi-C data are needed to shed light on the few problem- atic/complex cases and to precisely identify the inversion boundaries. It was not surprising that we detected poly- morphic inversions in the spruce bark beetle genome, as there are many well-known examples of polymorphic inver- sions in natural populations (Wellenreuther and Bernatchez 2018). What was striking, however, was the complexity of the genomic landscape of polymorphic inversions we found in this species. The spruce bark beetle has at least 27 large inversions covering a substantial part of the genome (28% of the analyzed part of the genome). Numerous (a dozen or more) polymorphic inversions have been described for sev- eral species (e.g. Faria, Chaube, et al. 2019; Tigano et al. 2021; Harringmeyer and Hoekstra 2022), but it remains an open question whether many polymorphic inversions within species are the exception or the rule. The exceptionally complex inversion architecture we found in the spruce bark beetle, with multiple adjacent and often overlapping inversions, resembles well-known examples from Heliconius butterflies (Jay et al. 2021) and fire ants (Wang et al. 2013). In these insects, multiple adja- cent inversions are the basis for mimicry phenotypes and complex social organization, respectively. The presence of clusters of adjacent inversions and inversion overlaps are consistent with theoretical expectations of stepwise exten- sion of recombination suppression on supergenes (Jay et al. 2022) and with a highly polygenic architecture of adapta- tion (Schaal et al. 2022). Genome-Wide Variation versus Inversion Region Variation and Its Geographic Structure Analyses based on the whole genome (including inversions and collinear regions) revealed a clear latitudinal structuring of genetic variation in the spruce bark beetle (Fig. 1). NGSadmix supported the presence of two distinct genetic groups corresponding to southern and northern popula- tions, with Polish populations showing varying degrees of admixture between the two groups (supplementary fig. S7, Supplementary Material online). Based on these results, and to be able to assess the differentiation between two genetic clusters, we divided the studied populations into a northern and southern group, excluding admixed Polish populations. Despite unambiguous NGSadmix division into two genetic clusters, the genome-wide genetic differentiation between the northern and southern groups was very low (FST = 0.021). This was true also within inversion regions where FST ranged from 0.005 to 0.05. Similarly, pairwise FST between populations showed low levels of differentiation, ranging from 0.000 to 0.035 (calculated excluding inversions; supplementary table S2, Supplementary Material online). Mean genome-wide nu- cleotide diversity was moderate (π = 0.0062) and per popu- lation π ranged from 0.0055 to 0.0066 (supplementary fig. S8, Supplementary Material online). There was a weak negative correlation between nucleotide diversity and latitude (r2 = 0.32, P = 0.033; supplementary fig. S9, Supplementary Material online), as northern populations had slightly lower genetic variation than southern popu- lations (πsouthern = 0.0065; πnorthern = 0.0061; πPolish = 0.0066). Southern populations had an excess of rare alleles and, consequently, had more negative Tajima’s D-values along the genome than northern populations (mean Tajima’s D was −0.458 and −0.062 in the southern and northern groups, respectively; supplementary fig. S8, Supplementary Material online). All these results are con- sistent with previous phylogeographic studies of the spruce bark beetle that analyzed a much smaller number of genet- ic markers. The data suggests high levels of connectivity among populations and a very recent differentiation into Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 7 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 two genetic clusters (Stauffer et al. 1999; Sallé et al. 2007; Bertheau et al. 2013; Mayer et al. 2015). More recent RADseq data confirm a very weak genetic structure in the spruce bark beetle across much of Sweden (Ellerstrand et al. 2022), as is expected in a species with high dispersal (Nilssen 1984) and/or recent divergence between popula- tions. Tajima’s D-values indicate a different demographic history of the southern and northern groups (e.g. recent demographic expansion of southern populations and more stable demography in northern populations; supplementary fig. S8, Supplementary Material online). Inversion regions in the spruce bark beetle did not struc- ture in the same way as the collinear part of the genome and most regions do not show any clustering into a southern and northern group (supplementary fig. S2, Supplementary Material online). Frequencies of two inver- sions were significantly correlated with latitude (Fig. 3; see further discussed below). Almost all identified inver- sions were polymorphic across the beetle’s European range, except for one inversion (Inv9) that was polymorphic only in northern populations. For three inversions (Inv2, Inv5, and Inv16.2 + Inv23.2), unambiguous genotyping was only possible across part of the species range (see above). The differentiation between haplotypes of some in- versions was moderate or high (mean FST range 0.15 to 0.64; Fig. 2; supplementary fig. S2, Supplementary Material online), and, according to our age estimates, the inversions originated before the Last Glacial Period. Many inversions may be several million years old and probably predate the within-species differentiation into a southern and northern group. This was not unexpected, as many known inversion polymorphisms have been segregating within species for hundreds of thousands or millions of years (see table in Wellenreuther and Bernatchez 2018), sometimes even persisting through multiple speciation events (Brelsford et al. 2020). Association between Inversions and Fitness-Related Traits Inversion polymorphism is often associated with the main- tenance of complex polymorphic phenotypes (Schwander et al. 2014; Lamichhaney et al. 2015). Although the spruce bark beetle does not exhibit easily identifiable phenotypes, such as distinct color patterns or mating strategies, we were able to test for associations between inversions and two complex traits of key importance for many insects: diapause and olfaction. Diapause is an effective strategy to increase fitness and avoid mortality by synchronizing seasonal occurrence with critical resources and mitigating harmful effects of adverse abiotic conditions (Denlinger 2022). Diapause is character- ized by suppressed development, metabolism, and repro- duction, accompanied by increased stress resistance (Denlinger 2022; Schebeck et al. 2024). While multiple en- vironmental factors can influence this complex process, many aspects of diapause, such as its induction and termin- ation (i.e. the decision to enter into or exit from diapause, respectively), are heritable (Roff 1996). These traits may be controlled by a small number of loci or be highly poly- genic (Paolucci et al. 2016; Koštál et al. 2017; Pruisscher et al. 2018; Kozak et al. 2019; Dowle et al. 2020; Denlinger 2022). The spruce bark beetle enters reproduct- ive diapause as an adult to respond to unfavorable winter conditions (Schebeck et al. 2017, 2024). This diapause is physiologically and behaviorally characterized by sup- pressed development and reproduction, a lack of flight ac- tivity, an increased resistance against cold temperatures, and movement to suitable overwintering habitats (Schopf 1985; Annila 1969; Schopf 1989; Doležal and Sehnal 2007; Dworschak et al. 2014; Schebeck et al. 2017). The spruce bark beetle exhibits two diapause phenotypes that could be associated with polymorphic inversions: a faculta- tive photoperiod-regulated diapause and an obligate photoperiod-independent diapause (Schebeck et al. 2022). Beetles enter facultative diapause when day-length drops below specific threshold values that vary along latitu- dinal gradients. When day length is above the threshold, beetles develop and reproduce normally (Doležal and Sehnal 2007; Schroeder and Dalin 2017; Schebeck et al. 2022). Obligate diapausing beetles, however, enter dia- pause in each generation as a regular part of their life cycle, independently of day-length cues. These beetles can there- fore produce only one generation per year and do not re- sume development and reproduction until the next year (Schebeck et al. 2022). Both diapause phenotypes are usually found in European populations, with facultative diapause prevailing in more southern latitudes and obligate diapause being most common in northern regions (Schebeck et al. 2022). We found no association between diapause phenotypes and inversion genotypes (exact G-test; supplementary table S3 and fig. S10, Supplementary Material online), nor did we find any genom- ic regions that were highly differentiated between faculta- tively and obligately diapausing individuals (supplementary fig. S11, Supplementary Material online). However, it should be noted that due to a small sample size, we were only able to test for complete association between phenotype and in- version genotype. Olfaction is another fitness-related trait in many insects, including bark beetles (Raffa et al. 2016). Insect odorant receptors (ORs) are encoded by a large and dynamically evolving gene family. Some ORs are evolutionarily con- served across species within insect orders, whereas many others are species or genus specific. For the spruce bark beetle, odorant detection is essential for host and mate finding, as well as for recognition and maintenance of symbiotic fungal associations (Andersson et al. 2009; Mykhailenko et al. GBE 8 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Kandasamy et al. 2019). We examined 73 antennally ex- pressed ORs (Yuvaraj et al. 2021) and found that 46% of these (33 of 71 ORs, which mapped to 26 IpsContigs) were located within inverted regions (Table 1). A permuta- tion test (10,000 iterations; randomizing inversion locations across 170 Mb) validated the overrepresentation of OR genes within inversions (P = 0.0325; i.e. in just 325 of the 10,000 permutations, the number of ORs located within in- versions was higher than or equal to the number observed in the real data). In addition, several Ips-specific ORs (five out of seven ORs from an Ips-specific OR clade; Hou et al. 2021) were located within inverted regions; specifically, four out of these five were predominantly on IpsContig13. The five Ips-specific ORs (ItypOR23, ItypOR27, ItypOR28, ItypOR29, and ItypOR49) have been functionally characterized and respond to compounds pro- duced primarily by beetles (pheromones), the host tree, or fungal symbionts (Hou et al. 2021; Powell et al. 2021). We found no difference in OR number between inversion haplotypes (i.e. there were no OR deletions). Most ORs located in inversions harbor multiple nonsynonymous variants segregating within the studied spruce bark beetle populations (supplementary table S4, Supplementary Material online) and nine ORs (supplementary table S5, Supplementary Material online) had at least one fixed or nearly fixed (dxy > 0.9) nonsynonymous variant between in- version haplotypes. These diverged ORs were located at Inv5 (ItypOR43JOI and ItypOR29), Inv13 (ItypOR23 and ItypOR27), Inv14.1 (ItypOR20NTE and ItypOR36), and Inv16.1 (ItypOR30, ItypOR31, and ItypOR16). Among these ORs only ItypOR30 had dN/dS > 1 (i.e. signature of positive selection; two nonsynonymous substitutions, McDonald– Kreitman test, P = 0.55). Olfactory receptor activity was one of the enriched gene ontology (GO) categories found in inversions showing sig- nificant latitudinal variation (Fig. 3; supplementary table S6, Supplementary Material online; 0.01 < P < 0.05) and one of these inversions harbored multiple OR genes (Fig. 3; Inv13, Table 1). In Drosophila pseudoobscura, ORs are also associated with inversions (Fuller et al. 2016, Fig. 3. Correlation between inversion haplotype frequency and latitude for chromosomal inversions detected in European I. typographus populations sampled along a latitudinal gradient. Significant correlations are indicated with an asterisks (*P < 0.05, after Hommel multiple testing correction). Inversions harboring ORs are indicated in green. Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 9 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 2017). Interestingly, Inv13 includes a gene-encoding ItypOR23, a receptor that has been previously shown to pri- marily respond to fungal volatiles (Hou et al. 2021). This OR was also one of the most divergent ORs between alternative inversion haplotypes (supplementary table S5, Supplementary Material online). We hypothesize that dif- ferent OR alleles associated with Inv13 may be involved in spruce bark beetle interactions with fungal associates pre- sent across Europe. As many other bark beetle species, the spruce bark beetle is associated with several spatially differentiated fungal symbionts that can help beetles to ex- haust and overcome tree defenses (Zhao et al. 2019). Interestingly, recent studies suggest that beetles can recog- nize volatile compounds emitted by certain beneficial fungi and use these volatiles to locate and pick up the fungi. Furthermore, individual beetles may also have preferences for different fungal species (Kandasamy et al. 2019; Zhao et al. 2019; Kandasamy et al. 2023). Spruce bark beetles in Germany are, for example, more attracted to fungal spe- cies that are common in Germany (such as Grosmannia penicillata and Endoconidiophora polonica) than to rarer species (Leptographium europhioides). Unpublished data also suggest that Swedish beetles are more attracted to L. europhioides, which is common in Sweden (personal communication, D. Kandasamy). Although preliminary, these observations suggest that beetle preferences may be tuned to the local fungal community, and we speculate that inversions may be involved in the recognition of region-specific fungal species. Diapause and olfaction are not the only polymorphic complex traits in I. typographus. Several other potentially interesting traits include the existence of pioneer individuals that are the first to infest host trees and re-emergence of some females after egg laying to establish so-called sister broods in new trees (Anderbrant and Lofqvist 1988; Anderbrant 1989; Lehmanski et al. 2023). Other less obvi- ous/visible phenotypes could be tested for their association with inversion polymorphisms. GO enrichment analysis pro- vided a list of gene categories of interest (supplementary table S7, Supplementary Material online), but comprehen- sive research is needed to understand the relationship be- tween spruce bark beetle phenotypes and the species’ complex inversion polymorphism landscape. Evolutionary Mechanisms Maintaining Inversion Polymorphism in the Spruce Bark Beetle Several nonmutually exclusive evolutionary mechanisms can maintain polymorphic inversions within species, the most important mechanisms being divergent and balancing selection (Wellenreuther and Bernatchez 2018; Faria, Johannesson, et al. 2019). The importance of divergent se- lection has been postulated based on allele frequency pat- terns and associations of polymorphic inversions with local adaptations that persist despite extensive intraspecific gene flow (Tigano and Friesen 2016). For example, a recent study of deer mice (Peromyscus maniculatus; Harringmeyer and Hoekstra 2022) identified multiple polymorphic inversions with clinal variation across environmental gradients in two distinct habitats. Such frequency changes have also been reported across hybrid zones (Faria, Chaube, et al. 2019) or latitudinal gradients (Mérot et al. 2021) in other species. Although the spruce bark beetle does not occupy distinct environmental niches, it inhabits forests across a very wide latitudinal gradient (spanning at least 16°). It is also a species with high dispersal capacity and extensive gene flow, as indicated by low FST across its range. We found a significant correlation between the frequency of in- version haplotypes in populations and geographic location (latitude) for two inversions (Fig. 3; Inv12 and Inv13; r2 0.68 and 0.62, respectively). Importantly, in both cases, r2 was greater than what we observed for random SNPs with simi- lar frequencies but situated outside inversions (P = 0.003 and 0.021 for Inv12 and Inv13, respectively). For multiple SNPs (59), climate and land cover variation (summarized using PCA) were associated with allele frequency changes; however, no such association was found for inversion gen- otypes (Fig. 4). The identified SNPs were within 14 genes, five of which had assigned GO categories (supplementary table S8, Supplementary Material online). All the above re- sults indicate that there may be selection across environ- mental gradients in the spruce bark beetle but that this is not the only, or even a major, force maintaining inversion polymorphism within the species. It also remains an open question whether some inversions are involved in local adaptation or not. Although “local adaptation” is a major hypothesis pro- posed to explain inversion polymorphism (Christmas et al. 2019; Akopyan et al. 2022; Harringmeyer and Hoekstra 2022), balancing selection and related mechanisms may also be important in maintaining polymorphic arrange- ments (Faria, Johannesson, et al. 2019; Mérot et al. 2020; Berdan et al. 2021). Such mechanisms include overdomi- nance, associative overdominance, frequency-dependent selection, and selection that varies spatially and/or tempor- ally. We found no support for overdominance playing a role in the spruce bark beetle, as we detected no excess of inver- sion heterozygotes compared to Hardy–Weinberg (HW) ex- pectations in any of the populations. The only significant deviations from HW expectations were found in Inv12 (a deficit of heterozygotes across the whole species range) and Inv22.4 (more heterozygotes than expected across the whole species range and southern, Polish, and northern populations; supplementary table S9, Supplementary Material online). We therefore looked more closely at mu- tation load, which, if high, can indicate the importance of associative overdominance in maintaining inversion polymorphism. Theory predicts that recessive deleterious Mykhailenko et al. GBE 10 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 mutations will accumulate on both inversion arrangements but that most of these mutations will be private to only one arrangement (Navarro et al. 2000; Faria, Johannesson, et al. 2019; Berdan et al. 2021). This would lead to associative overdominance, as in heterozygotes, the effects of deleteri- ous recessive alleles on one arrangement would be masked by the wild-type alleles on the other arrangement. The re- sult would be long-term maintenance of the inversion poly- morphism, resulting in strong divergence between inversion haplotypes (Navarro et al. 2000; Guerrero et al. 2012; Berdan et al. 2021). Interestingly, several stable evolutionary scenarios that maintain polymorphic inversions are possible (for details, see figure 4 in Berdan et al. 2021). These scenarios differ in the expected mutation load, fitness, and frequency of the corresponding genotypes. Given the haplotype fre- quencies we observed in spruce bark beetle inversions, two scenarios are likely. First, that minority arrangements experience higher mutation load (due to reduced recom- bination and lower population size) but are maintained in the population at low frequencies due to, e.g. associative overdominance. Such a mechanism would favor balanced inversion polymorphisms of intermediate to large sizes (har- boring hundreds of genes; Ohta 1971; Connallon and Olito 2022) and has been shown to play a role in maintaining polymorphic inversions in several insect species (Yang et al. 2002; Jay et al. 2021). Second, mutation load may ac- cumulate on one or both inversion arrangements but be mi- tigated by the haplotype structuring process, i.e. the existence of multiple diverged subhaplotypes among inver- sion homozygotes that reduces the mutation load within homozygotes. If this process operates within one or both inversion arrangements, it may result in more equal fre- quencies of alternative inversion haplotypes. Such a mech- anism can operate when genetic variation and mutation load are high and are theoretically possible in the spruce bark beetle system (Berdan et al. 2021). Contrary to these theoretical expectations, we observed no sign of increased mutation load (measured by the πN/πS ratio) in inversion regions compared to the collinear part of the spruce bark beetle genome (supplementary table S10, Supplementary Material online; Fig. 5). We also found no sign of haplotype structuring that could reduce the mutation load within inversion homozygotes Fig. 4. GEAs across 36 contigs in the I. typographus genome. The y axis shows log-transformed P-values from LFMM analysis testing association between SNPs and three principal components (PC1 to 3, facet labels) summarizing the environmental variation between beetle populations. Each dot represents a single SNP and different colors different contigs. The upper three panels show results for all SNPs and the lower three panels show results using data where each inversion (blue circles) was coded as a single locus. Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 11 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 (supplementary fig. S12, Supplementary Material online). The only significant within-homozygote clustering we ob- served was in a few inversion haplotypes. However, this clus- tering divided individuals into southern and northern clades. This suggested either neutral differentiation or divergent se- lection acting on one of the inversion arrangements, rather than haplotype structuring being associated with mutation load (supplementary fig. S12, Supplementary Material online). These results are consistent with observations in other species of no significant mutation load when inversion haplotypes are subjected to divergent selection resulting in geographic structuring (Harringmeyer and Hoekstra 2022; Huang et al. 2022). In such cases, inversions facilitate adap- tive divergence but do not tend to accumulate a mutation load. However, the geographic structuring of inversion poly- morphisms is weak in the spruce bark beetle, and many in- versions appear to have a slightly lower mutation load associated with the more common inversion haplotype (Fig. 5; supplementary table S11, Supplementary Material online). It is possible that the accumulation of a mutation load in the spruce bark beetle is mitigated by a high effective population size in this species, as many homozygous indivi- duals in populations may purge mutation load to some ex- tent. It is also possible that the πN/πS ratio does not fully capture the patterns of mutation load in the species and that further and more in-depth analyses are needed. Overall, the absence of heterozygote excess does not support a role of overdominance in the maintenance of Fig. 5. Mutation load analysis in I. typographus. a) The ratio of nonsynonymous to synonymous nucleotide diversity (πN/πS) computed for 200-kb windows along collinear parts of the genome (COL) and different inversion haplotypes (blue: major haplotype; yellow: minor haplotype; green: haplotype produced by double-crossover events between minor and major haplotypes). n = 4 to 55, 200-kb windows per inversion haplotype and n = 531 for COL. For clarity, πN/πS outliers above 0.5 are not shown (a total of 58 values, 45 of them in COL; supplementary fig. S16, Supplementary Material online shows the plot with these outliers). The lower and upper box hinges correspond to the first and third quartiles, whiskers show 1.5× the interquartile range. b) Distribution of πN/πS values per 200-kb window in inversions and collinear part of the genome (nonparametric two-sided Wilcoxon test). c) Correlation between median πN/πS for inversion haplotypes and their frequency across 18 spruce bark beetle populations. Only inversion haplotypes that were present in more than four individuals and that included four or more 200-kb windows, and more than five genes were included in the mutation load analyses. Mykhailenko et al. GBE 12 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 inversion polymorphism in the spruce bark beetle. Likewise, the absence of an elevated mutation load in inverted re- gions does not support a role of associative overdominance either. However, since we only have genomic data from a single season, we cannot rule out that other mechanisms, such as negative frequency-dependent selection or antag- onistic pleiotropy, could maintain balanced inversion poly- morphisms. Additional temporal data spanning multiple years (generations) are needed to test whether temporally varying selection affects the frequencies of inversion haplotypes in the spruce bark beetle. Further research is thus essential to determine the importance of different me- chanisms in maintaining inversion polymorphisms within this species. Far-Reaching Consequences of Having an Inversion-Rich Genome The presence of multiple polymorphic inversions can have significant consequences both for the evolution of a species and for evolutionary inferences based on genome-wide polymorphism data. Importantly, polymorphic inversions are a reservoir of genetic variation that can facilitate adap- tation to rapidly changing environments. Indeed, several studies have shown that polymorphic inversions support ra- pid adaptation to changing climatic conditions (Rane et al. 2015; Kapun and Flatt 2019; McCulloch and Waters 2023) or adaptive tracking of fluctuating environments (Nunez et al. 2023). Spruce bark beetle populations are subjected to seasonal weather changes and a rapidly changing envir- onment due to strong anthropogenic pressures (Hofmann et al. 2024). Warmer weather and drought periods have been associated with an intensification of bark beetle out- breaks (Biedermann et al. 2019; Bentz et al. 2021; Hlásny et al. 2021), which may act as a strong selection factor with- in bark beetle populations. Whether inversions are involved in rapid adaptations in the spruce bark beetle is an open question that requires further investigation. Abundant polymorphic inversions within the genome can potentially bias inferences about, e.g. a species’ demo- graphic history and selection, if it is not properly accounted for. This is mainly due to the suppression of recombination within inversion regions and selection that influence inver- sion frequencies and variation patterns. It is well known that nonequilibrium demography and selection can leave similar genomic signatures. Traditionally, demographic analyses have used noncoding parts of the genome, based on the assumption that directional selection mostly affects protein-coding regions. However, there is growing evi- dence for the importance of background selection in shap- ing genome-wide diversity, and this is moving the field toward incorporating linked selection into inferences of demographic history (Li et al. 2012; Johri et al. 2021). We believe that new analytical approaches should also consider the potential influence of polymorphic inversion land- scapes, because variation patterns of inverted regions can be shaped by different types of balancing selection. In add- ition, genomic scans for selection in inversion-rich genomes may yield biased results because recombination is reduced within inversion regions. Importantly, the effect of reduced recombination may extend outside inversions (Adrion et al. 2020; Koury 2023; Li et al. 2023). Quantifying how inver- sions and inversion-rich genomes influence various evolu- tionary inference methods may become a standard approach to test the robustness of these methods (Novo et al. 2023). Conclusions Facilitated by advanced sequencing technologies and a re- cent focus on structural variation, evolutionary biology is fa- cing a new era of exciting discoveries that will deepen our understanding of genome complexity, genome-wide pat- terns of variation, and the major evolutionary forces re- sponsible for shaping these patterns (Lynch and Conery 2003; Wellenreuther et al. 2019; Weissensteiner et al. 2020; Sirén et al. 2021; Berdan et al. 2023; Wold et al. 2023). While extensive genome rearrangements are com- mon among species across the tree of life (Bhutkar et al. 2008; Escudero et al. 2023; Höök et al. 2023), and perhaps more so than previously thought (Lewin et al. 2024), accu- mulating genomic data suggests that such rearrangements are also common at the intraspecific level (Faria, Chaube, et al. 2019; Faria, Johannesson, et al. 2019; Todesco et al. 2020; Harringmeyer and Hoekstra 2022; Porubsky et al. 2022; Reeve et al. 2023). Here, we characterized one of the most complex polymorphic inversion landscapes described to date and found no support for any of several mechanisms that often are associated with the presence of polymorphic inversions (including directional selection, overdominance and associative overdominance, and strong divergent selection). These results suggest that inversions are either neutral or maintained by the combined action of multiple evolutionary forces and influence evolutionary processes in complex ways. Complete neutrality of large in- versions is unlikely (Connallon and Olito 2022), especially in species with large Ne, such as the spruce bark beetle, in which selection should be particularly efficient. Several studies have convincingly linked polymorphic inversions to adaptation (Lowry and Willis 2010; Twyford and Friedman 2015; Kirubakaran et al. 2016), but as more spe- cies are systematically screened for structural variation, data are accumulating on species in which (at least some) SVs show no clear adaptive advantage or disadvantage (this study and Harringmeyer and Hoekstra 2022). These find- ings raise questions about how prevalent inversions and other SVs are at different levels of the tree of life and, im- portantly, about their consequences for individual fitness Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 13 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 and the role they play in genome evolution and reorienta- tion of the evolutionary process (Marroni et al. 2014; Berden et al. 2023). Materials and Methods The Spruce Bark Beetle Genome and Its Quality The spruce bark beetle karyotype is n = 14 + XYp. The spe- cies genome assembly consists of 272 contigs, and the 36 largest contigs (>1 Mb) constitute the most of the assembly (78% of the 236.8-Mb-long assembly, N50 = 6.65 Mb; Powell et al. 2021). Telomeric motifs (sequences present at the end of chromosomes) were identified at the end of eight contigs (including the five largest contigs). BUSCO analysis indicated that 99.5% of the genes present in the insect database (insects_odb9) were also present among the predicted spruce bark beetle genes. To further investi- gate the quality of the spruce bark beetle assembly, we compared it to the chromosome-level assembly of the reference genome of the congener I. nitidus (NCBI: GCA_018691245.2; 1.5 to 2.3-Mya divergence from the spruce bark beetle; Wang et al. 2023). For clarity, we only ex- tracted the 36 and 16 longest contigs/pseudochromosomes from I. typographus and I. nitidus genomes, respectively, and used these in our synteny analysis. Syntenic blocks were de- tected using ntSynt version 1.0.1 (Coombe et al. 2024), with divergence set to 2.5% and using default values of the re- maining parameters. Links longer than 10 kb were processed using scripts provided with ntSynt and visualized using ggge- nomes version 1.0.0 (Hackl et al. 2024; https://github.com/ thackl/gggenomes) in R. SeqKit version 2.6.1 (Shen et al. 2024) was used for operations on fasta files, such as extrac- tion or reverse complement of sequences. Sampling of Beetles Adult spruce bark beetles were collected with pheromone- baited traps in the spring and summer of 2020. In total, we sampled 18 populations throughout Europe, using 13 to 14 individuals per locality (253 individuals in total) (Fig. 2; supplementary table S5, Supplementary Material online). Throughout the text, we use the term “population” to refer to an individual site or a collection of closely situated sites (within about 50 km). In Austria, we pooled individuals from three localities that were up to 120 km apart because of small sample sizes (five beetles or less per site). Populations from Fennoscandia are referred to as northern populations (or the northern group) and populations from central Europe are referred to as southern populations (or the southern group). Polish populations are considered sep- arately from other central European populations because downstream analysis revealed high admixture proportions from the northern group. Beetles were brought alive to the laboratory, starved on a paper diet for several days, dissected, sexed based on genitalia morphology, and subjected to DNA extraction (described below). DNA Extraction and Genome Resequencing DNA was extracted from the whole body of dissected beetles using the Wizard Genomic DNA Purification Kit (Promega). The concentration of extracted DNA was estimated using a Qubit fluorometer (Thermo Fisher Scientific). Genomic librar- ies were prepared with NEBNext Ultra II FS DNA Library Prep with Beads (New England Biolabs), with single indexes. Individual libraries were combined into three pools and 2 × 150 bp paired-end sequenced in three lanes of a S4 flowcell using the NovaSeq 6000 instrument and v1 sequen- cing chemistry (Illumina Inc.). Sequencing was done by the National Genomics Infrastructure, SNP&SEQ Technology Platform (Uppsala, Sweden). To assess the overall genotyp- ing error, we prepared and sequenced duplicate libraries for nine individuals. Data Preparation and Filtering Details of raw data processing and filtering are described in the Supplementary Material. In short, raw reads were mapped to the reference genome (Powell et al. 2021) using Bowtie 2 version 2.4.2 (Langmead and Salzberg 2012). Duplicated reads were removed using Picard MarkDuplicates version 2.24.1 (https://broadinstitute. github.io/picard/). To detect and correct systematic errors in base quality scores, recalibration was done using the GATK version 4.1.9.0, BaseRecalibrator, and ApplyBQSR (McKenna et al. 2010; Depristo et al. 2011). Variant calling and genotyping were done using GATK HaplotypeCaller, CombineGVCFs, and GenotypeGVCFs. GATK VariantRecalibrator and ApplyVQSR were used to calculate and filter (by variant) quality score log-odds (VQSLOD). Bcftools version 1.11 (Danecek et al. 2021) was used to re- move insertions and deletions (indels) as well as any poly- morphisms of five bases upstream and downstream. GATK VariantFiltration was applied to mask all genotypes with low sequencing depth or low genotype quality (McKenna et al. 2010; Depristo et al. 2011). Variants that were not biallelic SNPs or did not meet the recommended hard filtering thresholds (GATK Team; see Supplementary Material) were filtered out. To filter out polymorphisms that could come from duplicated regions, we removed var- iants located within repeat-masked regions of the genome (Powell et al. 2021), variants with excessive overall coverage (mean + 1 SD), and variants with heterozygote excess. Variants for which genotypes could be detected in less than half of the individuals were also removed. We used PLINK version 1.90b6.18 (Purcell et al. 2007) to detect sam- ple contaminations, swaps and duplications, and familial relationships (e.g. sibling pairs present in the data), which might bias downstream analyses. Individuals with excessive Mykhailenko et al. GBE 14 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 coverage were removed, as these could be a result of hu- man errors during library preparation or pooling. We only analyzed contigs longer than 1 Mb in downstream analysis. These constituted 78% of the genome assembly, i.e. a total of 186 Mb. Since part of IpsContig33 had high similarity to mtDNA, this contig was not included in the downstream analysis. Genotyping error was assessed using GATK Genotype Concordance. Genome-Wide Genetic Variation and Its Geographic Structuring Genome-wide genetic structuring was explored by PCA using PLINK. The most likely number of genetic clusters and admixture proportions was estimated using NGSadmix (Skotte et al. 2013). The analysis was run for five different K-values (1 to 5; ten replicates per K-value), using a minor allele frequency (MAF) filter of 0.05 and 10,000 iterations, and a SNP data set that was pruned for LD using PLINK (–indep-pairwise 50 10 0.1 option). To choose the most likely number of genetic clusters, the results were examined using CLUMPAK (http://clumpak.tau.ac.il/index.html). To examine the influence of inversions on genetic clustering and to facili- tate inversion genotyping, NGSadmix was run separately for (i) all autosomal contigs without potential inversions and (ii) each potential inversion. To assess genetic differentiation among different spruce bark beetle populations Weir and Cockerham (1984) FST was estimated using VCFtools version 0.1.16 (Danecek et al. 2011). FST was calculated between population groups identified by NGSadmix, as well as among all population pairs. Additionally, we summarized FST values in 100-kb overlapping windows (using 20-kb steps) along the contigs. Window-based analyses were done for all contigs and all possible population pairs. In addition, absolute sequence divergence (dxy), nucleotide diversity, and Tajima’s D statis- tic were estimated and summarized for 100-kb nonoverlap- ping windows using ANGSD version 0.935-44-g02a07fc (Korneliussen et al. 2014). These three statistics were calcu- lated for each population separately and were based on a maximum likelihood estimate of the folded site frequency spectrum. We excluded sites with mapping quality below 30 phred, base quality score below 20, and coverage less than three times the population sample size and more than three times the average coverage, following the ap- proach used in Delmore et al. (2018). The ANGSD calcula- tions were based on allele frequencies estimated from genotype likelihoods (Li 2011) and ngsPopGen scripts (https://github.com/mfumagalli/ngsPopGen). Identification of Inversions, Their Geographic Distribution, and Variation Patterns To identify inversions, we followed the standard population genetic approach (see e.g. Huang et al. 2020; Mérot et al. 2021; Reeve et al. 2023). Potential chromosomal inversion regions were identified based on (i) per contig PCAs, (ii) local PCAs, (iii) patterns of heterozygosity, and (iv) LD clustering. Per contig PCAs were performed as a first screen for inversions using PLINK (Purcell et al. 2007). Local PCAs were performed to identify the approximate pos- ition of inversions using the lostruct R package (Li and Ralph 2019). Analysis was performed using two inde- pendent approaches: localPCA analysis using all analyzed contigs and localPCA analysis for each contig independ- ently. The first analysis was run following the approach described in Huang et al. (2020). Windows were consid- ered outliers if their multidimensional scaling (MDS) scores were more than 2 SD greater than the mean across all windows. The maximum number of windows between what was considered separate inversions was set to 20. In the second analysis, separate contigs were subjected to the k-means clustering algorithm to define groups of windows that formed a single inversion. k-means clustering was performed on the first MDS scores. Different k-values (number of window clusters) were tested for each contig so that structurally different parts of the contig were isolated into separate clusters. For both approaches, different window sizes (1, 10, and 100 kb) were tested. LD among SNPs (thinned by selecting one SNP every 10 kb; MAF > 5%) was calculated for each contig using PLINK. We considered a genomic region to be an inversion region if (i) local PCA analysis identified the region as an outlier, (ii) the region exhibited high LD (most SNPs having r2 > 0.4), and (iii) PCA performed on SNPs from this region separated individuals into three distinct groups with hetero- zygosity patterns matching the expectation of the middle group of the PCA presenting higher heterozygosity. Genotyping of individual beetles with respect to the in- version haplotypes they carried was done based on NGSadmix clustering (with k = 2 inversion heterozygotes having mixed ancestry in approximate 50/50 proportions; supplementary fig. S13, Supplementary Material online). In a few more complex cases (putative overlapping inver- sions and double-crossover events) genotyping was done by clustering individuals based on the PCA groups visible along PC1 and PC2 (for more details, see supplementary fig. S3, Supplementary Material online and related text). Contigs with <10,000 SNPs provided ambiguous results and were excluded from genotyping. Approximate inversion boundaries were defined based on local PCAs and detection of sharp borders in LD clus- ters. The population genetic approach used in this study is a standard approach to detect putative large poly- morphic inversions but is not able to provide details on in- version breakpoints. As further evidence for inversions, we estimated popula- tion recombination rates (rho) for each putative inversion Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 15 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 region using individuals with a particular inversion geno- type (separately for homozygous and heterozygous indivi- duals). We followed the approach used in Jones et al. (2019) and sampled 20 individuals from each genotype group. For genotype groups with less than ten individuals, rho was not calculated, but rho was calculated for groups with fewer than 20 but more than ten individuals. Watterson theta estimates were used to create a custom likelihood lookup table using the Ldhat 2.2 program complete (McVean et al. 2002). The interval program was used to estimate the population recombination rate across investigated inversion (and adjacent regions in the same contig). The interval algorithm was run for 2 million itera- tions and the chain was sampled every 10,000 iterations with a burn-in of 100,000 generations (using the Ldhat package’s stat program, block penalty = 5). Population re- combination rates were summarized in nonoverlapping 100-kb windows using a custom perl script. Inversion genotype and haplotype frequencies were calculated using an in-house R-script. Frequencies were calculated (i) within each population, (ii) within the south- ern and northern groups, and (ii) for all sampled popula- tions combined. Deviations from HW equilibrium were estimated for all three data sets. Inversions with only two haplotypes were tested for HW deviations using Fisher’s exact tests. Inversions with more than two haplo- types (including recombinant haplotypes between two in- version arrangements) were tested using a permutation test, and sex chromosome inversions were tested as described in Graffelman and Weir (2018). All tests were done using the R package HardyWeinberg (Graffelman 2015). To investigate if inversion haplotypes differed in frequency along environmental gradients, Pearson correl- ation between inversion haplotype frequencies and latitude/longitude was calculated. To assess whether cor- relations were stronger than expected based on the col- linear part of the genome, we performed a permutation test. We generated sets of 1,000 randomly selected SNPs with MAFs similar to (±0.05) the frequencies of in- versions showing significant correlations. We first calcu- lated the MAF of each SNP for each population and then calculated Pearson’s correlation between SNP fre- quencies and population latitudinal coordinates. The re- sulting r2 distribution was used to determine the number of loci outside inversions that exceeded the r2 ob- tained from the correlation between inversion haplotype frequency and latitude and to calculate the P-value of the test. To assess levels of genetic differentiation between inversion haplotypes, FST and dxy between alternative in- version haplotypes (AA and BB) were estimated following the approach described above. Deletions present in alter- native inversion haplotypes were not included in dxy calculations but were summarized separately using custom scripts. Testing for Enrichment of GO Terms in Inversions We used gene annotations from Powell et al. (2021) and tested for overrepresentation of GO terms among genes lo- cated within inversions using R package topGO (Alexa and Rahnenführer 2024), applying Fisher’s exact test and the “weight01” algorithm (Alexa et al. 2006) to deal with the GO graph structure; only GO categories with at least ten members among SNP-associated genes were considered. Inversion Age Estimation Absolute sequence divergence between alternative inver- sion haplotypes was used to calculate the approximate time of divergence of inverted and noninverted haplotypes. We used the equation T = dxy/2μ, where T is the divergence time in generations, μ is the mutation rate per site per gen- eration, and dxy is a mean dxy calculated based on per SNP values estimated in ANGSD. Since mutation rates of the spruce bark beetle are unknown, we used a range of muta- tion rate estimates available for three diploid, sexually re- producing insects: Drosophila melanogaster (Krasovec 2021), Heliconius melpomene (Keightley et al. 2015), and Chironomus riparius (Oppold and Pfenninger 2017). These per-generation and per-basepair mutation rate esti- mates varied from 2.1 to 11.7 × 10−9. This approach could only give rough inversion age estimates due to the uncer- tainty of the mutation rate estimates, probable intraspecific variation in mutation rate (Krasovec 2021), and a (likely) substantial influence of selection and gene flux (Charlesworth 2023). Mutation Load Estimation To estimate mutation load, we calculated the ratio of nu- cleotide diversity at nonsynonymous sites (πN) versus syn- onymous sites (πS). Mutation load (πN/πS) was calculated separately for each inversion homozygote and for the col- linear part of the genome (for all individuals). We computed nucleotide diversity for each site using SNPGenie (Nelson et al. 2015). The πN/πS ratio was estimated in 200-kb win- dows using an in-house R script. To account for the fact that inversions can greatly suppress recombination in sur- rounding parts of the genome (Koury 2023), the collinear part of the genome was divided into two groups: (i) a group including all collinear 200-kb windows outside inversions and (ii) a group including all collinear windows outside in- versions but excluding windows that came from contigs with inversions (so-called strict filtering). Both collinear data sets were used to test for overall differences in muta- tion load between inversions and the collinear part of the genome (using two-sided t-test). One-sided t-tests were used to test whether minor (less frequent) homokaryotypes had higher mutation loads than major (more frequent) homokaryotypes. Homokaryotypes that contained fewer Mykhailenko et al. GBE 16 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 than four 200-kb windows and were present in few indivi- duals (two thresholds were tested: <4 and <10 individuals) were excluded from the analysis. Additionally, windows with a small number of genes were excluded (two thresh- olds were tested: <5 and <10 genes). Haplotype structuring, i.e. the presence of two or more distinct subhaplotypes among inversion homozygote hap- lotypes, can prevent fitness degeneration on one or both in- version haplotypes by carrying partially complementary sets of deleterious recessive alleles (Charlesworth and Charlesworth 1997; Berdan et al. 2021). To check if any in- version homozygotes exhibited haplotype structuring, we first phased the data using Beagle 5.2 (using default set- tings; Browning and Browning 2007). Next, for each in- verted region and homokaryotype, we filtered out all variants with MAF < 0.1 and used PGDSpider 2.1.1.0 (Lischer and Excoffier 2012) to convert VCFs to full-length sequences. Finally, we constructed neighbor-joining trees for alleles within haplotypes using MEGA7 and investigated if topologies showed the presence of distinct clusters within homokaryotype groups (Kumar et al. 2016). Phenotype–Genotype Associations To test whether inversion polymorphisms were associated with diapause phenotypes, we resequenced the whole genome of 18 female spruce bark beetle individuals with defined individual diapause phenotypes (ten beetles ex- pressing facultative diapause and eight beetles expressing obligate diapause). These individuals came from a spruce bark beetle diapause study by Schebeck et al. (2022). In brief, diapause phenotypes were determined under con- trolled conditions by exposing beetles from three locations (northern Scandinavia, Central Europe low-elevation site, and Central Europe high-elevation site) to diapause- inducing or diapause-averting photoperiodic conditions. Diapause expression was assessed by studying gonad devel- opment and timing of emergence from experimental logs (two reliable indicators of diapause expression in the spruce bark beetle). The individuals with the best DNA extractions were selected for sequencing using the DNBseq platform (BGI Tech Solutions, Poland) to a mean coverage of 20×. Data were processed in the same way as described above (raw data processing and inversion genotyping). To test for differentiation in inversion haplotype frequencies be- tween diapause phenotypes, we ran an exact G-test using Genepop 4.1.2 (Raymond and Rousset 1995). FST between individuals expressing facultative and obligate diapause was estimated using VCFtools version 0.1.16 (Danecek et al. 2011) and summarized in 100-kb overlapping windows (20-kb steps). To check if spruce bark beetle OR genes were associated with inversions, we examined 73 OR genes recently anno- tated by Yuvaraj et al. (2021). We used OR sequences from Yuvaraj et al. (2021) and mapped them to the spruce bark beetle reference genome using minimap2 (Li 2018). We located 71 out of 73 ORs in the 170 Mb covered by the 26 contigs we analyzed (since contigs with <10,000 SNPs were excluded from inversion detection and genotyp- ing). Three ORs mapped to more than one contig. ItypOR9 and ItypOR58 mapped to the end of IpsContig16 and 23, suggesting possible assembly error and duplication of end-of-contig sequences, which are difficult to assemble. ItypOR59NTE mapped to three nearby (within 7 kb) loca- tions on IpsContig6, suggesting either assembly error or recent duplication. For these three ORs, we used one ran- domly chosen location in downstream analyses. To test if inversions were significantly enriched in OR genes, we ran a permutation test (10,000 iterations; permutating inver- sion locations over the 170-Mb genome consisting of 26 contigs, while keeping OR locations fixed). We tested how many permutations yielded a higher (or equal) number of ORs located within inversions than that observed in the real data. To check if alternative inversion arrangements harbored different numbers of OR genes (e.g. if one ar- rangement carried a deletion), we compared sequence coverage within ORs in individual beetles identified as inver- sion homozygotes. Additionally, we counted nonsynon- ymous segregating variants present in ORs situated within inversions and checked if any of these SNPs were highly di- verged between inversion haplotypes (by calculating dxy). Finally, to investigate whether natural selection was in- volved in shaping nonsynonymous variation patterns in ORs between inversion haplotypes, we calculated dN/dS ra- tio for ORs with multiple highly diverged nonsynonymous variants and performed MacDonald–Kreitman test for ORs with dN/dS > 1. dxy was calculated using GATK VariantsToTable outputs and custom perl scripts, dN/dS was calculated using the pairwise distances option in MEGA11 (syn–nonsynonymous substitution model, Kumar method), and the MacDonald–Kreitman test was performed in DnaSP v6 (Rozas et al. 2017). Genotype–environment Association Genotype–environment association (GEA) analyses were done to test if allele frequency changes in SNPs were asso- ciated with the beetle populations’ local environment. The analysis was performed using (i) all SNPs present in the ana- lyzed part of the genome and (ii) excluding SNPs present within inversion regions and coding each inversion as single SNPs. Many different environmental variables were sum- marized along two principal components (see below). To control for the confounding effect of the overall genetic dif- ferentiation, we used latent factor mixed models (LFMMs; Caye et al. 2019) as implemented in the lfmm2() function from the R package LEA (Gain and François 2021). We used only SNPs with <20% missing data, MAF ≥ 0.1, and Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 17 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 that occurred in individuals with <30% missing data. Because LFMM cannot handle missing data, we imputed missing genotypes using impute() from LEA. We ran LFMM for (i) all SNPs that passed the filters described above and (ii) for a data set where each inversion was represented as a single “SNP” inversion genotype. We used five latent factors (k = 5) in lfmm2(). P-values were calculated using lfmm2.test() from LEA and false discovery rate (FDR) cor- rected using the p.adjust() R function with method = “fdr” (Benjamini and Hochberg 1995). Each beetle population’s local environment was charac- terized according to climate and land cover data. We used all 19 bioclimatic variables from WorldClim version 2.1. (Fick and Hijmans 2017) with a resolution of ∼1 km2. These bioclimatic variables included annual mean tempera- ture, mean diurnal range, isothermality, temperature sea- sonality, maximum temperature of warmest month, minimum temperature of coldest month, temperature an- nual range, mean temperature of wettest quarter, mean temperature of driest quarter, mean temperature of warm- est quarter, mean temperature of coldest quarter, annual precipitation, precipitation of wettest month, precipitation of driest month, precipitation seasonality, precipitation of wettest quarter, precipitation of driest quarter, precipita- tion of warmest quarter, and precipitation of coldest quar- ter. These variables are averaged over the years 1970 to 2000. Proportions of forest, cropland, and built-up areas across the geographical range of the spruce bark beetle were downloaded from https://lcviewer.vito.be/ for 2015 with a spatial resolution of ∼100 m2 (Buchhorn et al. 2019). These global land cover maps are part of the Copernicus Land Service, derived from PROBA-V satellite observations, and have an accuracy of 80% as measured by the CEOS land product validation subgroup. We also in- cluded the proportion of land area covered by spruce trees (genus Picea) at a resolution of ∼1 km2, as obtained by Brus et al. (2012) using a statistical mapping approach. All biocli- matic variables were reprojected to a final resolution of ∼1 km2 using the Lambert azimuthal equal area method. We then extracted mean values for all environmental vari- ables (supplementary fig. S14, Supplementary Material on- line) within a 50-km radius from each population location using the R package terra (Hijmans 2022). Finally, PCA was used to summarize the multiscale environmental vari- ation among populations. The first two PCA components (PC1 and PC2) explained 25% and 23% of the environ- mental variation, respectively, and were used as the final in- put for the GEA analyses. PC1 represented environmental variation mainly related to latitude, with northern popula- tions showing higher values indicative of higher tempera- ture variation and lower temperatures during the coldest months. PC2 represented environmental variation mostly related to temperature and amount of cropland, with high- er values representing localities with higher temperatures during the warmest months and a higher proportion of cropland (and conversely less forest cover and spruce; supplementary fig. S15, Supplementary Material online). Additionally, we identified genes closest to SNPs that showed significant association with each PC. Supplementary Material Supplementary material is available at Genome Biology and Evolution online. Acknowledgments We thank A. Hietala, E. Stengel, K. Zub, Å. Lindelöw, O. Langvall, M. Holmlund, E. Kristensen, U. Johansson, R. Modlinger, J. Reisenberger, K. Szreder, W. Skowroński L. Stanecki, and M. Ahlström for the help in sampling spruce bark beetles across Europe. T. Mokrzycki helped with beetle sexing and identification. We thank members of the Genomics and Experimental Evolution Group at Jagiellonian University for their help in improving this manuscript. We thank Tomasz Gaczorek for the help in writing scripts and optimizing data analysis. Author Contributions A.M., P.Z., and K.N.-B. conceived the study, performed the main analyses, and wrote the manuscript. A.B. mana- ged sample shipment, sexed beetles, isolated DNA, and prepared samples for sequencing. F.S., P.K., and M.N.A. provided an unpublished version of the spruce bark beetle genome and insights on the species’ ecology and sensory biology. M.M., J.M., Z.B., M.S., P.K., and C.S. organized sampling and provided beetles for analysis and insights on sampled populations. B.A. and W.B. performed the GEA analysis. Z.N. analyzed the diapause data. J.M. phased the data. M.S. provided the diapause samples. J.Z performed the I. nitidus analysis. W.B. helped in data interpretation and provided feedback on all manuscript versions. All authors read and approved the final manuscript. Funding This work was funded by a Polish National Science Center 2018/30/E/NZ8/00105 grant to K.N.-B and the Foundation in Memory of Oscar and Lili Lamm to M.N.A. Data Availability All DNA sequences have been deposited to the National Center for Biotechnology Information Sequence Read Archive under the BioProject ID PRJNA1013983. Scripts can be found https://github.com/AnstMykh/Ips_ inversions_MS/tree/main. Mykhailenko et al. GBE 18 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Literature Cited Adrion JR, Galloway JG, Kern AD. Predicting the landscape of recom- bination using deep learning. Mol Biol Evol. 2020:37(6): 1790–1808. https://doi.org/10.1093/molbev/msaa038. Akopyan M, Tigano A, Jacobs A, Wilder AP, Baumann H, Therkildsen NO. Comparative linkage mapping uncovers recombination sup- pression across massive chromosomal inversions associated with local adaptation in Atlantic silversides. Mol Ecol. 2022:31(12): 3323–3341. https://doi.org/10.1111/mec.16472. Alexa A, Rahnenführer J. topGO: enrichment analysis for gene ontol- ogy. https://bioconductor.org/packages/topGO. Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006:22(13):1600–1607. https://doi. org/10.1093/bioinformatics/btl140. Anderbrant O. Reemergence and second brood in the bark beetle Ips typographus. Holarct Ecol. 1989:12:494–500. https://doi.org/10. 1111/j.1600-0587.1989.tb00927.x. Anderbrant O, Lofqvist J. Relation between first and second brood pro- duction in the bark beetle Ips typographus (Scolytidae). Oikos. 1988:53(3):357–365. https://doi.org/10.2307/3565536. Andersson MN, Larsson MC, Schlyter F. Specificity and redundancy in the olfactory system of the bark beetle Ips typographus: single-cell responses to ecologically relevant odors. J Insect Physiol. 2009: 55(6):556–567. https://doi.org/10.1016/j.jinsphys.2009.01.018. Annila E. Influence of temperature upon the development and the voltinism of Ips typographus L. (Coleoptera, Scolytidae). Ann Zool Fennici. 1969:6:161–208. https://www.jstor.org/stable/2373 1366. Ayala D, Acevedo P, Pombi M, Dia I, Boccolini D, Costantini C, Simard F, Fontenille D. Chromosome inversions and ecological plasticity in the main African malaria mosquitoes. Evolution. 2017:71(3): 686–701. https://doi.org/10.1111/evo.13176. Ayala D, Ullastres A, González J. Adaptation through chromosomal in- versions in Anopheles. Front Genet. 2014:5:129. https://doi.org/ 10.3389/fgene.2014.00129. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a prac- tical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995:57(1):289–300. https://doi.org/10.1111/j. 2517-6161.1995.tb02031.x. Bentz BJ, Hansen EM, Davenport M, Soderberg D. Complexities in pre- dicting mountain pine beetle and spruce beetle response to cli- mate change. In: Gandhi KJK, Hofstetter RW, editors. Bark beetle management, ecology and climate change. Cambridge, MA: Academic Press; 2021. p. 31–54. Berdan EL, Barton NH, Butlin R, Charlesworth B, Faria R, Fragata I, Gilbert KJ, Jay P, Kapun M, Lotterhos KE. How chromosomal inver- sions reorient the evolutionary process. J Evol Biol. 2023:36: 1761–1782. https://doi.org/10.1111/jeb.14242. Berdan EL, Barton NH, Butlin R, Charlesworth B, Faria R, Fragata I, Gilbert KJ, Jay P, Kapun M, Lotterhos KE, et al. How chromosomal inversions reorient the evolutionary process. J Evol Biol. 2023:36(12):1761–1782. https://doi.org/10.1111/jeb.14242. Berdan EL, Blanckaert A, Butlin RK, Bank C. Deleterious mutation accu- mulation and the long-term fate of chromosomal inversions. PLoS Genet. 2021:17(3):e1009411. https://doi.org/10.1371/journal. pgen.1009411. Bertheau C, Schuler H, Arthofer W, Avtzis DN, Mayer F, Krumböck S, Moodley Y, Stauffer C. Divergent evolutionary histories of two sympatric spruce bark beetle species. Mol Ecol. 2013:22(12): 3318–3332. https://doi.org/10.1111/mec.12296. Bhutkar A, Schaeffer SW, Russo SM, Xu M, Smith TF, Gelbart WM. Chromosomal rearrangement inferred from comparisons of 12 Drosophila genomes. Genetics. 2008:179(3):1657–1680. https:// doi.org/10.1534/genetics.107.086108. Biedermann PHW, Müller J, Grégoire JC, Gruppe A, Hagge J, Hammerbacher A, Hofstetter RW, Kandasamy D, Kolarik M, Kostovcik M, et al. Bark beetle population dynamics in the Anthropocene: challenges and solutions. Trends Ecol Evol. 2019:34(10):914–924. https://doi.org/10.1016/j.tree.2019. 06.002. Brelsford A, Purcell J, Avril A, Tran Van P, Zhang J, Brütsch T, Sundström L, Helanterä H, Chapuisat M. An ancient and eroded social supergene is widespread across formica ants. Curr Biol. 2020:30(2):304–311.e4. https://doi.org/10.1016/j.cub.2019. 11.032. Browning SR, Browning BL. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. Am J Hum Genet. 2007:81(5):1084–1097. https://doi.org/10.1086/ 521987. Brus DJ, Hengeveld GM, Walvoort DJJ, Goedhart PW, Heidema AH, Nabuurs GJ, Gunia K. Statistical mapping of tree species over Europe. Eur J For Res. 2012:131(1):145–157. https://doi.org/10. 1007/s10342-011-0513-5. Buchhorn M, Smets B, Bertels L, Lesiv M, Tsendbazar NE, Herold M, Fritz S. Copernicus global land service: land cover 100m: epoch 2015: globe. Version V2. 0.2. 2019. Caye K, Jumentier B, Lepeule J, François O. LFMM 2: fast and accurate inference of gene-environment associations in genome-wide stud- ies. Mol Biol Evol. 2019:36(4):852–860. https://doi.org/10.1093/ molbev/msz008. Charlesworth B. The effects of inversion polymorphisms on patterns of neutral genetic diversity. Genetics. 2023:224(4):1–21. https://doi. org/10.1093/genetics/iyad116. Charlesworth B, Charlesworth D. Rapid fixation of deleterious alleles can be caused by Muller’s ratchet. Genet Res. 1997:70(1):63–73. https://doi.org/10.1017/S0016672397002899. Cheng C, Tan JC, Hahn MW, Besansky NJ. Systems genetic analysis of inversion polymorphisms in the malaria mosquito Anopheles gam- biae. Proc Natl Acad Sci U S A. 2018:115(30):E7005–E7014. https://doi.org/10.1073/pnas.1806760115. Christmas MJ, Wallberg A, Bunikis I, Olsson A, Wallerman O, Webster MT. Chromosomal inversions associated with environmental adap- tation in honeybees. Mol Ecol. 2019:28(6):1358–1374. https://doi. org/10.1111/mec.14944. Connallon T, Clark AG. Balancing selection in species with separate sexes: insights from Fisher’s geometric model. Genetics. 2014:197(3):991–1006. https://doi.org/10.1534/genetics.114. 165605. Connallon T, Olito C. Natural selection and the distribution of chromo- somal inversion lengths. Mol Ecol. 2022:31(13):3627–3641. https://doi.org/10.1111/mec.16091. Coombe L, Kazemi P, Wong J, Birol I, Warren RL. Multi-genome synteny detection using minimizer graph mappings. bioRxiv 579356. https:// doi.org/10.1101/2024.02.07.579356, 13 February 2024, preprint: not peer reviewed. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. The variant call format and VCFtools. Bioinformatics. 2011:27(15): 2156–2158. https://doi.org/10.1093/bioinformatics/btr330. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021:10(2): giab008. https://doi.org/10.1093/gigascience/giab008. Delmore KE, Lugo Ramos JS, Van Doren BM, Lundberg M, Bensch S, Irwin DE, Liedvogel M. Comparative analysis examining patterns Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 19 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 of genomic differentiation across multiple episodes of population divergence in birds. Evol Lett. 2018:2(2):76–87. https://doi.org/10. 1002/evl3.46. Denlinger D. Insect diapause. Cambridge, England: Cambridge University Press; 2022. Depristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, Del Angel G, Rivas MA, Hanna M, et al. A frame- work for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011:43(5):491–501. https:// doi.org/10.1038/ng.806. Dobzhansky T. Genetics of the evolutionary process. New York, USA: Columbia University Press; 1970. Dobzhansky T, Sturtevant AH. Inversions in the chromosomes of Drosophila pseudoobscura. Genetics. 1938:23(1):28–64. https:// doi.org/10.1093/genetics/23.1.28. Doležal P, Sehnal F. Effects of photoperiod and temperature on the de- velopment and diapause of the bark beetle Ips typographus. J Appl Entomol. 2007:131(3):165–173. https://doi.org/10.1111/j.1439- 0418.2006.01123.x. Dowle EJ, Powell THQ, Doellman MM, Meyers PJ, Calvert MB, Walden KKO, Robertson HM, Berlocher SH, Feder JL, Hahn DA, et al. Genome-wide variation and transcriptional changes in diverse de- velopmental processes underlie the rapid evolution of seasonal adaptation. Proc Natl Acad Sci U S A. 2020:117(38):23960–23969. https://doi.org/10.1073/pnas.2002357117. Dworschak K, Gruppe A, Schopf R. Survivability and post-diapause fitness in a scolytid beetle as a function of overwintering developmental stage and the implications for population dynam- ics. Ecol Entomol. 2014:39(4):519–526. https://doi.org/10.1111/ een.12127. Ellerstrand SJ, Choudhury S, Svensson K, Andersson MN, Kirkeby C, Powell D, Schlyter F, Jönsson AM, Brydegaard M, Hansson B, et al. Weak population genetic structure in Eurasian spruce bark beetle over large regional scales in Sweden. Ecol Evol. 2022:12(7):e9078. https://doi.org/10.1002/ece3.9078. Escudero M, Marques A, Lucek K, Hipp AL. Genomic hotspots of chromosome rearrangements explain conserved synteny despite high rates of chromosome evolution in a holocentric lineage. Mol Ecol. 2023:33:e17086. https://doi.org/10.1111/mec.17086. Faria R, Chaube P, Morales HE, Larsson T, Lemmon AR, Lemmon EM, Rafajlovic´ M, Panova M, Ravinet M, Johannesson K, et al. Multiple chromosomal rearrangements in a hybrid zone between Littorina saxatilis ecotypes. Mol Ecol. 2019:28(6):1375–1393. https://doi. org/10.1111/mec.14972. Faria R, Johannesson K, Butlin RK, Westram AM. Evolving inversions. Trends Ecol Evol. 2019:34(3):239–248. https://doi.org/10.1016/j. tree.2018.12.005. Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017:37(12): 4302–4315. https://doi.org/10.1002/joc.5086. Fuller ZL, Haynes GD, Richards S, Schaeffer SW. Genomics of natural populations: how differentially expressed genes shape the evolu- tion of chromosomal inversions in Drosophila pseudoobscura. Genetics. 2016:204(1):287–301. https://doi.org/10.1534/genetics. 116.191429. Fuller ZL, Haynes GD, Richards S, Schaeffer SW. Genomics of natural populations: evolutionary forces that establish and maintain gene arrangements in Drosophila pseudoobscura. Mol Ecol. 2017:26(23):6539–6562. https://doi.org/10.1111/mec.14381. Fuller ZL, Koury SA, Phadnis N, Schaeffer SW. How chromosomal rear- rangements shape adaptation and speciation: case studies in Drosophila pseudoobscura and its sibling species Drosophila persimilis. Mol Ecol. 2019:28(6):1283–1301. https://doi.org/10. 1111/mec.14923. Gain C, François O. LEA 3: factor models in population genetics and ecological genomics with R. Mol Ecol Resour. 2021:21(8): 2738–2748. https://doi.org/10.1111/1755-0998.13366. Graffelman J. Exploring diallelic genetic markers: the HardyWeinberg package. J Stat Softw. 2015:64(3):1–23. https://doi.org/10. 18637/jss.v064.i03. Graffelman J, Weir BS. Multi-allelic exact tests for Hardy–Weinberg equilibrium that account for gender. Mol Ecol Resour. 2018:18(3):461–473. https://doi.org/10.1111/1755-0998.12748. Guerrero RF, Rousset F, Kirkpatrick M. Coalescent patterns for chromosomal inversions in divergent populations. Philos Trans R Soc Lond B Biol Sci. 2012:367(1587):430–438. https://doi.org/ 10.1098/rstb.2011.0246. Gutiérrez-Valencia J, Hughes PW, Berdan EL, Slotte T. The genomic architecture and evolutionary fates of supergenes. Genome Biol Evol. 2021:13(5):evab057. https://doi.org/10.1093/gbe/evab057. Hackl T, Ankenbrand MJ, van Adrichem B, Haslinger K. gggenomes: a grammar of graphics for comparative genomics. 2024. [Accessed 2024 Sep 9]. Available from: https://github.com/thackl/ gggenomes. Harringmeyer OS, Hoekstra HE. Chromosomal inversion polymorph- isms shape the genomic landscape of deer mice. Nat Ecol Evol. 2022:6(12):1965–1979. https://doi.org/10.1038/s41559- 022-01890-0. Hijmans R. terra: spatial data analysis. R package version 1.7-18. 2022. [accessed 2022 Dec]. Available from: https://CRAN.R-project.org/ package=terra. Hlásny T, König L, Krokene P, Lindner M, Montagné-Huck C, Müller J, Qin H, Raffa KF, Schelhaas MJ, Svoboda M, et al. Bark beetle out- breaks in Europe: state of knowledge and ways forward for man- agement. Curr For Reports. 2021:7:138–165. https://doi.org/10. 1007/s40725-021-00142-x. Hlásny T, Krokene P, Liebhold A, Montagné-Huck C, Müller J, Qin H, Raffa K, Schelhaas M-J, Seidl R, Svoboda M. Living with bark bee- tles: impacts, outlook and management options. From Science to Policy 8. European Forest Institute; 2019. Hofmann S, Schebeck M, Kautz M. Diurnal temperature fluctuations improve predictions of developmental rates in the spruce bark bee- tle Ips typographus. J Pest Sci. 2024:97(4):1839–1852. https://doi. org/10.1007/s10340-024-01758-1. Höök L, Näsvall K, Vila R, Wiklund C, Backström N. High-density link- age maps and chromosome level genome assemblies unveil direc- tion and frequency of extensive structural rearrangements in wood white butterflies (Leptidea spp.). Chromosome Res. 2023:31(1):2. https://doi.org/10.1007/s10577-023-09713-z. Hou XQ, Yuvaraj JK, Roberts RE, Zhang DD, Unelius CR, Löfstedt C, Andersson MN. Functional evolution of a bark beetle odorant re- ceptor clade detecting monoterpenoids of different ecological origins. Mol Biol Evol. 2021:38(11):4934–4947. https://doi.org/ 10.1093/molbev/msab218. Huang K, Andrew RL, Owens GL, Ostevik KL, Rieseberg LH. Multiple chromosomal inversions contribute to adaptive divergence of a dune sunflower ecotype. Mol Ecol. 2020:29(14):2535–2549. https://doi.org/10.1111/mec.15428. Huang K, Ostevik KL, Elphinstone C, Todesco M, Bercovich N, Owens GL, Rieseberg LH. Mutation load in sunflower inversions is nega- tively correlated with inversion heterozygosity. Mol Biol Evol. 2022:39(5):msac101. https://doi.org/10.1093/molbev/msac101. Jay P, Chouteau M, Whibley A, Bastide H, Parrinello H, Llaurens V, Joron M. Mutation load at a mimicry supergene sheds new light on the evolution of inversion polymorphisms. Nat Genet. 2021:53(3):288–293. https://doi.org/10.1038/s41588-020-00771-1. Jay P, Tezenas E, Véber A, Giraud T. Sheltering of deleterious mutations explains the stepwise extension of recombination suppression on Mykhailenko et al. GBE 20 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 sex chromosomes and other supergenes. PLoS Biol. 2022:20(7): e3001698. https://doi.org/10.1371/journal.pbio.3001698. Johri P, Riall K, Becher H, Excoffier L, Charlesworth B, Jensen JD. The impact of purifying and background selection on the infer- ence of population history: problems and prospects. Mol Biol Evol. 2021:38(7):2986–3003. https://doi.org/10.1093/molbev/ msab050. Jones JC, Wallberg A, Christmas MJ, Kapheim KM, Webster MT. Extreme differences in recombination rate between the genomes of a solitary and a social bee. Mol Biol Evol. 2019:36(10): 2277–2291. https://doi.org/10.1093/molbev/msz130. Jönsson AM, Harding S, Bärring L, Ravn HP. Impact of climate change on the population dynamics of Ips typographus in southern Sweden. Agric For Meteorol. 2007:146:70–81. Joron M, Frezal L, Jones RT, Chamberlain NL, Lee SF, Haag CR, Whibley A, Becuwe M, Baxter SW, Ferguson L, et al. Chromosomal rearran- gements maintain a polymorphic supergene controlling butterfly mimicry. Nature. 2011:477(7363):203–206. https://doi.org/10. 1038/nature10341. Kandasamy D, Gershenzon J, Andersson MN, Hammerbacher A. Volatile organic compounds influence the interaction of the Eurasian spruce bark beetle (Ips typographus) with its fungal sym- bionts. ISME J. 2019:13(7):1788–1800. https://doi.org/10.1038/ s41396-019-0390-3. Kandasamy D, Zaman R, Nakamura Y, Zhao T, Hartmann H, Andersson MN, Hammerbacher A, Gershenzon J. Conifer-killing bark beetles locate fungal symbionts by detecting volatile fungal metabolites of host tree resin monoterpenes. PLoS Biol. 2023:21(2):e3001887. https://doi.org/10.1371/journal.pbio.3001887. Kapun M, Fabian DK, Goudet J, Flatt T. Genomic evidence for adaptive inversion clines in Drosophila melanogaster. Mol Biol Evol. 2016: 33(5):1317–1336. https://doi.org/10.1093/molbev/msw016. Kapun M, Flatt T. The adaptive significance of chromosomal inversion polymorphisms in Drosophila melanogaster. Mol Ecol. 2019:28(6): 1263–1282. https://doi.org/10.1111/mec.14871. Keightley PD, Pinharanda A, Ness RW, Simpson F, Dasmahapatra KK, Mallet J, Davey JW, Jiggins CD. Estimation of the spontaneous mu- tation rate in Heliconius melpomene. Mol Biol Evol. 2015:32(1): 239–243. https://doi.org/10.1093/molbev/msu302. Kim KW, De-Kayne R, Gordon IJ, Omufwoko KS, Martins DJ, Ffrench-Constant R, Martin SH. Stepwise evolution of a butterfly supergene via duplication and inversion. Philos Trans R Soc Lond B Biol Sci. 2022:377(1856):20210207. https://doi.org/10.1098/ rstb.2021.0207. Kirubakaran TG, Grove H, Kent MP, Sandve SR, Baranski M, Nome T, De Rosa MC, Righino B, Johansen T, Otterå H, et al. Two adjacent inversions maintain genomic differentiation between migratory and stationary ecotypes of Atlantic cod. Mol Ecol. 2016:25(10): 2130–2143. https://doi.org/10.1111/mec.13592. Koch EL, Morales HE, Larsson J, Westram AM, Faria R, Lemmon AR, Lemmon EM, Johannesson K, Butlin RK. Genetic variation for adaptive traits is associated with polymorphic inversions in Littorina saxatilis. Evol Letters. 2021:5(3):196–213. https://doi. org/10.1002/evl3.227. Korneliussen TS, Albrechtsen A, Nielsen R. ANGSD: analysis of next generation sequencing data. BMC Bioinformatics. 2014:15(1): 356. https://doi.org/10.1186/s12859-014-0356-4. Koštál V, Stetina T, Poupardin R, Korbelová J, Bruce AW. Conceptual framework of the eco-physiological phases of insect diapause de- velopment justified by transcriptomic profiling. Proc Natl Acad Sci U S A. 2017:114(32):8532–8537. https://doi.org/10.1073/pnas. 1707281114. Koury SA. Predicting recombination suppression outside chromosomal inversions in Drosophila melanogaster using crossover interference theory. Heredity (Edinb). 2023:130(4):196–208. https://doi.org/ 10.1038/s41437-023-00593-x. Kozak GM, Wadsworth CB, Kahne SC, Bogdanowicz SM, Harrison RG, Coates BS, Dopman EB. Genomic basis of circannual rhythm in the European corn borer moth. Curr Biol. 2019:29(20):3501–3509.e5. https://doi.org/10.1016/j.cub.2019.08.053. Krasovec M. The spontaneous mutation rate of Drosophila pseudoobs- cura. G3 (Bethesda). 2021:11(7):jkab151. https://doi.org/10.1093/ g3journal/jkab151. Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016:33(7):1870–1874. https://doi.org/10.1093/molbev/ msw054. Lamichhaney S, Fan G, Widemo F, Gunnarsson U, Thalmann DS, Hoeppner MP, Kerje S, Gustafson U, Shi C, Zhang H, et al. Structural genomic changes underlie alternative reproductive strategies in the ruff (Philomachus pugnax). Nat Genet. 2015:48(1):84–88. https://doi.org/10.1038/ng.3430. Lange H, Økland B, Krokene P. Thresholds in the life cycle of the spruce bark beetle under climate change. Inter J Complex Syst. 2006:1648:1–10. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012:9(4):357–359. https://doi.org/10.1038/ nmeth.1923. Lehmanski LMA, Kandasamy D, Andersson MN, Netherer S, Alves EG, Huang J, Hartmann H. Addressing a century-old hypothesis—do pioneer beetles of Ips typographus use volatile cues to find suitable host trees? New Phytol. 2023:238(5):1762–1770. https://doi.org/ 10.1111/nph.18865. Lewin TD, Liao IJ, Luo Y. Conservation of animal genome structure is the exception not the rule. bioRxiv 606322. https://doi.org/10. 1101/2024.08.02.606322, 6 August 2024, preprint: not peer reviewed. Li H. A statistical framework for SNP calling, mutation discovery, asso- ciation mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011:27(21):2987–2993. https://doi.org/10.1093/bioinformatics/btr509. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018:34(18):3094–3100. https://doi.org/10.1093/ bioinformatics/bty191. Li H, Berent E, Hadjipanteli S, Galey M, Muhammad-Lahbabi N, Miller DE, Crown KN. Heterozygous inversion breakpoints suppress mei- otic crossovers by altering recombination repair outcomes. PLoS Genet. 2023:19(4):e1010702. https://doi.org/10.1371/journal. pgen.1010702. Li H, Ralph P. Local PCA shows how the effect of population structure differs along the genome. Genetics. 2019:211(1):289–304. https://doi.org/10.1534/genetics.118.301747. Li J, Li H, Jakobsson M, Li S, Sjödin P, Lascoux M. Joint analysis of dem- ography and selection in population genetics: where do we stand and where could we go? Mol Ecol. 2012:21(1):28–44. https://doi. org/10.1111/j.1365-294X.2011.05308.x. Lischer HEL, Excoffier L. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics. 2012:28(2):298–299. https://doi.org/10.1093/ bioinformatics/btr642. Lohse K, Clarke M, Ritchie MG, Etges WJ. Genome-wide tests for intro- gression between cactophilic Drosophila implicate a role of inver- sions during speciation. Evolution. 2015:69(5):1178–1190. https://doi.org/10.1111/evo.12650. Lowry DB, Willis JH. A widespread chromosomal inversion polymorph- ism contributes to a major life-history transition, local adaptation, and reproductive isolation. PLoS Biol. 2010:8(9):e1000500. https:// doi.org/10.1371/journal.pbio.1000500. Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 21 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Lynch M, Conery JS. The origins of genome complexity. Science. 2003:302(5649):1401–1404. https://doi.org/10.1126/science. 1089370. Marroni F, Pinosio S, Morgante M. Structural variation and genome complexity: is dispensable really dispensable? Curr Opin Plant Biol. 2014:18:31–36. https://doi.org/10.1016/j.pbi.2014.01.003. Matschiner M, Barth JMI, Tørresen OK, Star B, Baalsrud HT, Brieuc MSO, Pampoulie C, Bradbury I, Jakobsen KS, Jentoft S. Supergene origin and maintenance in Atlantic cod. Nat Ecol Evol. 2022:6(4): 469–481. https://doi.org/10.1038/s41559-022-01661-x. Mayer F, Piel FB, Cassel-Lundhagen A, Kirichenko N, Grumiau L, Økland B, Bertheau C, Grégoire JC, Mardulyn P. Comparative mul- tilocus phylogeography of two Palaearctic spruce bark beetles: influence of contrasting ecological strategies on genetic variation. Mol Ecol. 2015:24(6):1292–1310. https://doi.org/10.1111/mec. 13104. McClung CE. The chromosome complex of orthopteran spermatocytes. Biol Bull. 1905:9(5):304–340. https://doi.org/10.2307/1535568. McCulloch GA, Waters JM. Rapid adaptation in a fast-changing world: emerging insights from insect genomics. Glob Chang Biol. 2023:29(4):943–954. https://doi.org/10.1111/gcb.16512. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next- generation DNA sequencing data. Genome Res. 2010:20(9): 1297–1303. https://doi.org/10.1101/gr.107524.110. McVean G, Awadalla P, Fearnhead P. A coalescent-based method for detecting and estimating recombination from gene sequences. Genetics. 2002:160(3):1231–1241. https://doi.org/10.1093/ genetics/160.3.1231. Mérot C, Berdan EL, Cayuela H, Djambazian H, Ferchaud AL, Laporte M, Normandeau E, Ragoussis J, Wellenreuther M, Bernatchez L. Locally adaptive inversions modulate genetic variation at different geographic scales in a seaweed fly. Mol Biol Evol. 2021:38(9): 3953–3971. https://doi.org/10.1093/molbev/msab143. Mérot C, Llaurens V, Normandeau E, Bernatchez L, Wellenreuther M. Balancing selection via life-history trade-offs maintains an inver- sion polymorphism in a seaweed fly. Nat Commun. 2020:11(1): 670. https://doi.org/10.1038/s41467-020-14479-7. Müller M, Olsson P-O, Eklundh L, Jamali S, Ardö J. Features predispos- ing forest to bark beetle outbreaks and their dynamics during drought. For Ecol Manage. 5232022:523:120480. Navarro A, Barbadilla A, Ruiz A. Effect of inversion polymorphism on the neutral nucleotide variability of linked chromosomal regions in Drosophila. Genetics. 2000:155(2):685–698. https://doi.org/ 10.1093/genetics/155.2.685. Nelson CW, Moncla LH, Hughes AL. SNPGenie: estimating evolution- ary parameters to detect natural selection using pooled next- generation sequencing data. Bioinformatics. 2015:31(22): 3709–3711. https://doi.org/10.1093/bioinformatics/btv449. Nilssen AC. Long-range aerial dispersal of bark beetles and bark wee- vils (Coleoptera, Scolytidae and Curculionidae) in northern Finland. Ann Entomol Fenn. 1984:50:37–42. Novo I, Ordás P, Moraga N, Santiago E, Quesada H, Caballero A. Impact of population structure in the estimation of recent historical effective population size by the software GONE. Genet Sel Evol. 2023:55(1):86. https://doi.org/10.1186/s12711-023-00859-2. Nunez JCB, Lenhart BA, Bangerter A, Murray CS, Mazzeo GR, Yu Y, Nystrom TL, Tern C, Erickson PA, Bergland AO. A cosmopolitan in- version facilitates seasonal adaptation in overwintering Drosophila. Genetics. 2023:226(2):iyad207. https://doi.org/10. 1093/genetics/iyad207. Öhrn P, Berlin M, Elfstrand M, Krokene P, Jönsson AM. Seasonal vari- ation in Norway spruce response to inoculation with bark beetle- associated bluestain fungi one year after a severe drought. For Ecol Manag. 2021:496:119443. Ohta T. Associative overdominance caused by linked detrimental mu- tations. Genet Res. 1971:18(3):277–286. https://doi.org/10.1017/ S0016672300012684. Oppold AM, Pfenninger M. Direct estimation of the spontaneous mu- tation rate by short-term mutation accumulation lines in Chironomus riparius. Evol Lett. 2017:1(2):86–92. https://doi.org/ 10.1002/evl3.8. Paolucci S, Salis L, Vermeulen CJ, Beukeboom LW, van de Zande L. QTL analysis of the photoperiodic response and clinal distribution of period alleles in Nasonia vitripennis. Mol Ecol. 2016:25(19): 4805–4817. https://doi.org/10.1111/mec.13802. Porubsky D, Höps W, Ashraf H, Hsieh PH, Rodriguez-Martin B, Yilmaz F, Ebler J, Hallast P, Maria Maggiolini FA, Harvey WT, et al. Recurrent inversion polymorphisms in humans associate with gen- etic instability and genomic disorders. Cell. 2022:185(11): 1986–2005.e26. https://doi.org/10.1016/j.cell.2022.04.017. Powell D, Groβe-Wilde E, Krokene P, Roy A, Chakraborty A, Löfstedt C, Vogel H, Andersson MN, Schlyter F. A highly-contiguous genome assembly of the Eurasian spruce bark beetle, Ips typographus, pro- vides insight into a major forest pest. Commun Biol. 2021:4(1): 1059. https://doi.org/10.1038/s42003-021-02602-3. Pruisscher P, Nylin S, Gotthard K, Wheat CW. Genetic variation under- lying local adaptation of diapause induction along a cline in a butterfly. Mol Ecol. 2018:27(18):3613–3626. https://doi.org/10. 1111/mec.14829. Purcell J, Brelsford A, Wurm Y, Perrin N, Chapuisat M. Convergent genetic architecture underlies social organization in ants. Curr Biol. 2014:24(22):2728–2732. https://doi.org/10.1016/j.cub. 2014.09.071. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, de Bakker PIW, Daly MJ, et al. PLINK: a tool set for whole-genome association and population-based linkage ana- lyses. Am J Hum Genet. 2007:81(3):559–575. https://doi.org/10. 1086/519795. Raffa KF, Andersson MN, Schlyter F. Host selection by bark beetles: playing the odds in a high-stakes game. In: Tittiger C, Blomquist GJ, editors. Advances in insect physiology. Cambridge, MA: Academic Press; 2016. p. 1–74. Rane RV, Rako L, Kapun M, Lee SF, Hoffmann AA. Genomic evidence for role of inversion 3RP of Drosophila melanogaster in facilitating climate change adaptation. Mol Ecol. 2015:24(10):2423–2432. https://doi.org/10.1111/mec.13161. Raymond M, Rousset F. An exact test for population differentiation. Evolution. 1995:49(6):1280–1283. . https://doi.org/10.1111/j. 1558-5646.1995.tb04456.x. Reeve J, Butlin RK, Koch EL, Stankowski S, Faria R. Chromosomal inver- sion polymorphisms are widespread across the species ranges of rough periwinkles (Littorina saxatilis and L. arcana). Mol Ecol. 2023:33:e17160. https://doi.org/10.1111/mec.17160. Roesti M, Gilbert KJ, Samuk K. Chromosomal inversions can limit adap- tation to new environments. Mol Ecol. 2022:31(17):4435–4439. https://doi.org/10.1111/mec.16609. Roff DA. The evolution of threshold traits in animals. Q Rev Biol. 1996:71(1):3–35. https://doi.org/10.1086/419266. Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guirao-Rico S, LibradoP, Ramos-Onsins SE, Sanchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017: 34(12):3299–3302. https://doi.org/10.1093/molbev/msx248. Saitou M, Masuda N, Gokcumen O. Similarity-based analysis of allele frequency distribution among multiple populations identifies adaptive genomic structural variants. Mol Biol Evol. 2022:39(3): msab313. https://doi.org/10.1093/molbev/msab313. Mykhailenko et al. GBE 22 Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024 Sallé A, Arthofer W, Lieutier F, Stauffer C, Kerdelhué C. Phylogeography of a host-specific insect: genetic structure of Ips typographus in Europe does not reflect past fragmentation of its host. Biol J Linn Soc. 2007:90(2):239–246. https://doi.org/10. 1111/j.1095-8312.2007.00720.x. Schaal SM, Haller BC, Lotterhos KE. Inversion invasions: when the gen- etic basis of local adaptation is concentrated within inversions in the face of gene flow. Philos Trans R Soc Lond B Biol Sci. 2022: 377(1856):20210200. https://doi.org/10.1098/rstb.2021.0200. Schebeck M, Dobart N, Ragland GJ, Schopf A, Stauffer C. Facultative and obligate diapause phenotypes in populations of the European spruce bark beetle Ips typographus. J Pest Sci. 2022:95(2):889–899. https://doi.org/10.1007/s10340-021- 01416-w. Schebeck M, Hansen EM, Schopf A, Ragland GJ, Stauffer C, Bentz BJ. Diapause and overwintering of two spruce bark beetle species. Physiol Entomol. 2017:42(3):200–210. https://doi.org/10.1111/ phen.12200. Schebeck M, Lehmann P, Laparie M, Bentz BJ, Ragland GJ, Battisti A, Hahn DA. Seasonality of forest insects: why diapause matters. Trends Ecol Evol. 2024:39(8):757–770. https://doi.org/10.1016/j. tree.2024.04.010. Schopf A. Zum einfluß der photoperiode auf die entwicklung und Kälteresistenz des Buchdruckers, Ips typographus L. (Col., Scolytidae). J Pest Sci. 1985:58:73–75. https://doi.org/10.1007/ BF01903228. Schopf A. Die wirkung der photoperiode auf die induktion der imagi- naldiapause von Ips typographus (L.) (Col., Scolytidae). J Appl Entomol. 1989:107(1-5):275–288. https://doi.org/10.1111/j. 1439-0418.1989.tb00257.x. Schroeder M, Dalin P. Differences in photoperiod-induced diapause plas- ticity among different populations of the bark beetle Ips typographus and its predator Thanasimus formicarius. Agric For Entomol. 2017:19(2):146–153. https://doi.org/10.1111/afe.12189. Schwander T, Libbrecht R, Keller L. Supergenes and complex pheno- types. Curr Biol. 2014:24(7):R288–R294. https://doi.org/10.1016/ j.cub.2014.01.056. Shen W, Sipos B, Zhao L. SeqKit2: a Swiss army knife for sequence and alignment processing. Imeta. 2024:3(3):e191. https://doi.org/10. 1002/imt2.191. Sirén J, Monlong J, Chang X, Novak AM, Eizenga JM, Markello C, Sibbesen JA, Hickey G, Chang PC, Carroll A, et al. Pangenomics en- ables genotyping of known structural variants in 5202 diverse gen- omes. Science. 2021:374(6574):abg8871. https://doi.org/10. 1126/science.abg8871. Skotte L, Korneliussen TS, Albrechtsen A. Estimating individual admix- ture proportions from next generation sequencing data. Genetics. 2013:195(3):693–702. https://doi.org/10.1534/genetics.113. 154138. Stauffer C, Lakatos F, Hewitt GM. Phylogeography and postglacial col- onization routes of Ips typographus L. (Coleoptera, Scolytidae). Mol Ecol. 1999:8(5):763–773. https://doi.org/10.1046/j.1365- 294X.1999.00626.x. Sturtevant AH. Linkage variation and chromosome maps. Proc Natl Acad Sci U S A. 1921:7(7):181–183. https://doi.org/10.1073/ pnas.7.7.181. Thompson MJ, Jiggins CD. Supergenes and their role in evolution. Heredity (Edinb). 2014:113(1):1–8. https://doi.org/10.1038/hdy. 2014.20. Tigano A, Friesen VL. Genomics of local adaptation with gene flow. Mol Ecol. 2016:25(10):2144–2164. https://doi.org/10.1111/mec. 13606. Tigano A, Jacobs A, Wilder AP, Nand A, Zhan Y, Dekker J, Therkildsen NO. Chromosome-level assembly of the Atlantic silverside genome reveals extreme levels of sequence diversity and structural genetic variation. Genome Biol Evol. 2021:13(6):evab098. https://doi.org/ 10.1093/gbe/evab098. Todesco M, Owens GL, Bercovich N, Légaré JS, Soudi S, Burge DO, Huang K, Ostevik KL, Drummond EBM, Imerovski I, et al. Massive haplotypes underlie ecotypic differentiation in sunflowers. Nature. 2020:584(7822):602–607. https://doi.org/10.1038/ s41586-020-2467-6. Twyford AD, Friedman J. Adaptive divergence in the monkey flower Mimulus guttatus is maintained by a chromosomal inversion. Evolution. 2015:69(6):1476–1486. https://doi.org/10.1111/evo. 12663. Van’t Hof AE, Campagne P, Rigden DJ, Yung CJ, Lingley J, Quail MA, Hall N, Darby AC, Saccheri IJ. The industrial melanism mutation in British peppered moths is a transposable element. Nature. 2016:534(7605):102–105. https://doi.org/10.1038/nature17951. Vega FE, Hofstetter RW. Bark beetles: biology and ecology of native and invasive species. Cambridge, MA: Academic Press; 2014. Wang J, Wurm Y, Nipitwattanaphon M, Riba-Grognuz O, Huang YC, Shoemaker D, Keller L. A Y-like social chromosome causes alterna- tive colony organization in fire ants. Nature. 2013:493(7434): 664–668. https://doi.org/10.1038/nature11832. Wang Z, Liu Y, Wang H, Roy A, Liu H, Han F, Zhang X, Lu Q. Genome and transcriptome of Ips nitidus provide insights into high-altitude hypoxia adaptation and symbiosis. iScience. 2023:26(10):107793. https://doi.org/10.1016/j.isci.2023.107793. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984:38(6):1358–1370. https:// doi.org/10.1111/j.1558-5646.1984.tb05657.x. Weissensteiner MH, Bunikis I, Catalán A, Francoijs KJ, Knief U, Heim W, Peona V, Pophaly SD, Sedlazeck FJ, Suh A, et al. Discovery and population genomics of structural variation in a songbird genus. Nat Commun. 2020:11(1):3403. https://doi.org/10.1038/s41467- 020-17195-4. Wellenreuther M, Bernatchez L. Eco-evolutionary genomics of chromosomal inversions. Trends Ecol Evol. 2018:33(6):427–440. https://doi.org/10.1016/j.tree.2018.04.002. Wellenreuther M, Mérot C, Berdan E, Bernatchez L. Going beyond SNPs: the role of structural genomic variants in adaptive evolution and species diversification. Mol Ecol. 2019:28(6):1203–1209. https://doi.org/10.1111/mec.15066. Wold JR, Guhlin JG, Dearden PK, Santure AW, Steeves TE. The promise and challenges of characterizing genome-wide structural variants: a case study in a critically endangered parrot. Mol Ecol Resour. 2023: Early View, https://doi.org/10.1111/1755-0998.13783. Yang YY, Lin FJ, Chang HY. Comparison of recessive lethal accumula- tion in inversion-bearing and inversion-free chromosomes in Drosophila. Zool Stud. 2002:41:271–282. Yuvaraj JK, Roberts RE, Sonntag Y, Hou XQ, Grosse-Wilde E, Machara A, Zhang DD, Hansson BS, Johanson U, Löfstedt C, et al. Putative ligand binding sites of two functionally characterized bark beetle odorant receptors. BMC Biol. 2021:19(1):16. https://doi.org/10. 1186/s12915-020-00946-6. Zhao T, Kandasamy D, Krokene P, Chen J, Gershenzon J, Hammerbacher A. Fungal associates of the tree-killing bark beetle, Ips typographus, vary in virulence, ability to degrade conifer phe- nolics and influence bark beetle tunneling behavior. Fungal Ecol. 2019:38:71–79. https://doi.org/10.1016/j.funeco.2018.06.003. Associate editor: Tatiana Giraud Complex Genomic Landscape of Inversion Polymorphism GBE Genome Biol. Evol. 16(12) https://doi.org/10.1093/gbe/evae263 Advance Access publication 4 December 2024 23 D ow nloaded from https://academ ic.oup.com /gbe/article/16/12/evae263/7916417 by N atural R esources Institute Finland (Luke) user on 31 D ecem ber 2024