# Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species

Eric Lombaert, Aurélie Blin, Barbara Porro, Thomas Guillemaud, Julio S Bernal, Gary Chang, Natalia Kirichenko, Thomas W Sappington, Stefan Toepfer, Emeline Deleury

► To cite this version:

Eric Lombaert, Aurélie Blin, Barbara Porro, Thomas Guillemaud, Julio S Bernal, et al.. Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species. Peer Community Journal, 2025, 5, pp.e36. 10.24072/pcjournal.539. hal-05246288

HAL Id: hal-05246288

https://hal.science/hal-05246288v1

Submitted on 9 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Research article

Published
2025-04-02

**Cite as**
Eric Lombaert, Aurélie Blin, Barbara Porro, Thomas Guillemaud, Julio S. Bernal, Gary Chang, Natalia Kirichenko, Thomas W. Sappington, Stefan Toepfer and Emeline Deleury (2025) *Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species*, Peer Community Journal, 5: e36.

**Correspondence**
eric.lombaert@inrae.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Evolutionary Biology,
https://doi.org/10.24072/pci.
evolbiol.100845

<figure><img src="image_2.png" /></figure>

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species

Eric Lombaert<sup>1,1</sup>, Aurélie Blin<sup>1</sup>, Barbara Porro<sup>1,2</sup>, Thomas Guillemaud<sup>1,1</sup>, Julio S. Bernal<sup>1,3</sup>, Gary Chang<sup>1,4</sup>, Natalia Kirichenko<sup>1,5,6,7</sup>, Thomas W. Sappington<sup>1,8,9</sup>, Stefan Toepfer<sup>10,11</sup>, and Emeline Deleury<sup>1,1</sup>

Volume 5 (2025), article e36

https://doi.org/10.24072/pcjournal.539

## Abstract

Many invasive species undergo a significant reduction in genetic diversity, i.e. a genetic bottleneck, in the early stages of invasion. However, this reduction does not necessarily prevent them from achieving considerable ecological success and becoming highly efficient colonizers. Here we investigated the purge hypothesis, which suggests that demographic bottlenecks may facilitate conditions (e.g., increased homozygosity and inbreeding) under which natural selection can purge deleterious mutations, thereby reducing genetic load. We used a transcriptome-based exome capture protocol to identify thousands of SNPs in coding regions of native and invasive populations of two highly successful invasive insect species, the western corn rootworm (Chrysomelidae: *Diabrotica virgifera* virgifera) and the harlequin ladybird (Coccinellidae: *Harmonia axyridis*). We categorized and polarized SNPs to investigate changes in genetic load between invasive populations and their sources. Our results differed between species. In *D. virgifera virgifera*, although there was a general reduction in genetic diversity in invasive populations, including that associated with genetic load, we found no clear evidence for purging of genetic load, except marginally for highly deleterious mutations in one European population. Conversely, in *H. axyridis*, the reduction in genetic diversity was minimal, and we detected signs of genetic load fixation in invasive populations. These findings provide new insights into the evolution of genetic load during invasions, but do not offer a definitive answer to the purge hypothesis. Future research should include larger genomic datasets and a broader range of invasive species to further elucidate these dynamics.

<footnote n="1">¹INRAE, Université Côte d'Azur, ISA, 400 Route des Chappes, 06903 Sophia-Antipolis, France, ²Génomique Métabolique, Genoscope, Institut François Jacob, CEA, CNRS, Univ Evry, 2 Rue Gaston Crémieux, 91000 Évry-Courcouronnes, France, ³Department of Entomology, Texas A&M University, 370 Olsen Blvd., College Station, TX 77843-2475, USA, ⁴Biology Department, Gonzaga University, 502 East Boone Ave., Spokane, WA 99258-0102, USA, ⁵Sukachev Institute of Forest, Siberian Branch of the Russian Academy of Sciences, Federal Research Center "Krasnoyarsk Science Center SB RAS", Akademgorodok 50/28, 660036 Krasnoyarsk, Russia, ⁶Institute of Ecology and Geography, Siberian Federal University, Svobodny pr. 79, 660041 Krasnoyarsk, Russia, ⁷All-Russian Plant Quarantine Center, Krasnoyarsk branch, Zhelyabova ut. 6/6, 660020 Krasnoyarsk, Russia, ⁸USDA-Agricultural Research Service, Corn Insects & Crop Genetics Research Unit, 1012 Crop Genome Informatics Lab, 819 Wallace Road, Ames, IA 50011-4014, USA, ⁹Iowa State University, Department of Plant Pathology, Entomology, and Microbiology, 1344 Advanced Teaching and Research Building, 2213 Pammel Drive, Ames, IA 50011, USA, ¹⁰CABI Switzerland, c/o Plant Protection and Soil Conservation Directorate, Rarosi ut. 110, H 6800 Hodmezovasarhely, Hungary, ¹¹MARA-CABI Joint Laboratory for Biosafety, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China</footnote>
Introduction

Biological invasions represent a significant aspect of global change, profoundly impacting
biodiversity through the alteration of species distributions worldwide, particularly in recent times
due to the significant increase in human-assisted dispersal (Seebens et al., 2017). The key factors
determining the success of invasive species are not yet fully understood, although numerous
hypotheses have been proposed (Enders et al., 2018; Sherpa & Després, 2021; Daly et al., 2023).
Most proposals fail to account for one or more consistent characteristics of successful invasions:
(i) the rarity of successful invasions resulting from introductions (Williamson & Fitter, 1996), (ii) the
lag time between initial introduction and invasion (Sakai et al., 2001), (iii) the frequent reduction in
genetic diversity due to demographic bottlenecks (Nei et al., 1975), and (iv) the prevalence of
multiple invasions originating from an initial invasive population (i.e. the bridgehead effect;
Lombaert et al., 2010). These characteristics suggest that the success of invasions may stem from
partially stochastic biological processes spanning multiple generations and combining both
demographic and genetic mechanisms.

A hypothesis that aligns with these characteristics is that genetic load is purged in the initial stages of invasions, i.e. deleterious alleles that cause inbreeding depression are eliminated from the introduced population (Estoup et al., 2016). During the demographic bottleneck following introduction of a few individuals into a new environment, genetic drift intensifies and homozygosity increases. While genetic drift may randomly eliminate some deleterious alleles, thereby reducing part of the genetic load, it also contributes to the transformation of the masked load (i.e., the load which may become express in future generations; Bertorelle et al., 2022) into a realized load (i.e., the load which reduces fitness in the current generation; Bertorelle et al., 2022). This shift may result in a "mutational meltdown" (Lynch et al., 1995; Simberloff, 2009), where expression of deleterious mutations increases risk of extinction, ultimately leading to failure of the introduced population to establish. Conversely, the increased homozygosity resulting from bottlenecks exposes recessive deleterious alleles to selection. In the short term, this exposure can reduce mean fitness, but if the population persists, purifying selection may progressively remove some highly deleterious alleles over multiple generations, potentially reducing genetic load (Dussex et al., 2023). Theoretically, such purges can occur under specific demographic (i.e., moderate reduction of population size) and genetic (i.e., highly deleterious, recessive alleles) conditions, optimizing exposure to natural selection (Crow, 1970; Charlesworth et al., 1990; Glémin, 2003; Robinson et al., 2023). The potential for such purging in introduced populations is particularly interesting for its role in facilitating successful invasions by reducing inbreeding depression and enabling inbred individuals to maintain high fitness levels.

Purging has been demonstrated empirically by measuring the evolution of various traits in
artificially bottlenecked populations (Crnokrak & Barrett, 2002; Ávila et al., 2010). However,
measuring life history traits can be challenging in the context of naturally occurring bottlenecks
during biological invasions, and evidence of purging has been documented in only a few invasive
species (Parisod et al., 2005; Mullarkey et al., 2013; Fountain et al., 2014; Marchini et al., 2016).
One notable formal test of this hypothesis was conducted on the invasive Asian ladybird *Harmonia
axyridis*, where measurement of life history traits revealed that invasive populations showed no
evidence of the inbreeding depression observed in native ones, suggesting that deleterious alleles
were purged during the invasion process (Facon et al., 2011). Overall, case studies focusing on
the dynamics of genetic load during biological invasions have either examined a few life history
traits (e.g. Facon et al., 2011) or used a single locus approach (e.g. Zayed et al., 2007), making it
difficult to generalize the results.

Advances in population genomics over the past decade have provided promising avenues for
investigating genetic load on large scales. Initially explored in humans (Lohmueller et al., 2008;
Henn et al., 2015), these approaches have been applied in the fields of domestication (e.g.
Schubert et al., 2014; Marsden et al., 2016; Makino et al., 2018; Wang et al., 2021) and
conservation (Xue et al., 2015; Robinson et al., 2016; Grossen et al., 2020; Dussex et al., 2021;
Ochoa & Gibbs, 2021). These studies typically involve SNP-calling within coding regions, determination of ancestral SNP states, categorization of the severity of fitness reduction caused by derived alleles, and population comparisons while accounting for genetic drift using synonymous and/or intergenic polymorphisms. High quality genomic resources are essential for such studies, which may explain the limited application of these methods to invasive species, which are mostly non-model organisms. However, advances in genome sequencing technologies and bioinformatics have significantly reduced costs and improved accessibility, making these methods increasingly routine and affordable (Bertorelle et al., 2022).

In this study, we examined the dynamics of genetic load during the invasion of two successful insect species, namely the western corn rootworm, *Diabrotica virgifera virgifera*, and the harlequin ladybird, *Harmonia axyridis*, by directly measuring/estimating genetic load using genomic data from feral native and invasive populations. Importantly, our study was not designed to assess the instrumental role of purging in invasion success, but rather to investigate its occurrence during these invasions. A broader study of both successful and failed invasions across many species would be required to draw conclusions on purging's role in invasion success. We used a pool-seq transcriptome-based exome capture protocol previously developed for non-model species (Deleury et al., 2020) to identify SNPs within coding sequences and categorize them as synonymous, moderately deleterious, or highly deleterious. Our results offer insights into the fate of the genetic load and provide valuable perspectives on the purge hypothesis in the context of biological invasions.

## Methods

### Species choice

The two species, *Diabrotica virgifera virgifera* (hereafter DVV) and *Harmonia axyridis* (HA), are good candidates for testing the purging of genetic load hypothesis. Both species are highly successful invaders with extensive invasive ranges (Gray et al., 2009; Roy et al., 2016). Additionally, their invasion routes are well-documented and supported by robust analyses using diverse datasets and methodologies (Miller et al., 2005; Ciosi et al., 2008; Lombaert et al., 2010, 2011, 2014, 2018). Both species exhibit bridgehead effects, meaning that certain invasive populations have played a crucial role in global dissemination (Guillemaud et al., 2011). Additionally, purging of genetic load was detected in HA in a laboratory study focusing on life-history traits (Facon et al., 2011).

### Design of exome capture probes for target enrichment

Design of the exome capture probes was performed as described in Deleury et al. (2020). We searched for peptide-coding sequences using FRAMEDP (v1.2.2; Gouzy et al., 2009) on the *de novo* transcriptomes described in Coates et al. (2021) and Vogel et al. (2017) for DVV and HA, respectively. BLASTX results (e-value ≤ 1e⁻⁷) of transcripts were used against the insect proteomes of *Tribolium castaneum*, *Anoplophora glabripennis*, *Dendroctonus ponderosae*, *Drosophila melanogaster* and the SwissProt database (v2016-02-17) for training. From the obtained coding sequences (CDS), we eliminated (i) *Wolbachia* and other putative endosymbiont sequences, (ii) CDS with > 1% missing nucleotide bases (Ns) or with more than four consecutive Ns, (iii) CDS with a GC% below 25 or above 75 and (iv) CDS with lengths < 400 bp or > 3500 bp. From the remaining 7,132 and 12,739 CDS for DVV and HA respectively, we drew at random c.a. 5.5 Mb for each species.

Probes based on the selected CDS were designed and manufactured by NimbleGen. In the case of DVV, repetitiveness of the probes was assessed based on the highest 15-mer frequency among the genomes of *Tribolium castaneum* (GCA_000002335.3), *Dendroctonus ponderosae* (GCA_000355655.1 and GCA_000346045.2) and DVV (the one available at the time, GCA_003013835.2). Probes with more than five close matches in the DVV genome were discarded. Close matches were defined as no more than five single-base insertions, deletions or
substitutions using the SSAHA algorithm (Ning et al., 2001). Probes that matched sequences in the mitochondrial genome were also discarded. Random nucleotides were used to replace residual Ns in target sequences. For HA, the method used to ensure probe uniqueness is described in Deleury *et al.* (2020).

Overall, the final probe sets corresponded to a total of 4,151 CDS (5,282,603 bases across 12,017 regions of overlapping probes) and 5,717 CDS (5,347,461 bases across 6,400 regions of overlapping probes; Deleury et al., 2020) for DVV and HA respectively. This final set of probes was manufactured in the form of biotinylated DNA oligomers. One capture reaction contained 2,100,000 overlapping probes of 50 to 99 bp in length (mean length of 73.86 ± 4.46 bp for DVV, and 74.71 ± 4.92 bp for HA).

## Sample collection

For both species, the choice of populations to be sampled was based on previously known invasion routes (Figure S1; Miller et al., 2005; Lombaert et al., 2010, 2018), so that each invasive population could be compared to its source. We sampled adult DVV at 6 sites: two in the native area (Mexico), two in North America (Colorado and Pennsylvania, respectively corresponding to the core and front of the first invasive population; Lombaert et al., 2018), and two in Europe (Hungary, referred to the Central-Southeastern European population, and north-western Italy, which are independently derived from the same source area of eastern North America; Miller et al., 2005). Adult HA were sampled at four sites: two in the native area (Russia (Siberia) and China), and two in North America (Pennsylvania and Washington, corresponding to two independent outbreaks from the native area; Lombaert et al., 2010, 2014). None of the selected invasive populations resulted from multiple introductions. The east North American population of HA was previously hypothesized to be an admixture between two populations (Lombaert et al., 2011), but our analysis using ABC (approximate Bayesian computation; Beaumont et al., 2002) with synonymous SNPs from the current dataset indicates that this is unlikely (see Appendix S1 for details).

To determine the ancestral and derived alleles for each SNP (i.e. to polarize the alleles), outgroup species were also sampled. For DVV, we selected the closely related species *Diabrotica adelpha*, as well as the more phylogenetically distant chrysomelid *Cerotoma trifurcata* (Eben & de los Monteros, 2013). In the case of HA, from phylogenetically closest to farthest (Tomaszewska et al., 2021), the sampled outgroup species were *Harmonia yedoensis*, *Harmonia conformis* and *Harmonia quadripunctata*. Complete information about samples is provided in Table 1 and Table S1.

## DNA extraction, exome capture and pool sequencing

For each population, 4 legs from each of 40 individuals (DVV) or 2 legs from each of 36 individuals (HA) were pooled for DNA extraction with the Qiagen DNeasy Blood & Tissue kit, in accordance with the manufacturer's recommendations. For the outgroup species, the same kit was used to extract DNA from a single individual each for the three *Harmonia* species, and from a pool of two individuals each for the *Diabrotica* and *Cerotoma* species.

Genomic libraries were prepared using NimbleGen SeqCap EZ HyperCap Library v2.0 and NimbleGen SeqCap EZ Library v5.0 for DVV and HA, respectively. In brief, for each of the population and outgroup samples, DNA (2 µg in 100 µl) was mechanically sheared to an average size of 200 bp using a Covaris S2 E210 device (6 cycles of 30 seconds each). In the case of DVV, the fragmented DNA was then divided into three technical replicates for each population. Subsequently, the fragments were subjected to end-repair, A-tailing and indexing (with one unique index per sample) using the KAPA Library Preparation kit designed for Illumina platforms. Following the ligation of Illumina adapters and indexes, only fragments falling within the size range of 250 to 450 bp were retained. A PCR amplification step consisting of 7 cycles was performed using standard Illumina paired-end primers, and the resulting amplicons were purified using AMPure XP beads (Beckman). The length, quality, and concentration of the prepared DNA
fragments were assessed using a BioAnalyzer with Agilent High Sensitivity DNA Assay, along with a Qubit.

<table><caption>**Table 1** - Description of *Diabrotica virgifera virgifera* (DVV) and *Harmonia axyridis* (HA) samples, including outgroup species (see Table S1 for BioSample accessions). Sources of invasive populations were determined based on known invasion routes (Figure S1; Miller et al., 2005; Lombaert et al., 2010, 2018), except for H-I-PEN for which we retraced invasion routes in this study (see Appendix S1 for details). Because the two native DVV populations are genetically almost indistinguishable, we selected D-N-MX1 as the source of D-I-COL based on the smallest pairwise $F_{ST}$ value between native and invasive populations (see Results section).</caption><thead><tr><th>Species</th><th>Population code name</th><th>Status (source)</th><th>Sampling site (Lat.; long.)</th><th>Sampling date</th><th>Haploid sample size</th></tr></thead><tbody><tr><td rowspan="7">Diabrotica virgifera virgifera</td><td>D-N-MX1</td><td>Native</td><td>Canatlán, Mexico (24.554; -104.741)</td><td>Oct. 2015</td><td>80</td></tr><tr><td>D-N-MX2</td><td>Native</td><td>Durango, Mexico (23.999; -104.614)</td><td>Oct. 2015</td><td>80</td></tr><tr><td>D-I-COL</td><td>Invasive (D-N-MX1)</td><td>Fort Morgan, CO, USA (40.218; -103.869)</td><td>Aug. 2015</td><td>80</td></tr><tr><td>D-I-PEN</td><td>Invasive (D-I-COL)</td><td>Landisville, PA, USA (40.119; -76.432)</td><td>Aug. 2015</td><td>80</td></tr><tr><td>D-I-HUN</td><td>Invasive (D-I-PEN)</td><td>Kondoros, Hungary (46.736; 20.816)</td><td>July 2015</td><td>80</td></tr><tr><td>D-I-ITA</td><td>Invasive (D-I-PEN)</td><td>Cuneo, Italy (44.463; 7.571)</td><td>July 2015</td><td>80</td></tr><tr><td>H-N-CHI</td><td>Native</td><td>Beijing, China (40.057; 116.540)</td><td>Oct. 2015</td><td>72</td></tr><tr><td rowspan="4">Harmonia axyridis</td><td>H-N-RUS</td><td>Native</td><td>Krasnoyarsk, Russia (55.992; 92.757)</td><td>Sept. 2015</td><td>72</td></tr><tr><td>H-I-PEN</td><td>Invasive (H-N-CHI)</td><td>Landisville, PA, USA (40.119; -76.432)</td><td>Aug. 2015</td><td>72</td></tr><tr><td>H-I-WAS</td><td>Invasive (H-N-CHI)</td><td>Spokane, WA, USA (47.665; -117.403)</td><td>Nov. 2015</td><td>72</td></tr><tr><td>D-OG-ADE</td><td>Outgroup DVV</td><td>Oaxaca State, Mexico (15.926; 97.151)</td><td>Jan. 2017</td><td>4</td></tr><tr><td>Cerotoma trifurcata</td><td>D-OG-TRI</td><td>Outgroup DVV</td><td>Weldon, CA, USA (35.660; -118.330)</td><td>Oct. 2016</td><td>4</td></tr><tr><td>Harmonia yedoensis</td><td>H-OG-YED</td><td>Outgroup HA</td><td>Beijing, China (40.031; 116.265)</td><td>Feb. 2015</td><td>2</td></tr><tr><td>Harmonia conformis</td><td>H-OG-CON</td><td>Outgroup HA</td><td>Mouans-Sartoux, France (43.620; 6.934)</td><td>Oct. 2009</td><td>2</td></tr><tr><td>Harmonia quadripunctata</td><td>H-OG-QUA</td><td>Outgroup HA</td><td>Rearing, Valbonne, France (NA)</td><td>July 2007</td><td>2</td></tr></tbody></table>

For each capture (five for DVV and two for HA), we used a total of approximately 1 µg of amplified DNA for exome enrichment, combining multiple samples in proportions that allowed for equimolarity between population samples (see Table S2 and Table S3). This enrichment was performed using the capture probes described above, following the guidelines of either the NimbleGen SeqCap EZ HyperCap Library Protocol v2.0 or the SeqCap EZ Library Protocol v5.0. Following each capture, we conducted two parallel PCRs, each comprising 14 cycles, on the elution solution. In the case of HA, the resulting PCR products were combined. Subsequently, all PCR products underwent purification using AMPure XP beads. The length, quality, and concentration of the final DNA fragments were assessed using a BioAnalyzer equipped with the Agilent High Sensitivity DNA Assay and a Qubit fluorometer.

For sequencing, we used one lane of an Illumina HiSeq3000 sequencer per species, following the manufacturer's instructions, in paired-end mode for 150 cycles. After sequencing, the data
were demultiplexed and exported as FastQ files, and the libraries were processed independently.
For DVV, the FastQ files of technical replicates were merged prior to subsequent analysis.

Mapping, SNP calling and annotation

Sequence quality assessment was conducted using FastQC v0.11.5 (Andrews, 2010).
Subsequently, adapter sequences were removed, and low-quality base pairs were eliminated
using *Trimmomatic* v0.35 (Bolger et al., 2014), with the following parameter settings:
ILLUMINACLIP:TruSeq-file.fa:2:30:10 LEADING:25 TRAILING:25 SLIDINGWINDOW:5:30
MINLEN:75. We used almost the same parameters for the outgroup species sequences, but with
less stringency on quality of reads (SLIDINGWINDOW:5:20).

Filtered reads were mapped onto the genome assemblies PGI_DIABVI_V3a (GenBank identifier: GCA_917563875.2) and icHarAxyr1.1 (GenBank identifier: GCA_914767665.1), for DVV and HA respectively, using default options of the *bwa-mem* aligner v0.7.15 (Li, 2013). We used the *SAMtools* v1.15.1 software package (Li et al., 2009) and its *fixmate*, *view*, and *markdup* tools to perform the following operations sequentially: (1) removal of unmapped reads and secondary alignments, (2) filtering out read alignments with a mapping quality Phred-score <20 and improperly paired reads, and (3) identification and removal of PCR duplicates. Processing filtered reads for the outgroup species followed the same steps, including mapping to the genome of the corresponding focal species, except that PCR duplicates were not removed.

For each species, we conducted variant calling on the resulting *bam* alignment files using
Freebayes v1.3.6 software (Garrison & Marth, 2012), with a focus on exonic regions (options -t
exons-bed --pooled-continuous --min-alternate-count 1 --min-alternate-fraction 0.001 --min-
alternate-total 1 --min-coverage 80 --use-best-n-alleles 3). The resulting vcf file was subsequently
filtered using the *view* tool within bcftools v1.13 software (Danecek et al., 2021) and an in-house
script to retain only true bi-allelic SNPs. Annotation was then performed with the *SnpEff* program
v5.0 (Cingolani et al., 2012).

Allele polarization

SNP positions were extracted from both vcf files using the bcftools query tool. Subsequently, for each outgroup species, we identified the nucleotides (or absence of typing) at positions matching those in their respective focal species using the *samtools mpileup* tool on the previously generated *bam* alignment files.

To polarize the SNPs, we used est-sfs software v2.03 (Keightley & Jackson, 2018) with the
Kimura 2-parameter model. The input files contained data from one native population of the focal
species, together with the corresponding outgroup species, ensuring a consistent phylogenetic
topology for the software. Given that est-sfs computes the probability of the most frequent allele
being ancestral solely for polymorphic SNPs within the native population, we used probabilities
from a parallel analysis for monomorphic loci, incorporating an extra haplotype for each allele. The
complete process was repeated twice for each focal species, considering the availability of two
native populations for each. This resulted in two est-sfs probabilities per focal species, both
reflecting the likelihood of the most frequent allele being the ancestral one.

To consider a SNP as polarized, we applied the following rules. If both probabilities were above
0.5 for the same allele, and at least one of the probabilities exceeded 0.75, we considered the
most frequent allele as the ancestral one. If both probabilities were below 0.5 for the same allele,
with at least one below 0.25, and if there were no more than two distinct nucleotides present in
total across the focal and outgroup species combined, we considered the less frequent allele to be
the ancestral one. SNPs that did not meet these criteria were considered as not polarized but were
retained for computations that did not require polarization. Note that, because the est-sfs
probabilities were predominantly close to 1 (Figure S2), the threshold value had minimal impact on
the number of polarized SNPs (Figure S3).
## Additional filters and categorization of deleterious mutations

Key information was extracted from the vcf file using the *SnpSift* program v5.0 (Cingolani et al., 2012), including SNP coordinates, reference and alternative alleles, allele depths, predicted effects, impact annotations, and protein-level changes. Additional filtering was performed with an in-house R script (R Core Team, 2021). First, we retained only biallelic SNPs with coverage greater than 50 reads, falling below the 95th percentile of overall coverage in each pool, and exhibiting a minor allele frequency exceeding 0.01 in at least one population. Second, we retained SNPs located in coding regions with unambiguous annotations. Finally, we excluded SNPs located on the X-chromosome in the specific case of HA, for which this information was available.

*SnpEff* annotations were used to categorize the severity of fitness loss due to mutations. Synonymous variants were used as proxies for neutral polymorphism. Non-synonymous mutations were categorized as either “missense” (considered as potentially moderately deleterious, involving amino acid changes) or “LoF” (loss-of-function, considered as potentially highly deleterious, involving gain or loss of stop codons).

## Genetic diversity and genetic load analyses

We computed several descriptive statistics on the whole set of SNPs. At the inter-population level, we assessed genetic differentiation by calculating pairwise *F*<sub>ST</sub> values using the R package *poolfstat* v2.1.1 (Gautier et al., 2022). At the intra-population level, we used in-house R scripts to compute the synonymous expected heterozygosity *H*<sub>eS</sub>, as a measure of diversity, and the ratio of non-synonymous to synonymous expected heterozygosity *H*<sub>eN/H*H*<sub>eS</sub></sub>, which provides an indirect measure of the efficacy of selection. We used a 100-block jackknife resampling approach to estimate means and standard errors. Additionally, we estimated bottleneck intensities for all invasive populations of both species using ABC analyses (See Appendix S1 for details).

Using polarized SNPs, we reported for each population and each severity category the proportion of positions for which the derived allele is absent (frequency *f* = 0), rare (*f* > 0 and *f* ≤ 0.1), common (*f* > 0.1 and *f* < 1), or fixed (*f* = 1). The proportion of rare derived alleles serves as a proxy for the masked load, whereas common and fixed alleles serve as a proxy for the realized load. Additionally, within each severity category and population, we calculated the mean derived allele frequency and estimated standard errors via a 100-block jackknife resampling approach. Within each species and severity category, mean derived allele frequencies were compared among populations using z-tests with false discovery rate (FDR; Benjamini & Hochberg, 1995) correction for multiple comparisons.

Finally, polarized SNPs were used to compute the *R*<sub>XY</sub> statistic, which allows assessment of relative excess or deficit of derived alleles within specific categories of deleterious mutations in one population compared to another (Appendix S2; Do et al., 2015; Xue et al., 2015). This statistic effectively corrects for variation due to demography using neutral mutations (synonymous SNPs in our case) as a reference point. In our analyses, population X is an invasive population, and population Y is its source population (provided in Table 1). An *R*<sub>XY</sub> value of 1 indicates an equivalent relative genetic load in both populations, whereas values greater or less than 1 suggest an excess (including fixation) or a deficit (including purge) of the genetic load in the invasive population X, respectively. We assessed the significance of any deviation from 1 using a z-score two-tailed test based on the 100-block jackknife resampling method.

# Results

## SNP-calling and genomic variation

In the target enrichment experiments, we obtained mean numbers of raw read sequences of 50,196,834 and 63,360,676 for DVV and HA, respectively (BioProject PRJNA1079689; Table S1). After trimming, 85.0% and 93.8% of the sequences were retained for DVV and HA respectively. For outgroup species, the total number of raw reads was highly variable, with a mean of 2,631,775 of which 97.8% were conserved after trimming (see detailed information for all libraries, including
outgroup species, in Table S4). Mapping, calling and filtering identified 66,274 SNPs (of which 62,034 could be polarized) and 169,755 SNPs (of which 169,102 could be polarized) within coding sequences for DVV and HA respectively (see Table S5 for details).

In the case of DVV, pairwise $F_{ST}$ computed on the full set of SNPs showed patterns very close to those expected from results of previous analyses of microsatellite datasets (Figure S4). The differentiation between the two native Mexican populations was almost negligible with a mean $F_{ST}$ estimate below 0.001. Given that the $F_{ST}$ estimates between the invasive populations and the native sample D-N-MX1 (mean $F_{ST}$ of 0.121) were consistently lower than those between the invasive populations and D-N-MX2 sample (mean $F_{ST}$ of 0.128), D-N-MX1 was therefore considered the source of the D-I-COL sample (see Table 1). The Pennsylvania and Colorado populations had a low pairwise $F_{ST}$ of 0.012, consistent with an Eastern expansion from Colorado (Figure S1), with plausible continuous gene flow. Conversely, all other $F_{ST}$ estimates were relatively high, particularly between the native and European populations (mean $F_{ST}$ value = 0.152). The highest value was observed between the two European populations at 0.169, confirming previous results indicating these populations originated from two independent introductions (Miller et al., 2005).

In the case of HA, $F_{ST}$ estimates were moderate and fairly homogeneous, ranging from 0.020 to 0.050 (Figure S4). The smallest values were found between the eastern native population H-N-CHI and both invasive populations. The $F_{ST}$ between invasive populations was slightly larger, consistent with previous results indicating independent origins for the two North American outbreaks (Lombaert et al., 2010). The highest values were observed between the western native population H-N-RUS and all other populations.

Within each species, invasive populations were characterized by reduced synonymous heterozygosity compared to their native counterparts, reflecting demographic bottlenecks encountered during the invasion process (Figure 1, x axis). One exception was the western native population of HA (H-N-RUS), which exhibited the lowest diversity within this species. Overall, the loss of diversity was more pronounced in DVV invasive populations, consistent with our ABC results, which indicate that bottlenecks were generally more intense in DVV than in HA (see Appendix S1). The ratio of non-synonymous to synonymous expected heterozygosity $H_e/N/H_s$ is typically expected to increase in populations experiencing pronounced drift, as the reduced efficacy of selection allows more non-synonymous mutations, including potentially deleterious ones, to persist. In line with this expectation, we observed a higher ratio in invasive populations (Figure 1, y axis), and a negative correlation between this ratio and synonymous expected heterozygosity within each species (Pearson's *r* = -0.95, *P* < 10⁻² for DVV; Pearson's *r* = -0.90, *P* < 10⁻¹ for HA; Figure 1). Overall, the differences between native and invasive populations were more pronounced in DVV, with invasive ranges showing lower diversities and higher $H_e/N/H_s$ ratios. This could suggest a reduced efficacy of selection, although other factors may also contribute to the observed patterns.

## Genetic load on polarized SNPs

In all populations studied and for each species, derived alleles were mostly rare (with frequencies below 0.1) or absent (Figure 2). The highest derived allele frequencies consistently showed a higher prevalence in synonymous positions compared to non-synonymous positions. In DVV, we observed a substantial loss of non-synonymous derived alleles in invasive populations compared to native ones. Although this was characterized by a reduction in both masked load (proportion of non-synonymous SNPs with rare derived alleles ranging from 0.496 to 0.523 in native populations and from 0.145 to 0.188 in invasive populations) and realized load (proportion of non-synonymous SNPs with frequent or fixed derived alleles ranging from 0.075 to 0.077 in native populations and from 0.046 to 0.051 in invasive populations), it also came with slightly higher proportions of fixed derived alleles. This trend was less pronounced for mutations inferred to be highly deleterious (LoF) than for those inferred to be moderately deleterious (missense). For HA, the disparities between native and invasive populations were considerably less pronounced (Figure 2). Compared to their native Chinese source, HA invasive populations displayed a slight
reduction in the proxy for masked load (proportion of non-synonymous SNPs with rare derived alleles: 0.455 in the Chinese population, and 0.347-0.381 in invasive populations) and a slight increase in the proxy for realized load (proportion of non-synonymous SNPs with frequent or fixed derived alleles: 0.019 in the Chinese population, and 0.024-0.027 in invasive populations). Again, this trend was less pronounced for highly deleterious mutations (Figure 2).

<figure><img src="image_3.png" /><figcaption>**Figure 1** - Ratio of non-synonymous to synonymous expected heterozygosity *H*<sub>eN</sub>/<i>H</i><sub>eS</sub> vs. the synonymous expected heterozygosity *H*<sub>eS</sub> of each *Diabrotica virgifera* *virgifera* (DVV) and *Harmonia axyridis* (HA) population. Means and standard errors (error bars) were determined from 100-block jackknife resampling. See Table 1 for population codes (underlined population codes correspond to the native populations).</figcaption></figure>

As expected, the mean derived allele frequencies within populations decline with increasing putative severity in both species (Figure 3). For DVV, synonymous and missense mean derived allele frequencies were consistently significantly lower in the invasive populations than in the native populations ($P < 10^{-4}$ for all comparisons after FDR correction), whereas no significant differences were observed between native populations or between invasive populations ($P > 0.3$ for all comparisons after FDR correction). Within each of the three severity categories, mean derived allele frequencies were virtually identical in the four HA populations ($P > 0.02$ for all comparisons after FDR correction). No differences in LoF derived allele frequencies were significant ($P > 0.2$ for all comparisons after FDR correction), for either species. Overall, derived allele frequencies were lower in HA than in DVV.
<figure><img src="image_4.png" /><figcaption>**Figure 2** - Proportion of derived alleles found in each population of *Diabrotica virgifera virgifera* (DVV) and *Harmonia axyridis* (HA) within categories of putative severity, classified by their frequency. See Table 1 for population codes (underlined population codes correspond to the native populations).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3** - Mean derived allele frequencies in all *Diabrotica virgifera virgifera* (DVV) and *Harmonia axyridis* (HA) populations in all three severity categories. Means and standard errors (error bars) were determined from 100-block jackknife resampling. See Table 1 for population codes (underlined population codes correspond to the native populations).</figcaption></figure>

Finally, the *R*<sub>XY</sub> ratio revealed no significant differences in relative frequencies of missense derived alleles between invasive populations and their respective sources in DVV (Figure 4). For putatively highly deleterious loss-of-function alleles, *R*<sub>XY</sub> values sometimes deviate sharply from 1, but the number of SNPs is small (Table S5), and most differences are not statistically significant. The only exception is found in the population from Hungary (D-I-HUN), which shows reduced relative frequencies within the putatively highly deleterious category (LoF) of alleles (Figure 4). Regarding HA, slightly higher loads were significant in three of four comparisons across severity categories between the two invasive populations (H-I-PEN and H-I-WAS) and the native Chinese population (Figure 4).

Discussion

Exploring the evolution of genetic load using genomic data has long been restricted to model species. Here, we successfully identified putatively deleterious derived alleles within coding regions in two non-model invasive insect species using a pool-seq transcriptome-based exome capture protocol (Deleury et al., 2020). By comparing native and invasive populations, we tested the hypothesis that genetic load was purged during the invasion process, which may enhance invasion success of certain populations (Sherpa & Després, 2021; Daly et al., 2023). Given that bottlenecks amplify genetic drift and weaken selection effectiveness (Crow, 1970; Glémin et al., 2003), the likelihood of purging is expected to increase for recessive, highly deleterious mutations, while moderately deleterious mutations may still experience some degree of drift, leading to higher frequencies (Agrawal & Whitlock, 2011). Our study exploring the dynamics of the genetic load in DVV and HA yielded mixed results, underscoring the nuanced interplay of evolutionary forces in invasive species.
<figure><img src="image_6.png" /></figure>

**Figure 4** - Relative mutation load estimated as the *R*<sub>XY</sub> ratio of derived alleles for moderate (missense) and high (LoF) severity categories in populations of *Diabrotica virgifera virgifera* (D) and *Harmonia axyridis* (H). A ratio below or above 1 indicates, respectively, a relative frequency deficit or excess of derived alleles of indicated severity category in the invasive population (listed first) compared to its population source (in parentheses). Means and standard errors (error bars) were determined from 100-block jackknife resampling. A star indicates that the corresponding *R*<sub>XY</sub> value is significantly different from 1 (z-score two-tailed test). For population codes, refer to Table 1 (underlined population codes correspond to native populations).

Evolution of the genetic load in *Diabrotica virgifera virgifera*

DVV exhibits marked genetic differences between native and invasive populations. Consistent
with previous studies (Ciosi et al., 2008; Lombaert et al., 2018), invasive populations, including the
long-established population in Colorado, display sharp declines in genetic diversity, consistent with
the strong bottlenecks inferred from ABC analyses. Compared to their native counterparts, all
invasive populations generally exhibit lower frequencies of derived alleles, regardless of their
potential fitness impacts, thus confirming the substantial role of genetic drift in the invasion process
of the species. While this overall reduction in genetic load within invasive populations is apparent,
it is worth emphasizing the presence of a slightly larger fixed load.

In North America, despite a strong initial loss of genetic diversity, invasive populations showed
no significant relative deficit or excess of deleterious mutations compared to their source, although
a trend toward an excess of loss-of-function variants was observed.

Interestingly, the oldest European outbreak in Central-Southeastern Europe (here a sample
from Hungary) exhibited a significant deficit of “loss-of-function” deleterious mutations compared
to its source population, while “missense” deleterious mutations did not differ significantly. This
suggests that the demographic and selective constraints in this population were effective at
reducing highly deleterious mutations through a combination of genetic drift and purifying selection,
while moderately deleterious mutations have been less affected. Such differences in evolutionary
trajectories according to mutation severity are theoretically expected (Whitlock, 2002; Caballero et
al., 2017; Dussex et al., 2023) and have been observed in other species experiencing severe
bottlenecks (Xue et al., 2015; Khan et al., 2021; Ochoa & Gibbs, 2021; Wang et al., 2023).
Finally, the invasive northwestern Italian population, despite originating from an independent introduction (Miller et al., 2005), showed patterns similar to that of the Central-Southeastern European population. However, it did not exhibit significant signals of deficit or excess of deleterious alleles, regardless of the mutation severity considered.

## Evolution of the genetic load in *Harmonia axyridis*

In HA, the overall loss of genetic diversity in both invasive populations was relatively low, in agreement with previous findings (Lombaert et al., 2010, 2011) and consistent with our ABC results suggesting that bottlenecks were relatively mild in HA invasive populations. While the putative masked genetic load decreased slightly, the putative realized load, although consistently low in all populations, showed no sign of decreasing and even exhibited a slight increase. This observation aligns with a significant trend, although subtle, towards a relative excess of deleterious mutations in invasive populations compared to their native source, regardless of putative fitness impact. This tendency to fix alleles contributing to genetic load, likely driven by drift, is significant across all populations and severity categories, except for the highly deleterious mutations in western North America (sampled in Washington).

These findings diverge from a previous study based on life history traits, which demonstrated that invasive populations of HA experienced minimal inbreeding depression while maintaining fitness comparable to outbred native populations (Facon et al., 2011). That study has greatly contributed to the hypothesis of purging genetic load as a key mechanism driving successful biological invasions (Estoup et al., 2016; Sherpa & Després, 2021; Daly et al., 2023). Given the markedly different methodologies used in the study by Facon et al. (2011) and our study, it is challenging to ascertain why our results do not align. It is conceivable that a few key genes could have disproportionate effects on invasion success, an aspect difficult to quantify using our approach, particularly in the context of the overall low genetic load previously described in this species (Tayeh et al., 2013). This issue is further complicated by the imperfect nature of our categorization, where missense mutations likely cover a broad spectrum of severity levels, potentially obscuring the influence of key genes. Moreover, very highly deleterious mutations, which are critical contributors to inbreeding depression, are expected to segregate at extremely low frequencies in large populations and may thus escape detection in our sampling due to their rarity. Finally, our exome capture protocol only targets a fraction of the exome, approximately one-third, and a broader approach may yield qualitatively different results.

## Conclusions, limits and perspectives

We tested the hypothesis that genetic load was purged during the invasion of two insect species, a crop pest (DVV) and a predator (HA). At first glance, our results offer a nuanced perspective. In the case of DVV, we were unable to detect any significant evolution of the genetic load, except for the purge of highly deleterious mutations in the invasive Central-Southeastern European population. For HA, we observed subtle evidence suggesting a tendency toward fixation, rather than toward purge, of the genetic load. Although these results are not entirely contradictory, they may stem from differences in initial genetic composition, varying demographic history, and divergent ecological niches. For instance, the two species experienced markedly different bottleneck intensities, yet our data do not allow us to establish a direct link between these differences and the observed dynamics of genetic load.

While population genomics approaches to assess the evolution of genetic load have become increasingly popular in conservation and domestication biology (Moyers et al., 2018; Bertorelle et al., 2022), our study is, to our knowledge, the first to use these methods to specifically investigate the evolution of genetic load in the context of biological invasions. However, it does not yet provide a comprehensive and definitive answer to the question addressed. One limitation lies in our reliance on allelic frequencies derived from pool-seq data, which precludes access to genotypic frequencies and therefore limits our ability to precisely estimate inbreeding depression and explore the non-linear effects of genetic load (Bataillon & Kirkpatrick, 2000). While we used complementary
statistics to mitigate this issue, future studies based on individual-level sequencing would allow for
more accurate assessments of genetic variation and its relationship to genetic load. Additionally,
to further broaden our findings and enhance our statistical power, particularly for highly deleterious
mutations, it will be essential to consider the entire exome rather than only a fraction. Another
promising avenue involves categorizing mutation severity by quantifying evolutionary constraints
within a phylogenetic framework, which could enhance our understanding of genetic load evolution,
including in non-coding regions (Davydov et al., 2010). Furthermore, extending this research to a
wider range of species would enhance our ability to identify potential universal patterns in the
evolution of genetic load during biological invasions. Finally, to test not only the occurrence of
purging in invasions but also its role in the success of invasions, we must compare the load
dynamics in introductions that fail with those that succeed – a comparison that is rarely feasible
(Zenni & Nuñez, 2013), except in controlled invasions such as in classical biological control
(Fauvergue et al., 2012).

Acknowledgments

We thank our colleagues Pamela Bruno, Charlyne Martignier-Jaccard, Pascal Mignet and Wang Su for *Diabrotica adelpha*, *Harmonia quadripunctata* and *Harmonia yedoensis* samples. We also thank Emmanuelle Murciano-Germain for administrative assistance. Sequencing was performed at the GENOTOUL GeT platform (https://get.genotoul.fr/en/).

Preprint version 4 of this article has been peer-reviewed and recommended by Peer Community In Evolutionary Biology (https://doi.org/10.24072/pci.evolbiol.100845; Rougemont, 2025).

Funding

This work was funded by grants from ANR project GENLOADICS, and INRAe SPE department.
Sampling in Russia was supported by the basic project of the Sukachev Institute of Forest SB RAS
(no. FWES-2024-0029).

**Conflict of interest disclosure**

The authors declare they comply with the PCI rule of having no financial conflicts of interest. TG and EL are both recommenders at PCI Evolutionary Biology and PCI Ecology. TG is co-founder of Peer Community In.

Data availability

Target capture pool sequencing files and all sequencing data have been deposited in Sequence
Read Archive, National Center for Biotechnology Information, under project PRJNA1079689 (see
also Table S1): https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1079689. All scripts used in this
study have been deposited at Data INRAE: https://doi.org/10.57745/ESQFDB (Lombaert, 2024).
Appendices, supplementary tables, and supplementary figures are available in the supplementary
material: https://doi.org/10.5281/zenodo.15097352 (Lombaert et al., 2024).

Author contributions

EL and ED designed the study. JB, GC, NK, TS and ST managed the collection of samples. AB and ED prepared the samples and the libraries. ED and BP developed the bioinformatics pipelines. EL and ED analyzed the data. EL and TG wrote the paper. All authors have revised and approved the final manuscript.
References

Agrawal AF, Whitlock MC (2011) Inferences about the distribution of dominance drawn from yeast gene knockout data. *Genetics*, **187**, 553–566. https://doi.org/10.1534/genetics.110.124560

Andrews S (2010) FastQC: a quality control tool for high throughput sequence data.

Ávila V, Amador C, Garcia-Dorado A (2010) The purge of genetic load through restricted panmixia in a *Drosophila* experiment. *Journal of Evolutionary Biology*, **23**, 1937–1946.
https://doi.org/10.1111/j.1420-9101.2010.02058.x

Bataillon T, Kirkpatrick M (2000) Inbreeding depression due to mildly deleterious mutations in finite populations: Size does matter. *Genetical Research*, **75**, 75–81.
https://doi.org/10.1017/S0016672399004048

Beaumont MA, Zhang WY, Balding DJ (2002) Approximate Bayesian computation in population genetics. *Genetics*, **162**, 2025–2035. https://doi.org/10.1093/genetics/162.4.2025

Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate - a practical and powerful approach to multiple testing. *Journal of the Royal Statistical Society Series B-Methodological*, **57**, 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

Bertorelle G, Raffini F, Bosse M, Bortoluzzi C, Iannucci A, Trucchi E, Morales HE, van Oosterhout C (2022) Genetic load: genomic estimates and applications in non-model animals. *Nature Reviews Genetics*, **23**, 492–503. https://doi.org/10.1038/s41576-022-00448-x

Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: A flexible trimmer for Illumina sequence data. *Bioinformatics*, **30**, 2114–2120. https://doi.org/10.1093/bioinformatics/btu170

Caballero A, Bravo I, Wang J (2017) Inbreeding load and purging: implications for the short-term survival and the conservation management of small populations. *Heredity*, **118**, 117–185.
https://doi.org/10.1038/hdy.2016.80

Charlesworth D, Morgan MT, Charlesworth B (1990) Inbreeding Depression, Genetic Load, and the Evolution of Outcrossing Rates in a Multilocus System With No Linkage. *Evolution*, **44**, 1469–1489.
https://doi.org/10.1111/j.1558-5646.1990.tb03839.x

Cingolani P, Patel VM, Coon M, Nguyen T, Land SJ, Ruden DM, Lu X (2012) Using *Drosophila melanogaster* as a model for genotoxic chemical mutational studies with a new program, *SnpSift*. *Frontiers in Genetics*, **3**, 1–9.
https://doi.org/10.3389/fgene.2012.00035

Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM (2012) A program for annotating and predicting the effects of single nucleotide polymorphisms, *SnpEff*. *Fly*, **6**, 80–92.
https://doi.org/10.4161/fly.19695

Ciosi M, Miller NJ, Kim KS, Giordano R, Estoup A, Guillemaud T (2008) Invasion of Europe by the western corn rootworm, *Diabrotica virgifera virgifera*: Multiple transatlantic introductions with various reductions of genetic diversity. *Molecular Ecology*, **17**, 3614–3627.
https://doi.org/10.1111/j.1365-294X.2008.03866.x

Coates BS, Deleury E, Gassmann AJ, Hibbard BE, Meinke LJ, Miller NJ, Petzold-Maxwell J, French BW, Sappington TW, Siegfried BD, Guillemaud T (2021) Up-regulation of apoptotic- and cell survival-related gene pathways following exposures of western corn rootworm to *B. thuringiensis* crystalline pesticidal proteins in transgenic maize roots. *BMC Genomics*, **22**, 1–27.
https://doi.org/10.1186/s12864-021-07932-4

Crnokrak P, Barrett SCH (2002) Purging the genetic load: A review of the experimental evidence. *Evolution*, **56**, 2347–2358.
https://doi.org/10.1111/j.0014-3820.2002.tb00160.x

Crow JF (1970) Genetic loads and the cost of natural selection. In: *Mathematical topics in population genetics* (eds Krickeberg K, Lewontin RC, Neyman J, Schreiber M), pp. 128–177.
Springer Berlin Heidelberg, Berlin.
https://doi.org/10.1007/978-3-642-46244-3_5

Daly EZ, Chabrerie O, Massol F, Facon B, Hess MCM, Tasiemski A, Grandjean F, Chauvat M, Viard F, Forey E, Folcher L, Buisson E, Boivin T, Baltora-Rosset S, Ulmer R, Gibert P, Thiébaut G, Pantel JH, Heger T, Richardson DM, Renault D (2023) A synthesis of biological invasion hypotheses associated with the introduction–naturalisation–invasion continuum. *Oikos*, **e09645**, 1–29.
https://doi.org/10.1111/oik.09645
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM (2021) Twelve years of SAMtools and BCFtools. *GigaScience*, **10**, 1-4. https://doi.org/10.1093/gigascience/giab008

Davydov E V., Goode DL, Sirota M, Cooper GM, Sidow A, Batzoglou S (2010) Identifying a high fraction of the human genome to be under selective constraint using GERP++. *PLoS Computational Biology*, **6**, 1-13. https://doi.org/10.1371/journal.pcbi.1001025

Deleury E, Guillemaud T, Blin A, Lombaert E (2020) An evaluation of pool-sequencing transcriptome-based exon capture for population genomics in non-model species. *bioRxiv*, **10.1101/58**, v7 peer-reviewed and recommended by PCI Genomics. https://doi.org/10.1101/583534

Do R, Balick D, Li H, Adzhubei I, Sunyaev S, Reich D (2015) No evidence that selection has been less effective at removing deleterious mutations in Europeans than in Africans. *Nature Genetics*, **47**, 126-131. https://doi.org/10.1038/ng.3186

Dussex N, Morales HE, Grossen C, Dalén L, van Oosterhout C (2023) Purging and accumulation of genetic load in conservation. *Trends in Ecology and Evolution*, **38**, 961-969. https://doi.org/10.1016/j.tree.2023.05.008

Dussex N, van der Valk T, Morales HE, Wheat CW, Díez-del-Molino D, von Seth J, Foster Y, Kutschera VE, Guschanski K, Rhie A, Phillippy AM, Korlach J, Howe K, Chow W, Pelan S, Mendes Damas JD, Lewin HA, Hastie AR, Formenti G, Fedrigo O, Guhlin J, Harrop TWR, Le Lec MF, Dearden PK, Haggerty L, Martin FJ, Kodali V, Thibaud-Nissen F, Iorns D, Knapp M, Gemmell NJ, Robertson F, Moorhouse R, Digby A, Eason D, Vercoe D, Howard J, Jarvis ED, Robertson BC, Dalén L (2021) Population genomics of the critically endangered *kākāpō*. *Cell Genomics*, **100002**, 1-10. https://doi.org/10.1016/j.xgen.2021.100002

Eben A, de los Monteros AE (2013) Tempo and mode of evolutionary radiation in Diabroticina beetles (genera *Acalymma*, *Cerotoma*, and *Diabrotica*). *ZooKeys*, **332**, 207-231. https://doi.org/10.3897/zookeys.332.5220

Enders M, Hutt MT, Jeschke JM (2018) Drawing a map of invasion biology based on a network of hypotheses. *Ecosphere*, **9**, https://doi.org/10.1002/ECS2.2146

Estoup A, Ravigné V, Hufbauer R, Vitalis R, Gautier M, Façon B (2016) Is There A Genetic Paradox of Biological Invasion? *Annual Review of Ecology Evolution and Systematics*, **47**, 51-72. https://doi.org/10.1146/annurev-ecolsys-121415-032116

Façon B, Hufbauer RA, Tayeh A, Loiseau A, Lombaert E, Vitalis R, Guillemaud T, Lundgren JG, Estoup A (2011) Inbreeding depression is purged in the invasive insect *Harmonia axyridis*. *Current Biology*, **21**, 424-427. https://doi.org/10.1016/j.cub.2011.01.068

Fauvergue X, Vercken E, Malausa T, Hufbauer RA (2012) The biology of small, introduced populations, with special reference to biological control. *Evolutionary Applications*, **5**, 424-443. https://doi.org/10.1111/j.1752-4571.2012.00272.x

Fountain T, Duvaux L, Horsburgh G, Reinhardt K, Butlin RK (2014) Human-facilitated metapopulation dynamics in an emerging pest species, *Cimex lectularius*. *Molecular Ecology*, **23**, 1071-1084. https://doi.org/10.1111/mec.12673

Garrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. *arXiv*, **1207.3907**, 1-9. https://doi.org/10.48550/arXiv.1207.3907

Gautier M, Vitalis R, Flori L, Estoup A (2022) f-Statistics estimation and admixture graph construction with Pool-Seq or allele count data using the R package poolstat. *Molecular Ecology Resources*, **22**, 1394-1416. https://doi.org/10.1111/1755-0998.13557

Glémin S (2003) How are deleterious mutations purged? Drift versus nonrandom mating. *Evolution*, **57**, 2678-2687. https://doi.org/10.1111/j.0014-3820.2003.tb01512.x

Glémin S, Ronfort J, Bataillon T (2003) Patterns of Inbreeding Depression and Architecture of the Load in Subdivided Populations. *Genetics*, **165**, 2193-2212. https://doi.org/10.1093/genetics/165.4.2193

Gouzy J, Carrere S, Schiex T (2009) FrameDP: Sensitive peptide detection on noisy matured sequences. *Bioinformatics*, **25**, 670-671. https://doi.org/10.1093/bioinformatics/btp024
Gray ME, Sappington TW, Miller NJ, Moeser J, Bohn MO (2009) Adaptation and Invasiveness of Western Corn Rootworm: Intensifying Research on a Worsening Pest. *Annual Review of Entomology*, **54**, 303–321. https://doi.org/10.1146/annurev.ento.54.110807.090434

Grossen C, Guillaume F, Keller LF, Croll D (2020) Purging of highly deleterious mutations through severe bottlenecks in Alpine ibex. *Nature Communications*, **11**, 1001. https://doi.org/10.1038/s41467-020-14803-1

Guillemaud T, Ciosi M, Lombaert E, Estoup A (2011) Biological invasions in agricultural settings: Insights from evolutionary biology and population genetics. *Comptes Rendus Biologies*, **334**, 237–246. https://doi.org/10.1016/j.crvi.2010.12.008

Henn BM, Botigué LR, Bustamante CD, Clark AG, Gravel S (2015) Estimating the mutation load in human genomes. *Nature Reviews Genetics*, **16**, 333–343. https://doi.org/10.1038/nrg3931

Keightley PD, Jackson BC (2018) Inferring the Probability of the Derived vs. the Ancestral Allelic State at a Polymorphic Site. *Genetics*, **209**, 897–906. https://doi.org/10.1534/genetics.118.301120

Khan A, Patel K, Shukla H, Viswanathan A, van der Valk T, Borthakur U, Nigam P, Zachariah A, Jhala Y V., Kardos M, Ramakrishnan U (2021) Genomic evidence for inbreeding depression and purging of deleterious genetic variation in Indian tigers. *Proceedings of the National Academy of Sciences*, **118**, e2023018118. https://doi.org/10.1073/pnas.2023018118

Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. *arXiv*, **1303.3997**, 1–3. https://doi.org/10.48550/arXiv.1303.3997

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R (2009) The Sequence Alignment/Map format and SAMtools. *Bioinformatics*, **25**, 2078–2079. https://doi.org/10.1093/bioinformatics/btp352

Lohmueller KE, Indap AR, Schmidt S, Boyko AR, Hernandez RD, Hubisz MJ, Sninsky JJ, White TJ, Sunyaev SR, Nielsen R, Clark AG, Bustamante CD (2008) Proportionally more deleterious genetic variation in European than in African populations. *Nature*, **451**, 994–997. https://doi.org/10.1038/nature06611

Lombaert E (2024) Scripts used for data analyzing in the manuscript: “Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species.” *Recherche Data Gouv*, v1. https://doi.org/10.57745/ESQFDB

Lombaert, E., Blin, A., Porro, B., guillemaud, thomas, Bernal, J. S., Gary, C., Kirichenko, N., Sappington, T., Toepfer, S., & deleury, emeline. (2025). Supplementary material for the article "Unraveling genetic load dynamics during biological invasion: insights from two invasive insect species" published in Peer Community Journal. Zenodo. https://doi.org/10.5281/zenodo.15097352

Lombaert E, Ciosi M, Miller NJ, Sappington TW, Blin A, Guillemaud T (2018) Colonization history of the western corn rootworm (Diabrotica virgifera virgifera) in North America: insights from random forest ABC using microsatellite data. *Biological Invasions*, **20**, 665–677. https://doi.org/10.1007/s10530-017-1566-2

Lombaert E, Guillemaud T, Cornuet JM, Malausa T, Facon B, Estoup A (2010) Bridgehead effect in the worldwide invasion of the biocontrol harlequin ladybird. *PLoS One*, **5**, e9743. https://doi.org/10.1371/journal.pone.0009743

Lombaert E, Guillemaud T, Lundgren J, Koch R, Facon B, Grez A, Loomans A, Malausa T, Nedved O, Rhule E, Staverlokk A, Steenberg T, Estoup A (2014) Complementarity of statistical treatments to reconstruct worldwide routes of invasion: the case of the Asian ladybird *Harmonia axyridis*. *Molecular ecology*, **23**, 5979–5997. https://doi.org/10.1111/mec.12989

Lombaert E, Guillemaud T, Thomas CE, Lawson Handley LJ, Li J, Wang S, Pang H, Goryacheva I, Zakharov IA, Jousselin E, Poland RL, Migeon A, Van Lenteren J, De Clercq P, Berkvens N, Jones W, Estoup A (2011) Inferring the origin of populations introduced from a genetically structured native range by approximate Bayesian computation: case study of the invasive ladybird *Harmonia axyridis*. *Molecular Ecology*, **20**, 4654–4670. https://doi.org/10.1111/j.1365-294X.2011.05322.x
Lynch M, Conery IJ, Burger R (1995) Mutation Accumulation and the Extinction of Small Populations. *The American Naturalist*, **146**, 489–518. https://doi.org/10.1086/285812

Makino T, Rubin C-J, Carneiro M, Axelsson E, Andersson L, Webster MT (2018) Elevated proportions of deleterious genetic variation in domestic animals and plants. *Genome biology and evolution*, **10**, 276–290. https://doi.org/10.1093/gbe/evy004

Marchini GL, Sherlock NC, Ramakrishnan AP, Rosenthal DM, Cruzan MB (2016) Rapid purging of genetic load in a metapopulation and consequences for range expansion in an invasive plant. *Biological Invasions*, **18**, 183–196. https://doi.org/10.1007/s10530-015-1001-5

Marsden CD, Ortega-Del Vecchyo D, O'Brien DP, Taylor JF, Ramirez O, Vilà C, Marques-Bonet T, Schnabel RD, Wayne RK, Lohmueller KE (2016) Bottlenecks and selective sweeps during domestication have increased deleterious genetic variation in dogs. *Proceedings of the National Academy of Sciences*, **113**, 152–157. https://doi.org/10.1073/pnas.1512501113

Miller N, Estoup A, Toepfer S, Bourguet D, Lapchin L, Derridj S, Kim KS, Reynaud P, Furlan L, Guillemaud T (2005) Multiple transatlantic introductions of the western corn rootworm. *Science*, **310**, 992. https://doi.org/10.1126/science.1115871

Moyers BT, Morrell PL, McKay JK (2018) Genetic costs of domestication and improvement. *Journal of Heredity*, **109**, 103–116. https://doi.org/10.1093/jhered/esx069

Mullarkey A a., Byers DL, Anderson RC (2013) Inbreeding depression and partitioning of genetic load in the invasive biennial *Alliaria petiolata* (Brassicaceae). *American Journal of Botany*, **100**, 509–518. https://doi.org/10.3732/ajb.1200403

Nei M, Maruyama T, Chakraborty R (1975) Bottleneck Effect and Genetic Variability in Populations. *Evolution*, **29**, 1–10. https://doi.org/10.2307/2407137

Ning Z, Cox AJ, Mullikin JC (2001) SSAHA: a fast search method for large DNA databases. *Genome Research*, **11**, 1725–1729. https://doi.org/10.1101/gr.194201

Ochoa A, Gibbs HL (2021) Genomic signatures of inbreeding and mutation load in a threatened rattlesnake. *Molecular Ecology*, **30**, 5454–5469. https://doi.org/10.1111/mec.16147

Parisod C, Trippi C, Galland N (2005) Genetic variability and founder effect in the pitcher plant *Sarracenia purpurea* (Sarraceniaceae) in populations introduced into Switzerland: From inbreeding to invasion. *Annals of Botany*, **95**, 277–286. https://doi.org/10.1093/aob/mci023

R Core Team (2021) R: A language and environment for statistical computing.

Robinson J, Kyriazis CC, Yuan SC, Lohmueller KE (2023) Deleterious Variation in Natural Populations and Implications for Conservation Genetics. *Annual Review of Animal Biosciences*, **11**, 93–114. https://doi.org/10.1146/annurev-animal-080522-093311

Robinson JA, Ortega-Del Vecchyo D, Fan Z, Kim BY, Vonholdt BM, Marsden CD, Lohmueller KE, Wayne RK (2016) Genomic Flatlining in the Endangered Island Fox. *Current Biology*, **26**, 1183–1189. https://doi.org/10.1016/j.cub.2016.02.062

Rougemont Q (2025) The genetic load of invasive population: how little do we know? *Peer Community in Evolutionary Biology*, 100845. https://doi.org/10.24072/pci.evolbiol.100845

Roy HE, Brown PMJ, Adriaens T, Berkvens N, Borges I, Clusella-Trullas S, Comont RF, De Clercq P, Eschen R, Estoup A, Evans EW, Facon B, Gardiner MM, Gil A, Grez AA, Guillemaud T, Haelewaters D, Herz A, Honek A, Howe AG, Hui C, Hutchison WD, Kenis M, Koch RL, Kulfan J, Lawson Handley L, Lombaert E, Loomans A, Losey J, Lukashuk AO, Maes D, Magro A, Murray KM, Martin GS, Martinkova Z, Minnaar IA, Nedved O, Orlova-Bienkowskaja MJ, Osawa N, Rabitsch W, Ravn HP, Rondoni G, Rorke SL, Ryndevich SK, Saethre M-G, Sloggett JJ, Soares AO, Stals R, Tinsley MC, Vandereycken A, van Wielink P, Viglášová S, Zach P, Zakharov IA, Zaviezo T, Zhao Z (2016) The harlequin ladybird, *Harmonia axyridis*: global perspectives on invasion history and ecology. *Biological Invasions*, **18**, 997–1044. https://doi.org/10.1007/s10530-016-1077-6

Sakai AK, Allendorf FW, Holt JS, Lodge DM, Molofsky J, With KA, Baughman S, Cabin RJ, Cohen JE, Ellstrand NC, McCauley DE, O'Neil P, Parker IM, Thompson JN, Weller SG (2001) The population biology of invasive species. *Annual Review of Ecology and Systematics*, **32**, 305–332. https://doi.org/10.1146/annurev.ecolsys.32.081501.114037
Schubert M, Jónsson H, Chang D, Der Sarkissian C, Ermini L, Ginolhac A, Albrechtsen A, Dupanloup I, Foucal A, Petersen B, Fumagalli M, Raghavan M, Seguin-Orlando A, Korneliussen TS, Velazquez AM V, Stenderup J, Hoover CA, Rubin C-J, Alfarhan AH, Alquraishi SA, Al-Rasheid KAS, MacHugh DE, Kalbfleisch T, MacLeod JN, Rubin EM, Sicheritz-Ponten T, Andersson L, Hofreiter M, Marques-Bonet T, Gilbert MTP, Nielsen R, Excoffier L, Willerslev E, Shapiro B, Orlando L (2014) Prehistoric genomes reveal the genetic foundation and cost of horse domestication. *Proceedings of the National Academy of Sciences of the United States of America*, **111**, E5661-9. https://doi.org/10.1073/pnas.1416991111

Seebens H, Blackburn TM, Dyer E, Genovesi P, Hulme PE, Jeschke JM, Pagad S, Pyšek P, Winter M, Arianoutsou M, Bacher S, Blasius B, Brundu G, Capinha C, Celesti-Grapow L, Dawson W, Dullinger S, Fuentes N, Jäger H, Kartesz J, Kenis M, Kreft H, Kühn I, Lenzner B, Liebhold A, Mosena A, Moser D, Nishino M, Pearman D, Pergl J, Rabitsch W, Rojas-Sandoval J, Roques A, Rorke S, Rossinelli S, Roy HE, Scalera R, Schindler S, Štajerová K, Tokarska-Guzik B, Kleunen M van, Walker K, Weigelt P, Yamanaka T, Essl F (2017) No saturation in the accumulation of alien species worldwide. *Nature Communications*, **8**, 14435.
https://doi.org/10.1038/ncomms14435

Sherpa S, Després L (2021) The evolutionary dynamics of biological invasions: A multi-approach perspective. *Evolutionary Applications*, 1-23. https://doi.org/10.1111/eva.13215

Simberloff D (2009) The Role of Propagule Pressure in Biological Invasions. *Annual Review of Ecology Evolution and Systematics*, **40**, 81-102.
https://doi.org/10.1146/annurev.ecolsys.110308.120304

Tayeh A, Estoup A, Hufbauer RA, Ravigne V, Goryacheva I, Zakharov IA, Lombaert E, Facon B (2013) Investigating the genetic load of an emblematic invasive species: the case of the invasive harlequin ladybird *Harmonia axyridis*. *Ecology and Evolution*, **3**, 864-871.
https://doi.org/10.1002/ece3.490

Tomaszewska W, Escalona HE, Hartley D, Li J, Wang X, Li H Sen, Pang H, Ślipiński A, Zwick A (2021) Phylogeny of true ladybird beetles (Coccinellidae: Coccinellini) reveals pervasive convergent evolution and a rapid Cenozoic radiation. *Systematic Entomology*, **46**, 611-631.
https://doi.org/10.1111/syen.12479

Vogel H, Schmidtberg H, Vilcinskas A (2017) Comparative transcriptomics in three ladybird species supports a role for immunity in invasion biology. *Developmental and Comparative Immunology*, **67**, 452-456. https://doi.org/10.1016/j.dci.2016.09.015

Wang X, Peischl S, Heckel G (2023) Demographic history and genomic consequences of 10,000 generations of isolation in a wild mammal. *Current Biology*, **33**, 2051-2062.e4.
https://doi.org/10.1016/j.cub.2023.04.042

Wang MS, Zhang JJ, Guo X, Li M, Meyer R, Ashari H, Zheng ZQ, Wang S, Peng MS, Jiang Y, Thakur M, Suwannapoom C, Esmailizadeh A, Hirimuthugoda NY, Zein MSA, Kusza S, Kharrati-Koopae H, Zeng L, Wang YM, Yin TT, Yang MM, Li ML, Lu XM, Lasagna E, Ceccobelli S, Gunwardana HGTN, Senasig TM, Feng SH, Zhang H, Bhuiyan AKFH, Khan MS, Silva GLLP, Thuy LT, Mwai OA, Ibrahim MNM, Zhang G, Qu KX, Hanotte O, Shapiro B, Bosse M, Wu DD, Han JL, Zhang YP (2021) Large-scale genomic analysis reveals the genetic cost of chicken domestication. *BMC Biology*, **19**, 1-16. https://doi.org/10.1186/s12915-021-01052-x

Whitlock MC (2002) Selection, load and inbreeding depression in a large metapopulation. *Genetics*, **160**, 1191-1202. https://doi.org/10.1093/genetics/160.3.1191

Williamson M, Fitter A (1996) The varying success of invaders. *Ecology*, **77**, 1661-1666.
https://doi.org/10.2307/2265769

Xue Y, Prado-martinez J, Sudmant PH, Narasimhan V, Ayub Q, Szpak M, Frandsen P, Chen Y, Yngvadottir B, Cooper DN, Manuel M De, Hernandez-rodriguez J, Lobon I, Siegismund HR, Pagani L, Quail MA, Hvilsom C, Mudakikwa A, Eichler EE, Cranfield MR, Marques-bonet T (2015) Mountain gorilla genomes reveal the impact of long-term population decline and inbreeding. *Science*, **348**, 242-245. https://doi.org/10.1126/science.aaa3952

Zayed A, Constantin S a., Packer L (2007) Successful biological invasion despite a severe genetic load. *PLoS ONE*, **2**. https://doi.org/10.1371/journal.pone.0000868
Zenni RD, Nuñez MA (2013) The elephant in the room: The role of failed invasions in understanding invasion biology. *Oikos*, **122**, 801–815. https://doi.org/10.1111/j.1600-0706.2012.00254.x