Population genetics of the endangered narrowly endemic Island Marble butterfly (Euchloe ausonides insulanus)

Kara S. Jones1*, Aaron W. Aunins1, Colleen C. Young1, Robin L. Johnson1 and Cheryl L. Morrison1

  1. 1U.S. Geological Survey, Eastern Ecological Science Center at Leetown Research Laboratory, 11469 Leetown Road, Kearneysville, WV 25430

*Corresponding author:

Kara S. Jones, ksjones{at}usgs.gov

bioRxiv preprint DOI: https://doi.org/10.1101/2025.05.08.652852

Posted: May 14, 2025, Version 1

Copyright: This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available for use under a CC0 license

Abstract

The Island Marble butterfly (Euchloe ausonides insulanus) is an endangered species endemic to the San Juan Islands off the coast of Washington State, United States, and British Columbia, Canada. The species was thought to be extinct for ∼90 years before it was rediscovered at American Camp, San Juan Island National Historical Park in 1998. Here, we report the results of the first population genetic analyses for insulanus, using DNA collected non-invasively from individuals in the last known stronghold for the species. We used DNA extracted from meconium, larval exuviae, and natural mortalities to generate and test thirteen new microsatellite markers to estimate genetic diversity, population structure, and kinship. We assembled and annotated mitochondrial genomes, which were used alongside museum specimens of insulanus collected ∼100 years ago from Vancouver Island, and other members of the E. ausonides species complex, to infer the evolutionary history of the species. The results indicated that insulanus experiences low heterozygosity, a small effective population size (Ne), and low allelic diversity. High levels of inbreeding were found in some individuals, but inbreeding was uneven across the population. No population structure or partitioning of genetic variation by host plant was detected. The mitogenomes of extant insulanus were all identical and modern samples showed a loss of allelic diversity compared to insulanus from museums. Extant insulanus formed a clade with museum specimens and we identified multiple putatively diagnostic alleles to differentiate insulanus from other subspecies. Based on these results, we outline considerations for species management and genetic monitoring.

Introduction

The Island Marble butterfly (Euchloe ausonides insulanus), hereafter referred to as insulanus, is a small butterfly in the family Pieridae that is endemic to the San Juan Islands off the coast of Washington State, United States and British Columbia, Canada. The last specimen was collected from Gabriola Island, British Columbia, in 1908, and was believed to be extinct for ∼90 years before the butterfly was rediscovered at American Camp, San Juan Island National Historical Park (SAJH NHP) in 1998 (Shepard, 2000). Extensive survey efforts identified five areas of occupied habitat across two islands (San Juan Island and Lopez Island), but subsequent habitat loss has led to widespread declines in the distribution of the species (Lambert, 2011). By 2012, the species distribution had contracted to a single population mainly within the boundaries of American Camp. Due to the declining status of insulanus, an emergency captive-rearing effort to sustain the last remaining population began in 2013. In 2020, insulanus was listed as an endangered species (U.S. Fish and Wildlife Service 2020; 85 FR 26786) under the Endangered Species Act (ESA 1973, as amended), and 812 acres on the south end of San Juan Island were designated critical habitat (50 CFR 17.95(i)).

Several life history traits may contribute to the current precarity of insulanus. Based on Lambert’s (2011) detailed natural history of the species, oviposition occurs on three larval host plants in the mustard family (Brassicaceae), each of which primarily occurs in discrete habitats: Brassica rapa (grasslands), Sisymbrium altissimum (sand dunes), and Lepidium virginicum var. menziesii (tidal zones). Larvae develop through five instars, taking ∼38 days to progress from egg to pupa, after which pre-pupal larvae “wander” up to four meters to a different species of plant to pupate for 10-11 months. Over 50% of eggs are lost prior to hatching, approximately half from predation on host plants by deer. Larvae are vulnerable to death from starvation, particularly when feeding on host plants growing in marginal habitat, where the host plants may senesce before larvae reach the pre-pupal stage. Furthermore, insulanus is univoltine (i.e., has non-overlapping generations), making the species particularly vulnerable to the effects of storms and localized climate conditions during the long winter diapause as pupae (McDermott Long et al., 2017). The larval host plants S. altissimum and L. virginicum var. menziesii, are located primarily in sand dunes and tidal habitats, respectively, which increases the vulnerability of host plants to wind, wave action, and other disturbances. This can lead to fewer oviposition sites, higher larval mortality, or higher pupal mortality, depending on the time of year (Lambert, 2011). Hence the combined stressors of loss and degradation of habitat, storm surges, direct and incidental predation, and limited larval host plant availability, may have contributed to the extirpation of insulanus from Vancouver, Gabriola, and Lopez Islands and contraction of the species’ range on San Juan Island (Shepard 2000; U.S. Fish and Wildlife Service 2023a).

Previous reports examining the status of insulanus identified several information gaps that could be addressed by genetic data, including neutral genetic diversity, effective population size, and population structure (U.S. Fish and Wildlife Service 2023a, 2023b, 2023c). However, there are several challenges and limitations to consider when generating genetic markers for an imperiled species. For example, non-invasive sampling techniques are usually required to minimize harm to the organism being studied when collecting DNA samples. In the case of insulanus, DNA collection is mostly limited to lower-quality sources, such as frass or exuviae shed between larval instars, potentially resulting in fewer loci than data sets obtained by extracting DNA from fresh tissue samples (Schultz et al., 2022). Non- invasively sampled DNA can also lead to an increased number of genotyping errors due to allele dropout from partially degraded templates (Zhan et al., 2010), errors introduced from an increased number of amplification cycles (Arantes et al., 2023), and sequence contamination with DNA from non-target organisms and the surrounding environment (Balmori-de la Puente et al., 2023; Sittenthaler et al., 2023). Moreover, non-invasively collected samples are usually collected opportunistically, and therefore the individuals included in genetic analyses may not be representative of the population, which can bias results (Waits and Leberg, 2000). Hence, genetic studies of threatened and endangered species require careful consideration at each step to correct errors and reduce bias introduced through suboptimal sampling procedures.

In this study, we assessed the quality of DNA that could be recovered from available non- invasive sampling templates (e.g., meconium, chrysalides, exuviae, and frass) collected from insulanus being housed at the SAJH NHP Captive Rearing Facility. We developed a suite of microsatellite markers, which we assessed for information content and biases that may have been introduced by sampling constraints, before using the markers to estimate genetic diversity, levels of inbreeding, and effective population size. We used kinship analyses to identify cryptic relatedness and create a reduced- relatedness data set to explore whether the population was structured by larval host plant use. We also used mitochondrial genomes to explore the evolutionary history of the species in relation to other subspecies in the complex. Finally, we place the results in the broader conservation planning initiative for insulanus, particularly regarding captive rearing and translocation efforts.

Materials and Methods

Sample collection

Samples of insulanus used for this study were collected non-lethally by researchers at the SAJH NHP Captive Rearing Facility from 2013 through 2019, prior to insulanus’ listing as an endangered species. A variety of sample types were collected for DNA, including tissue from natural mortalities (e.g., eggs, larvae), and discarded solids (e.g., chrysalides, exuviae, frass) and liquids (e.g., meconium) (Table 1). Tissues and solid samples were collected using gloved hands and clean forceps, and samples were submerged in lysis buffer or 95% ethanol in 2 μl sterile screwcap tubes. Liquid meconium samples were spotted on FTA cards (Cytiva, Wilmington, Delaware), and the cards were stored in sealable plastic storage bags with desiccant. Samples were shipped with an ice pack to USGS EESC-LRL for further processing. Sample metadata are listed in Supplementary Table S1.

Table 1

Types of non-invasive samples collected for DNA. All samples were collected either dry, in 95% ethanol, or in Qiagen DNEasy Cell Lysis Buffer (Germantown, Maryland), except where noted.

Microsatellite marker development and genotyping

DNA was extracted from the meconium of a captive insulanus butterfly (Eai-S43) using the Qiagen DNeasy Blood and Tissue Kit (Germantown, Maryland) following the manufacturer’s protocol. A genomic DNA library was prepared for Illumina sequencing using the TruSeq Nano kit (San Diego, California), and 2×150 bp paired-end (PE) and sequenced on a NextSeq500 (San Diego, California) at the USGS EESC-LRL. The resulting sequences were quality trimmed using Qiagen CLC Genomics Workbench (CLCGW) v9.5.2 (Qiagen, 2016) using a quality limit value of 0.05 (corresponding to a Phred score of ∼13), and a minimum sequence length of 15 bp. Cleaned sequences were then de novo assembled (i.e., assembled without a reference) into contigs using CLCGW with the following parameters: automatic bubble and word size, minimum contig length of 200, scaffold assembly enabled, and auto-detection of paired distances.

We screened contigs for di-, tri-, tetra-, penta-, and hexanucleotide microsatellite repeat motifs with a minimum repeat length of five using the program QDD v3 (Meglécz et al., 2010). Loci were aligned to the GenBank Nucleotide database (Benson et al., 2017) using BLAST (Camacho et al., 2009) and discarded if they could be assigned to a probable contaminant (e.g., fungi or bacteria). Primer sequences for the microsatellite loci were identified using Primer 3 within QDD, with forward primers modified to contain a universal M13 tail to enable the incorporation of a fluorescent dye label during PCR for fragment analysis.

Seventy-four microsatellite primer pairs were chosen for testing. Microsatellite DNA amplification consisted of 5-40 ng of genomic DNA, 5X PCR buffer (10 mM Tris-HCl, pH 8.3, 50 mM KCl), 2 mM MgCl2, 0.25 mM dNTPs, 0.25 µM forward or long M13 tagged primer, 0.5 µM reverse or short primer, 0.1 µM M13 fluorescently labeled primer, 0.13 U Bovine Serum Albumin (BSA) and 0.1 U Taq DNA polymerase (Promega, Madison, WI, USA) in a final reaction volume of 15 µl. Each locus was amplified separately on a T100 Thermal Cycler (Bio-Rad, Hercules, California) using the following conditions: initial denaturing at 94°C for 15 min, 29 cycles of 94° C for 1 min, 60°C for 45 sec, 72 ° C for 45 sec, followed by 5 cycles of 94° C for 1 min, 53° C for 45 sec, 72 ° C for 45 sec and a final extension at 72 ° C for 10 min. Fragment length polymorphism was assessed using a ThermoFisher ABI 3130 Capillary Sequencer (Waltham, Massachusetts) and ThermoFisher’s GeneMapper v4.1 software (Applied Biosystems, 2009).

Thirteen loci amplified consistently across samples and were polymorphic. These primer sets were ordered with ABI fluorescent dyes and without M13 tails and were arranged into multiplex PCR reactions (one duplex, one triplex, two quadruplexes; Table S2) using Multiplex Manager (Holleley and Geerts 2009). These new multiplex PCRs were used to genotype additional samples, most of which were from meconium preserved on FTA cards, following the same scoring on the ABI 3130.

Summary statistics for each microsatellite locus, including the number of alleles (Na), effective number of alleles (Ae), observed heterozygosity (Ho), expected heterozygosity (He), fixation index (F), and tests for conformance to Hardy-Weinberg equilibrium (HWE), were calculated in GenAlEx (Peakall and Smouse, 2012).

Microsatellite kinship, inbreeding, and relatedness estimates

Population genetic analyses can be biased when using relatively few markers and/or markers with low information content. We therefore used the pic_calc function from the PopGenUtils package (Tourvas, 2024) in R 4.4.1 (R Core Team, 2024), to calculate the polymorphic information content (PIC) of each locus. We also used the pid_calc function to calculate the probability of identity for two unrelated (P(ID)) or two related (P(ID)sib) individuals. PIC measures how much variation occurs at a locus within a population, providing an estimate of whether a marker can detect genetic differences between individuals, and P(ID)/P(ID)sib estimates the probability that two individuals will have the same genotype (Waits et al., 2001Serrote et al., 2020). One locus (Locus-71) was dropped from all subsequent analyses because of low PIC and high P(ID)/P(ID)sib.

We estimated sibship among individuals, effective population size (Ne), and allelic error rates using COLONY 2.0.7 (Wang, 20092018Jones and Wang, 2010). Since the life cycle of insulanus lasts for one year and adults die after breeding, samples were grouped into three offspring-parent cohorts. Genotypes from 2017, 2018, and 2019 were each set as offspring and analyzed with candidate parent genotypes from the previous year (e.g., 2017 samples as offspring were analyzed with 2016 samples as parents). The sex of the samples was unknown, so the same set of genotypes was used for both maternal and paternal candidates. To ensure consistency across results, we ran each cohort with the full likelihood method using the “very high” precision setting and a “very long” run length (i.e., the most precise and longest lengths possible for the program). Each cohort was run for 10 replicates using the following settings: dioecious species (2), inbreeding allowed (1), polygamous males (1) and females (1), no clone inference (0), full sibship size scaling (1), no sibship size prior (0), unknown population allele frequency (0), allele frequency updates (1). A starting error rate of 1% was used for estimating allelic dropout and false alleles, and the probability that an actual mother or father was included in candidate parent genotypes was set at 1%.

We used EMIBD9 1.1.0.0 to estimate inbreeding coefficients (F; i.e., the probability that two alleles at a specific locus in an individual are identical by descent) for each individual and to calculate pairwise relatedness (r) for all individuals across all years (Wang, 2022). To determine whether significant variance in inbreeding was present among individuals, we calculated g2 with the g2_microsats function in the inbreedR package, using 10,000 permutations to obtain the p-value and 10,000 bootstraps to calculate 95% confidence intervals (Stoffel et al., 2016). To determine whether the insulanus population was being structured by host plant use, we evaluated the difference in mean pairwise relatedness between individuals found on the same host plants versus different host plants with two-sided t-tests using the Infer package in R (Couch et al., 2021).

Since the inclusion of highly related individuals can lead to clustering by kinship and obfuscation of other patterns in the data (Falush et al., 2007), we created a reduced-relatedness data set by removing all but one individual from each full sibship group with ≥0.9 inclusive probability assigned by COLONY. We performed a principal components analysis (PCA) (Patterson et al., 2006) on this reduced data set (N = 193) with the dudi.pca function of the adegenet 2.1.10 (Jombart, 2008) package in R to explore whether there was structure in the data that could be attributable to differential host plant use. Additionally, we ran Structure 2.3.4 (Pritchard et al., 2000) for population clusters (K) of 2, 3, and 4 using a burn-in of 100,000 steps and 500,000 MCMC steps. The MCMC chain was considered to have convergence when the likelihood values had minimal variance between steps. Individuals were split into groups by year, and each year and value of K was run ten times independently. Data was summed across runs using Clumpak (Kopelman et al., 2015). The PCA and Structure results were visualized in R using tools in the tidyverse set of packages (Wickham et al., 2019).

Mitogenome sequencing, assembly, and annotation

We compared shotgun sequencing with Long Range PCR (LR-PCR) to assess which method was more effective at recovering mitogenomes with low template quantity and/or quality (Ponce and Micol, 1992). For shotgun sequencing, seven insulanus samples and two samples from other Euchloe ausonides subspecies were prepared for shotgun genomic sequencing on an Illumina MiSeq (San Diego, California) using the New England Biolabs NEBNext Ultra II DNA Library Prep Kit for Illumina (Ipswich, Massachusetts). Samples were dual indexed with combinatorial indices and paired-end sequenced at EESC-LRL across two 500-cycle MiSeq cartridges. For LR-PCR, we used three insulanus and one E. a. subsp. from the same pool of individuals as above. LR-PCR primers were adapted from Park et al. (2012) to better match Euchloe sequences in GenBank (Table S3). Each mitogenomic fragment was amplified in a 50 µl PCR reaction consisting of 0.5 µl TaKaRa LA Taq (San Jose, California; 5U/µl stock concentration), 5 µl 10X LA PCR Buffer, 5 µl 25 nM MgCl2, 8 µl dNTP mixture (2.5 mM each stock concentration), 1 µl forward and 1 µl reverse primers (each at 10 µM stock concentration) 2 µl DNA template, and 27.5 µl PCR grade water. Thermal cycling parameters for each of the fragments were 94°C for 1 min, followed by 35 cycles of 98°C for 10 s, relevant annealing temperature (AT; degrees C) from Table S3 for 15 min, 72°C for 30 s, and a final extension at 72°C for 10 minutes. These fragments were purified and pooled in equimolar concentrations and then subjected to the same library preparation protocol as the shotgun libraries for Illumina sequencing (as described in the section Microsatellite marker development and genotyping).

Reads from WGS sequencing and LR-PCR were de novo assembled with MitoZ 3.6 (Meng et al., 2019) using Megahit (Li et al., 2015) for the assembler, and substrings of length k (k-mers) of 59, 79, 99, 119, 141. For comparison, a subset of individuals (EaC3, 17SI21, Eai-163) was also assembled with MitoFinder 1.4.1 (Allio et al., 2020), using the annotated Anthocharis bambusarum mitogenome (NCBI RefSeq: NC_025274) as reference (Ebdon et al., 2022). The resulting mitogenomes were annotated with MITOS2 2.1.9 and mitochondrial tRNA finder (MiTFi) on the public server Galaxy (usegalaxy.org), using a genetic code = 5 for translation and RefSeq89 Metazoa as a reference (Al Arab et al., 2017Afgan et al., 2018Donath et al., 2019). Additional tRNA annotation was performed using tRNAscan- SE 2.0 (Chan et al., 2021). Annotations were compared with existing mitochondrial genomes to ensure they had a similar length and structure. Protein coding genes were also verified by checking for the presence of start and stop codons, absence of stop codons within the gene, and similarity to genes found in other butterfly mitogenomes found in GenBank.

Mitogenome phylogeny

To determine how the insulanus mitogenome differed from those of other subspecies, we obtained whole mitogenome sequences that had been assembled for an unrelated Euchloe project (Zhang et al., University of Texas Southwestern Medical Center, unpublished data, 2023). The final data set had mitogenomes from six individuals of E. a. ausonides, six E. a. coloradensis, three E. a. mayi, four E. a. palaeoreios, eight E. a. transmontana, and five additional individuals of insulanus derived from museum specimens collected from Vancouver Island circa 1898-1904 (Table S4). All assemblies were rotated to the same starting position using Cyclic DNA Sequence Aligner (Fernandes et al., 2009), aligned using MUSCLE 3.8.31 (Edgar, 2004) and visually assessed for consistency in AliView 1.28 (Larsson, 2014). Single nucleotide polymorphisms were extracted from the aligned mitogenomes using SNP-sites (Page et al., 2016).

With the addition of a sequence from E. lotta to be used as an outgroup, we used the aligned mitogenomes to infer a phylogeny with IQTree2 2.2.2.6 (Minh et al., 2020). Mitochondrial sequences were partitioned by feature, with coding regions additionally partitioned by codon position and model selection chosen with the MFP+MERGE setting (Chernomor et al., 2016). Node support was assessed with ultrafast bootstraps and site concordance factors (sCF), both calculated using 1,000 iterations (Hoang et al., 2018Mo et al., 2023).

Results

Microsatellites

The sequencing run of Eai-S43 generated 94,216,390 raw read pairs, of which 93,522,078 read pairs remained after quality and length trimming. De novo assembly resulted in 345,212 contigs with an average length of 395 bp (N50 = 412 bp, min size = 200 bp, max size = 15,385 bp).

QDD identified 7,370 candidate microsatellite sequences from the assembled contigs of Eai-S43, of which 74 were chosen for testing. Nineteen of the 74 loci chosen for testing yielded two or more alleles, indicating they were potentially polymorphic for a larger number of alleles. Thirteen of these loci were optimized for multiplex PCR and genotyping (Table S2). DNA from the whole desiccated larvae provided the highest quality and quantity of DNA, followed by meconium preserved on FTA cards. Frass was found to be an unsuitable template due to poor DNA quality and quantity. Overall, 251 samples that were genotyped were distributed unevenly across years, with 14 individuals from 2016, 56 from 2017, 37 from 2018, and 44 from 2019 (Table S1).

The number of alleles per locus ranged from two to nine, but was low overall, with most loci only having two alleles (Table S5). The locus with the highest number of alleles was Locus-24 with seven. Testing for conformance to HWE indicated that 7 of the 13 loci in the 2018 cohort were out of HWE after Bonferroni correction, and observed heterozygosity (Ho) was lower than expected for all individuals. The average Ho across loci was 0.41 for the 2016 samples, 0.42 for the 2017 samples, and 0.40 for the 2018 samples indicating a relatively low level of heterozygosity in all three populations.

Polymorphic information content (PIC) ranged from 0.06-0.67, P(ID) from 0.13-0.87, P(ID)sib from 0.43-0.93 (Table S6). Locus L-71 scored particularly poorly across all metrics, having the lowest PIC and highest P(ID)/P(ID)sib and was therefore removed from subsequent analyses. After removal of L-71, Ho increased for 2016 and 2018 cohorts, and decreased for 2017 and 2019, whereas fixation index (FIS) estimates increased for 2017 and 2018 (Table 2).

Table 2

Number of alleles (A), effective number of alleles (Ae) observed heterozygosity (Ho), expected heterozygosity, (He) and fixation index (FIS) for each year using 12 microsatellite loci (after L-71 was removed due to low information content).

Using COLONY, individuals were assigned to 20 full sibship groups in the 2017 cohort, 28 in 2018, and 15 in 2019 (Table S7). Across all years, there were 26 full sibship groups with ≥0.9 inclusive probability (total membership = 78). The fixation index (FIS) per cohort was estimated at 0.09-0.17 and Ne was 31-47 (Table 3). Estimated false allele and allele dropout rates varied widely (0-50%) by year and marker (Figure S1), averaging 6.6% and 5%, respectively.

Table 3

Estimates from COLONY of inbreeding (i.e., the fixation index, FIS) and effective population size (Ne) and 95% confidence interval (in parentheses) for three sibling cohorts. Ne was calculated either assuming random or non-random mating.

Inbreeding coefficients varied widely across individuals (range of F = 0.006-0.903; Figure S2), and there was significant variance in inbreeding among individuals (g2 = 0.06 ± 0.018, p = 0.0001).

Mean pairwise relatedness between individuals was similar across years (r = 0.13-0.15) but there were more highly related pairs (r > 0.5; roughly full sibling or higher) in 2018 than other years (25% vs 2- 5%), which could have impacted within-year estimates of Ho and FIS, as well as clustering patterns in population structure analyses. Overall, individuals collected in the same year had significantly higher r values compared to those collected from different years (t = 23.6, df = 23,005, p = 4.7 x 10-122). There was no significant difference in relatedness between individuals from the same or different host plants (t = 1.39, df = 13,917, p = 0.164).

The first two axes of the PCA captured 27.5% of the variance. No combination of axes showed discernable clustering by year or host plant (Figure 1). Structure analyses showed no clear clustering by host plant at any values of K (Figure 2 and S3). Individuals in three years (2016, 2017, and 2019) showed equal proportions of ancestry assigned to each cluster, while 2018 showed some variability in cluster assignment among individuals, possibly due to the higher relatedness among individuals from that year class.

Fig. 1

Principal Component Analysis (PCA) using non-related individuals from all years, colored by larval host plant (BR = Brassica, LE = Lepidium, SI = Sisymbrium, UN = larval host plant not recorded). Percent variance explained by each axis shown in parentheses. There was no clear signal of structure by year or host plant.

Fig. 2

Comparison of estimated ancestry from Structure for population clusters (K) = 3. Each year was run through Structure separately using a reduced-relatedness data set created by removing all but one individual from each of the full sibship clusters (≥0.9 inclusive probability) output from COLONY. Individuals are grouped by larval host plant genus within each year and colors indicate estimated proportion of estimated ancestry to each cluster.

Mitogenomes

We assembled nearly complete circular mitogenomes for all mitogenomic samples except for Eai-162 and EaC8, which only had partial sequences for several genes and were therefore discarded. Mean read depth along the mitogenomes was 63-49,449 reads for shotgun sequencing and 10,833- 23,195 reads for LR PCR (Table S8). The only areas of the mitogenome that showed variation among insulanus assemblies were indels found in the AT-rich control region. However, there were also spikes in read depth in this region which exceeded the coverage across the rest of the mitogenome by >10x, indicating that mapping quality in the control region was likely poor due to repeats. We trimmed 2 bp indels at positions 13,426-13,427 and 13,171-13,172 (relative to the start of the Cytochrome c oxidase subunit 1 gene) to remove possible erroneous assemblies but the control region may still contain errors and/or be incomplete. All insulanus mitogenomes were 15,193 bp long, contained the expected 13 protein-coding genes, 22 tRNAs, 2 rRNAs, D-loop, and replication origin (OH), and had identical sequences after trimming (Figure S4).

We found nine SNPs that are putatively diagnostic for extant insulanus, three of which were not found in the museum specimens (Table 4). One SNP was unique to museum insulanus but not found in the extant samples, whereas three were variable in museum specimens but are now fixed in extant insulanus, potentially representing a loss of historical diversity. There was high node support (bootstrap [B] = 98%, site concordance factor [sCF] = 98%) for extant and museum insulanus forming a monophyletic clade (Figure 3). There was similarly high support (B = 98%, sCF = 98%) for the insulanus clade being placed as sister to E. a. mayi, rather than the geographically adjacent mainland subspecies E. a. transmontana. However, we note that relationships among subspecies were generally poorly resolved, which may reflect a lack of phylogenetic information in the mitogenomes and/or a mismatch between genetic divergence and current taxonomy.

Fig. 3

Mitogenome phylogeny for Euchloe ausinodes subspecies, including insulanus, using extant samples from San Juan Island and museum specimens from Vancouver Island. Diamonds represent node support with site concordance factor (sCF) ≥ 0.9. All insulanus specimens (extant and museum) formed a monophyletic clade with high node support (bootstrap = 98%, sCF = 98%), placed sister to E. a. mayi. The mitogenome for insulanus is available under GenBank accession PQ287242.1. At the time of publication, additional sequences used in the mitogenome phylogeny from the University of Texas Southwestern Medical Center were not available for publication.

Table 4

Mitochondrial single nucleotide polymorphisms (SNPs) found in Euchloe ausonides insulanus compared to other E. ausonides subspecies. Positions are relative to start of COX1 gene. All SNPs listed are potentially diagnostic for extant E. a. insulanus, except for position 6161, which was only found in museum insulanus specimens and is therefore likely to have been lost. Reference alleles are those found in other E. ausonides subspecies included in this study. * Potentially unique to insulanus on San Juan Island; ** Potentially diagnostic for extant insulanus population; † Allele likely lost in extant insulanus population. Mitogenome for insulanus available under GenBank accession PQ287242.1. At the time of publication, additional sequences used in the mitogenome phylogeny from the University of Texas Southwestern Medical Center were not available for publication.

Discussion

This study offers a first look into the genetic diversity of the endangered butterfly species Euchloe ausonides insulanus. We generated microsatellite markers from DNA collected non-invasively from exuviae, meconium, and other sources, showing the potential of monitoring genetic diversity in this species without relying on lethal sampling techniques or incidental mortalities. Additionally, the annotated mitochondrial genome presented here for insulanus sheds light on the genetic diversity within the species and its relationship to other members of the E. ausonides species complex. Our work also provides putatively diagnosable genetic markers to distinguish insulanus from other subspecies in the complex.

Our results indicate that insulanus has relatively low heterozygosity, a small effective population size (Ne), and low allelic diversity. These factors limit the available genetic variation the species has to draw from and reduce the likelihood that advantageous mutations will arise in the future, both of which decrease the species’ adaptive potential (Hoffmann et al., 2017Rousselle et al., 2020). The number of alleles per locus and observed heterozygosity estimated with microsatellites was lower than those found in most previous butterfly studies (as reviewed in Heffernan et al. 2024), and on par with some declining species, such as the endangered Poweshiek skipperling (Oarisma poweshiek) (Saarinen et al., 2016).

Furthermore, an Ne of 50-500 may be necessary to avoid extinction (Frankham et al., 2014; Pérez- Pereira et al., 2022), and our estimates for insulanus were well below that minimum. The lack of population structure and low genetic diversity also indicate that insulanus likely exists as a single population isolated on one island with no connectivity to other subspecies that could provide additional genetic diversity through gene flow. These circumstances may restrict insulanus’ access to the genetic diversity needed to assemble novel adaptive genotypes in the face of climate change and other stressors (Habel and Schmitt, 2009de Jong et al., 2023). Hence, insulanus may benefit from management approaches that focus on increasing genetic diversity through genetic rescue and expanding the number of populations through translocation of captive-reared individuals (refer to Considerations for species management section below).

Additionally, some individuals had high inbreeding coefficients (Figure S2), indicating there is a risk of inbreeding depression for the species if unconstrained inbreeding continues. Inbreeding depression in butterflies has been shown to reduce thermal tolerance (Dierks et al., 2012a), decrease fertility and larval survival (van Oosterhout et al., 2000; Cassel et al., 2001Nieminen et al., 2001), decrease adult lifespan (Saccheri et al., 1998De Nardin et al., 2016), reduce adaptive potential (Dierks et al., 2012b), and ultimately increase the risk of extinction (Saccheri et al., 1998Nonaka et al., 2019). The negative consequences of inbreeding depression largely result from the unmasking of deleterious alleles, when inbreeding leads to a higher number of individuals expressing recessive alleles with negative fitness consequences (Kyriazis et al., 2021), and additional studies to quantify inbreeding depression in insulanus may prove helpful to the conservation of the species. Furthermore, since inbreeding coefficients varied greatly among individuals and pairwise relatedness indicated that most individuals collected for this study were not closely related, careful selection of individuals during reintroduction and/or translocation could help to mitigate the effects of inbreeding.

Genetic diversity is not the only measure to consider when it comes to survival of the species. For example, low heterozygosity and low population size in the Uncompahgre fritillary (Boloria acrocnema) led Britten et al. (1994) to predict the species’ extinction as “inevitable.” Instead, the species persisted and maintained genetic diversity over time, rather than succumbing to the effects of genetic erosion (Monroe et al., 2016). Although low genetic diversity is a serious risk factor for extinction, it must be considered within the context of a species’ ecology and demographic history to determine the likely impact (Schultz et al., 2019Teixeira and Huber, 2021DiLeo et al., 2024). The next section provides more detail on how these results could be integrated into management considerations to improve the long-term outlook of the species.

Considerations for species management

Genetic rescue. Due to its small effective population size, low levels of genetic variation, and evidence of inbreeding, insulanus may be a good candidate for genetic rescue (Fitzpatrick et al., 2023), whereby new alleles are introduced into a population (i.e., gene flow) to ameliorate the deleterious effects of low genetic diversity and decrease the likelihood of extinction (Bell et al., 2019Hohenlohe et al., 2021). Increasing gene flow by improving connectivity between populations has been successful in slowing the erosion of genetic diversity in other butterfly species, such as the Rocky Mountain parnassian (Parnassius smintheus) (Jangjoo et al., 2016), but this strategy is impractical for insulanus, which is island-bound and only has one population. Instead, limited hybridization of insulanus with one of the closely related subspecies could be an option for increasing genetic diversity. This approach was successfully used with an isolated population of the Marsh Fritillary butterfly (Euphydryas aurinia), which benefited from an influx of new alleles but still maintained genetic distinctiveness (Davis et al., 2021). In addition to increasing genetic diversity, genetic rescue can increase population growth rates, which may be necessary to provide enough surplus individuals to translocate into new populations (Bell et al., 2019Berger-Tal et al., 2020). However, whereas there was only weak divergence between the mitogenomes of insulanus and some of the other subspecies (i.e., ausonides and mayi), a thorough analysis of the differences in nuclear genetic divergence between insulanus and the mainland subspecies would provide further information on whether genetic rescue through selective hybridization is feasible and identify potential genetic incompatibilities to reduce the risk of outbreeding depression.

Establishment of new populations. Currently, insulanus is only found in one small population, leaving it acutely vulnerable to stochastic events. We also identified multiple highly related (i.e., likely full sibling or higher) and inbred individuals which, when reintroduced into the same location after captive rearing, may contribute to an increase in inbreeding depression. One option to consider for mitigating these issues is to translocate highly related individuals to different locations where they could form new populations. Translocation, or movement of individuals from one location to another, can be an important tool for threatened and endangered species (Schultz et al., 2008Daniels et al., 2018). For example, the large blue (Maculinea arion), scarce large blue (Phengaris teleius), chequered skipper (Carterocephalus palaemon), and clouded Apollo (Parnassius mnemosyne) were all successfully translocated and able to establish new populations that persisted >5 years (Andersen et al., 2014Kuussaari et al., 2015Bourn et al., 2024Sánchez-García et al., 2024). The few studies that have tracked genetic diversity in translocated populations have been mixed; whereas M. arion showed evidence of new allelic diversity, P. teleius lost allelic diversity over time, likely due to founder effects (Andersen et al., 2014Sánchez-García et al., 2024). Founder effects and existing inbreeding load also likely led to deformed wings in a translocated population of the Apollo butterfly (Parnassius apollo) (Witkowski et al., 1997Adamski and Witkowski, 1999). Consequently, insulanus may benefit from additional research to quantify inbreeding load to reduce the likelihood of deleterious phenotypes arising in the new populations created through translocation.

Introduction of novel larval host plants. Larval host plant specificity has been identified as a major extinction risk for butterflies (Koh et al., 2004Palash et al., 2022). Guppy and Shepard (2001) suggested that hairy rock cress (Arabis eschscholtziana, previously A. hirsuta) may have been the original larval host plant for insulanus. In addition to a native pepperweed (Lepidium virginicum var. menziesii), insulanus included in this study were collected from the recently introduced tumble mustard (Sisymbrium altissimum) and field mustard (Brassica rapa). Lambert (2011) had some success using tower mustard (Turritis glabra) as a larval host plant for insulanus, providing further evidence of plasticity in host plant use in this species. Other research documenting host plant use by Euchloe ausonides subspecies show that they use a variety of weedy mustards (Brassicaceae) (Opler, 1974Guppy and Shepard, 2001Back et al., 2011), which could provide additional options for oviposition and larval food sources if made available. This flexibility in host plant use is underscored by our results, which show a lack of population structure and relatedness by host plant, suggesting that insulanus females do not exhibit strong natal fidelity when choosing ovipositing sites (Nice et al., 2002Downey and Nice, 2011).

Since adults lay eggs on unopened flower buds, the species may be susceptible to shifts in flowering phenology due to climate change, which could lead to asynchrony between flower budding and ovipositing (Hill et al., 2021). Maintaining access to a range of hosts may increase the likelihood that flower budding in at least one host species will continue to overlap with the timing of egg deposition (Navarro-Cano et al., 2015). Providing access to additional host plants could therefore prove useful in reducing the vulnerability of insulanus to environmental stressors. However, recently introduced non-native host plants can also act as ecological traps when novel phytochemicals negatively impact larval development or lead to phenotypic changes with knock-on effects to reproduction (Braga, 2023Jarrett and Miller, 2024). Thus, additional studies on the suitability of novel larval host plants for insulanus are warranted prior to their introduction to existing habitat or use as larval food during captive rearing.

Considerations for ongoing genetic monitoring

Marker choice. Molecular methods have improved substantially since this study began in 2014, and the adoption of more advanced techniques may be useful for ongoing genetic monitoring. For example, the potential issues of using capillary electrophoresis for genotyping microsatellites (i.e., difficulty transferring between laboratories, subjectivity in scoring, and labor intensiveness have been largely resolved by newer sequence-based microsatellite genotyping methods (such as SSRseq, also known as genotyping-by-sequencing [GBS]) (De Barba et al., 2017). Whereas capillary electrophoresis relies on visually scoring alleles by sequence length, SSRseq involves sequencing the entire microsatellite (Darby et al., 2016Vartia et al., 2016). Alleles can then be called bioinformatically to produce typical genotypes (Barbian et al., 2018Tibihika et al., 2019) or the entire sequence can be used to generate high-resolution haplotypes (Lepais et al., 2020), reducing the risk of false alleles being called. SSRseq also increases genotyping success with degraded DNA compared to capillary electrophoresis (De Barba et al., 2017), hence potentially resolving the issues with allele dropout seen with some loci in this study. Nevertheless, not all pre-existing microsatellites perform well with SSRseq, so implementing this method could require additional microsatellites to be developed from the candidate sequences generated for this study.

Source of genetic material. Our results indicate that meconium collected on FTA cards, larval exuviae, and pupae from natural mortalities all generally produced high enough quality DNA to be used for microsatellite genotyping, while frass did not. However, there were issues with allele dropout for the microsatellite data, which needed to be accounted for when assessing the results. Limiting genetic monitoring using opportunistically collected samples and natural mortalities could also bias the results, since most population genetic analyses assume that samples are randomly selected and representative of the population (Phillips et al., 2019). These biases may be reflected in the Structure results, which showed more clustering among individuals collected in 2018, possibly due to higher relatedness among individuals from non-random sampling skewing results (Schwartz and McKelvey, 2009Puechmaille, 2016).

To mitigate these issues, it could be advantageous to collect fresh tissue samples from the general population, rather than incidental samples from larvae in the captive rearing program. Sample design that limits the collection of too many individuals from one area, coupled with methods such as non-lethal wing clipping and/or leg removal, could provide higher quality DNA and reduce bias in the data. Wing clipping and leg removal have been shown to have no significant impact on behavior or survival in multiple species of butterfly (Hamm et al., 2010Koscinski et al., 2011Crawford et al., 2013). Wing clipping was previously approved by the USFWS for use in at least two endangered butterfly species (i.e., Mitchell’s satyr, Neonympha mitchellii mitchellii, and Karner blue, Plebejus samuelis) (Hamm et al., 2014Zhang et al., 2024), and thus may be an effective technique for obtaining higher quality DNA from adults that could be used to improve the accuracy of microsatellite-based monitoring or provide the material needed to more fine-scale genomic analyses (Webster et al., 2023).

Conclusion

Euchloe ausonides insulanus has lost most of its habitat over the last ∼100 years, leading to a marked reduction in population size. The results on this study underscore those demographic changes, providing genetic evidence that there is likely a single population remaining, with low genetic diversity and a small effective population size. However, there is often a temporal lag between population declines and genetic erosion, leading to an overestimate of genetic diversity measures immediately following demographic changes (Gargiulo et al., 2025). Hence, the full effects of the demographic events that have been affecting may not yet be reflected in the genetic diversity shown in this study and further declines in genetic diversity may continue to occur even if the population remains stable going forward. Given these concerns, conservation strategies such as genetic rescue and managed translocation may be necessary to mitigate inbreeding effects and bolster genetic diversity. The success of such interventions in other butterfly species suggests that similar approaches could increase adaptive potential and improve the long-term viability of insulanus in the face of ongoing environmental and demographic challenges.

Supporting information

Supplemental Tables and Figures[supplements/652852_file03.docx]

Statements and Declarations

Funding

This work was supported by the United States Geological Survey Science Support Program in FY2017-2020 and FY2021-FY2024.

Competing interests: The authors have no relevant financial or non-financial interests to disclose.

Author contributions

All authors contributed to the study conception and design. Material preparation, data collection, and lab work was performed by Aaron Aunins, Colleen Young, and Robin Johnson.

Data analyses were performed and interpreted by Kara Jones, Aaron Aunins, and Cheryl Morrison. The first draft of the manuscript was written by Kara Jones, Aaron Aunins, and Cheryl Morrison. All authors approved the final manuscript.

Data availability

Sequencing data and associated sample metadata were deposited to the National Center for Biotechnology Information (NCBI) under BioProject PRJNA1156227 (https://doi.org/10.5066/P136HXY4). The Euchloe ausonides insulanus annotated mitochondrial genome was deposited on GenBank under accession PQ287242. Microsatellite genotypes, ddRAD single nucleotide polymorphisms (SNPs), and additional sample metadata were deposited on Zenodo (https://doi.org/10.5281/zenodo.14502991). At the time of publication, additional sequences used in the mitogenome phylogeny from the University of Texas Southwestern Medical Center not been made available for publication.

Acknowledgments

We thank the following individuals for their contributions to this study: Jing Zhang, Qian Cong, and Nick Grishin (University of Texas Southwestern Medical Center) for providing mitochondrial sequence data and specimen metadata; Michael Eackles (USGS) for assisting with lab work, designing the multiplex PCRs, and microsatellite genotyping; Jenny Shrum (USFWS) for providing samples and metadata; Karen Reagan (USGS) for initial conceptualization of the study and procuring samples; Erin Sullivan (Woodland Park Zoo) for providing samples of additional E. ausonides subspecies; San Juan Island National Historical Park (NPS) staff for assisting with sample collection; Jaret Daniels (McGuire Center for Lepidoptera and Biodiversity, Florida Museum of Natural History) for assisting with initial study planning; Julie Combs (Washington Department of Fish and Wildlife) for additional information on the natural and management history of insulanus; and Erin Adams (USFWS) and two USGS peer reviewers for their helpful comments and suggestions on a draft of this manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Funder Information Declared

United States Fish and Wildlife Service, Science Support Program

Footnotes

References

  1. Adamski P, Witkowski Z (1999) Wing deformation in an isolated Carpathian population of Parnassius apollo. Nota lepidoptera 22: 67–73.
  2. Afgan E, Baker D, Batut B, van den Beek M, Bouvier D, Čech M, Chilton J, Clements D, Coraor N, Grüning BA, et al. (2018) The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Research 46: W537–W544.
  3. Al Arab M, Hönerzu Siederdissen C, Tout K, Sahyoun AH, Stadler PF, Bernt M (2017) Accurate annotation of protein-coding genes in mitochondrial genomes. Molecular Phylogenetics and Evolution 106: 209–216.
  4. Allio R, Schomaker-Bastos A, Romiguier J, Prosdocimi F, Nabholz B, Delsuc F (2020) MitoFinder: Efficient automated large-scale extraction of mitogenomic data in target enrichment phylogenomics. Molecular Ecology Resources 20: 892–905.
  5. Andersen A, Simcox DJ, Thomas JA, Nash DR (2014) Assessing reintroduction schemes by comparing genetic diversity of reintroduced and source populations: A case study of the globally threatened large blue butterfly (Maculinea arion). Biological Conservation 175: 34–41.
  6. Applied Biosystems (2009) GeneMapper.
  7. Arantes LS, Caccavo JA, Sullivan JK, Sparmann S, Mbedi S, Höner OP, Mazzoni CJ (2023) Scaling-up RADseq methods for large datasets of non-invasive samples: Lessons for library construction and data preprocessing. Molecular Ecology Resources 1755–0998.13859.
  8. Back W, Miller MA, Opler PA (2011) Genetic, phenetic, and distributional relationships of Nearctic Euchloe (Pieridae, Pierinae, Anthocharidini). lepi 65: 1–14.
  9. Balmori-de la Puente A, Escoda L, Fernández-González Á, Menéndez-Pérez D, González-Esteban J, Castresana J (2023) Evaluating the use of non-invasive hair sampling and ddRAD to characterize populations of endangered species: Application to a peripheral population of the European mink. Ecology and Evolution 13: e10530.
  10. Barbian HJ, Connell AJ, Avitto AN, Russell RM, Smith AG, Gundlapally MS, Shazad AL, Li Y, Bibollet-Ruche F, Wroblewski EE, et al. (2018) CHIIMP: An automated high-throughput microsatellite genotyping platform reveals greater allelic diversity in wild chimpanzees. Ecology and Evolution 8: 7946–7963.
  11. Bell DA, Robinson ZL, Funk WC, Fitzpatrick SW, Allendorf FW, Tallmon DA, Whiteley AR (2019) The exciting potential and remaining uncertainties of genetic rescue. Trends in Ecology & Evolution 34: 1070–1079.
  12. Benson DA, Cavanaugh M, Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW (2017) GenBank. Nucleic Acids Res 45: D37–D42.
  13. Berger-Tal O, Blumstein DT, Swaisgood RR (2020) Conservation translocations: a review of common difficulties and promising directions. Animal Conservation 23: 121–131.
  14. Bourn NAD, O’Riordan S, Maes D, Goffart P, Shadbolt T, Hordley L, Sainsbury AW, Bulman C, Hoare D, Field R, et al. (2024) The history, science and preliminary results from the reintroduction of the Chequered Skipper, Carterocephalus palaemon into Rockingham Forest, England. J Insect Conserv. doi:10.1007/s10841-024-00601-3
  15. Braga MP (2023) Are exotic host plants a life raft or a trap for butterflies? Current Opinion in Insect Science 58: 101074.
  16. Britten HB, Brussard PF, Murphy DD (1994) The pending extinction of the Uncompahgre fritillary butterfly. Conservation Biology 8: 86–94.
  17. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL (2009) BLAST+: architecture and applications. BMC Bioinformatics 10: 421.
  18. Cassel A, Windig J, Nylin S, Wiklund C (2001) Effects of population size and food stress on fitness- related characters in the scarce heath, a rare butterfly in Western Europe. Conservation Biology 15: 1667–1673.
  19. Chan PP, Lin BY, Mak AJ, Lowe TM (2021) tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes. Nucleic Acids Research 49: 9077–9096.
  20. Chernomor O, von Haeseler A, Minh BQ (2016) Terrace aware data structure for phylogenomic inference from supermatrices. Systematic Biology 65: 997–1008.
  21. Couch SP, Bray AP, Ismay C, Chasnovski E, Baumer BS, Çetinkaya-Rundel M (2021) infer: An R package for tidyverse-friendly statistical inference. Journal of Open Source Software 6: 3661.
  22. Crawford LA, Koscinski D, Watt KM, McNeil JN, Keyghobadi N (2013) Mating success and oviposition of a butterfly are not affected by non-lethal tissue sampling. J Insect Conserv 17: 859–864.
  23. Daniels JC, Nordmeyer C, Runquist E (2018) Improving standards for at-risk butterfly translocations. Diversity 10: 67.
  24. Darby BJ, Erickson SF, Hervey SD, Ellis-Felege SN (2016) Digital fragment analysis of short tandem repeats by high-throughput amplicon sequencing. Ecology and Evolution 6: 4502–4512.
  25. Davis ML, Barker C, Powell I, Porter K, Ashton P (2021) Combining modelling, field data and genetic variation to understand the post-reintroduction population genetics of the Marsh Fritillary butterfly (Euphydryas aurinia). J Insect Conserv 25: 875–886.
  26. De Barba M, Miquel C, Lobréaux S, Quenette PY, Swenson JE, Taberlet P (2017) High-throughput microsatellite genotyping in ecology: improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Molecular Ecology Resources 17: 492–507.
  27. de Jong M, van Rensburg AJ, Whiteford S, Yung CJ, Beaumont M, Jiggins C, Bridle J (2023) Rapid evolution of novel biotic interactions in the UK Brown Argus butterfly uses genomic variation from across its geographical range. Molecular Ecology 32: 5742–5756.
  28. De Nardin J, Da Silva L, De Araújo AM (2016) Comparative effects of inbreeding and outbreeding on the ontogeny of Heliconius erato phyllis (Lepidoptera: Nymphalidae). Entomological Science 19: 304–309.
  29. Dierks A, Baumann B, Fischer K (2012a) Response to selection on cold tolerance is constrained by inbreeding. Evolution 66: 2384–2398.
  30. Dierks A, Hoffmann B, Bauerfeind SS, Fischer K (2012b) Effects of inbreeding on life history and thermal performance in the tropical butterfly Bicyclus anynana. Popul Ecol 54: 83–90.
  31. DiLeo MF, Nair A, Kardos M, Husby A, Saastamoinen M (2024) Demography and environment modulate the effects of genetic diversity on extinction risk in a butterfly metapopulation. Proc Natl Acad Sci USA 121: e2309455121.
  32. Donath A, Jühling F, Al-Arab M, Bernhart SH, Reinhardt F, Stadler PF, Middendorf M, Bernt M (2019) Improved annotation of protein-coding genes boundaries in metazoan mitochondrial genomes. Nucleic Acids Research 47: 10543–10552.
  33. Downey MH, Nice CC (2011) Experimental evidence of host race formation in Mitoura butterflies (Lepidoptera: Lycaenidae). Oikos 120: 1165–1174.
  34. Ebdon S, Bisschop G, Lohse K, Saccheri I, Davies J, Wellcome Sanger Institute Tree of Life programme, Wellcome Sanger Institute Scientific Operations: DNA Pipelines collective, Tree of Life Core Informatics collective, Darwin Tree of Life Consortium (2022) The genome sequence of the orange-tip butterfly, Anthocharis cardamines (Linnaeus, 1758). Wellcome Open Res 7: 260.
  35. Edgar RC (2004) MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5: 113.
  36. Falush D, Stephens M, Pritchard JK (2007) Inference of population structure using multilocus genotype data: dominant markers and null alleles. Molecular Ecology Notes 7: 574–578.
  37. Fernandes F, Pereira L, Freitas AT (2009) CSA: An efficient algorithm to improve circular DNA multiple alignment. BMC Bioinformatics 10: 230.
  38. Fitzpatrick SW, Mittan-Moreau C, Miller M, Judson JM (2023) Genetic rescue remains underused for aiding recovery of federally listed vertebrates in the United States. J Hered 114: 354–366.
  39. Frankham R, Bradshaw CJA, Brook BW (2014) Genetics in conservation management: Revised recommendations for the 50/500 rules, Red List criteria and population viability analyses. Biological Conservation 170: 56–63.
  40. Gargiulo R, Budde KB, Heuertz M (2025) Mind the lag: understanding genetic extinction debt for conservation. Trends in Ecology & Evolution 40: 228–237.
  41. Guppy C, Shepard J (2001) Butterflies of British Columbia. UBC Press, British Columbia, Canada.
  42. Habel JC, Schmitt T (2009) The genetic consequences of different dispersal behaviours in Lycaenid butterfly species. Bulletin of Entomological Research 99: 513–523.
  43. Hamm CA, Aggarwal D, Landis DA (2010) Evaluating the impact of non-lethal DNA sampling on two butterflies, Vanessa cardui and Satyrodes eurydice. J Insect Conserv 14: 11–18.
  44. Hamm CA, Rademacher V, Landis DA, Williams BL (2014) Conservation Genetics and the Implication for Recovery of the Endangered Mitchell’s Satyr Butterfly, Neonympha mitchellii mitchellii. Journal of Heredity 105: 19–27.
  45. Heffernan E, Barkdull M, Brady N (2024) Microsatellites for butterfly conservation: historical challenges, current relevance, and a guide to implementation. Front Ecol Evol 12. doi:10.3389/fevo.2024.1344065
  46. Hill GM, Kawahara AY, Daniels JC, Bateman CC, Scheffers BR (2021) Climate change effects on animal ecology: butterflies and moths as a case study. Biological Reviews 96: 2113–2126.
  47. Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS (2018) UFBoot2: Improving the ultrafast bootstrap approximation. Molecular Biology and Evolution 35: 518–522.
  48. Hoffmann AA, Sgrò CM, Kristensen TN (2017) Revisiting adaptive potential, population size, and conservation. Trends in Ecology & Evolution 32: 506–517.
  49. Hohenlohe PA, Funk WC, Rajora OP (2021) Population genomics for wildlife conservation and management. Molecular Ecology 30: 62–82.
  50. Jangjoo M, Matter SF, Roland J, Keyghobadi N (2016) Connectivity rescues genetic diversity after a demographic bottleneck in a butterfly population network. Proceedings of the National Academy of Sciences 113: 10914–10919.
  51. Jarrett BJM, Miller CW (2024) Host plant effects on sexual selection dynamics in phytophagous insects. Annu Rev Entomol 69: 41–57.
  52. Jombart T (2008) adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24: 1403–1405.
  53. Jones OR, Wang J (2010) COLONY: a program for parentage and sibship inference from multilocus genotype data. Molecular Ecology Resources 10: 551–555.
  54. Koh LP, Sodhi NS, Brook BW (2004) Ecological Correlates of Extinction Proneness in Tropical Butterflies. Conservation Biology 18: 1571–1578.
  55. Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I (2015) Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Molecular Ecology Resources 15: 1179–1191.
  56. Koscinski D, Crawford LA, Keller HA, Keyghobadi N (2011) Effects of different methods of non-lethal tissue sampling on butterflies. Ecological Entomology 36: 301–308.
  57. Kuussaari M, Heikkinen RK, Heliölä J, Luoto M, Mayer M, Rytteri S, von Bagh P (2015) Successful translocation of the threatened Clouded Apollo butterfly (Parnassius mnemosyne) and metapopulation establishment in southern Finland. Biological Conservation 190: 51–59.
  58. Kyriazis CC, Wayne RK, Lohmueller KE (2021) Strongly deleterious mutations are a primary determinant of extinction risk due to inbreeding depression. Evolution Letters 5: 33–47.
  59. Lambert AM (2011) Natural History and Population Ecology of a Rare Pierid Butterfly, Euchloe Ausonides Insulanus Guppy and Shepard (Pieridae) (Dissertation). University of Washington, Seattle, WA.
  60. Larsson A (2014) AliView: a fast and lightweight alignment viewer and editor for large datasets. Bioinformatics 30: 3276–3278.
  61. Lepais O, Chancerel E, Boury C, Salin F, Manicki A, Taillebois L, Dutech C, Aissi A, Bacles CFE, Daverat F, et al. (2020) Fast sequence-based microsatellite genotyping development workflow. PeerJ 8: e9085.
  62. Li D, Liu C-M, Luo R, Sadakane K, Lam T-W (2015) MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31: 1674–1676.
  63. McDermott Long O, Warren R, Price J, Brereton TM, Botham MS, Franco AMA (2017) Sensitivity of UK butterflies to local climatic extremes: which life stages are most at risk? Journal of Animal Ecology 86: 108–116.
  64. Meglécz E, Costedoat C, Dubut V, Gilles A, Malausa T, Pech N, Martin J-F (2010) QDD: a user- friendly program to select microsatellite markers and design primers from large sequencing projects. Bioinformatics 26: 403–404.
  65. Meng G, Li Y, Yang C, Liu S (2019) MitoZ: a toolkit for animal mitochondrial genome assembly, annotation and visualization. Nucleic Acids Research 47: e63.
  66. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, Lanfear R (2020) IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era. Molecular Biology and Evolution 37: 1530–1534.
  67. Mo YK, Lanfear R, Hahn MW, Minh BQ (2023) Updated site concordance factors minimize effects of homoplasy and taxon sampling. Bioinformatics 39: btac741.
  68. Monroe EM, Alexander KD, Britten HB (2016) Still here after all these years: the persistence of the Uncompahgre fritillary butterfly. J Insect Conserv 20: 305–313.
  69. Navarro-Cano JA, Karlsson B, Posledovich D, Toftegaard T, Wiklund C, Ehrlén J, Gotthard K (2015) Climate change, phenology, and butterfly host plant utilization. AMBIO 44: 78–88.
  70. Nice CC, Fordyce JA, Shapiro AM, Ffrench-Constant R (2002) Lack of evidence for reproductive isolation among ecologically specialised lycaenid butterflies. Ecological Entomology 27: 702– 712.
  71. Nieminen M, Singer MC, Fortelius W, Schöps K, Hanski I (2001) Experimental confirmation that inbreeding depression increases extinction risk in butterfly populations. The American Naturalist 157: 237–244.
  72. Nonaka E, Sirén J, Somervuo P, Ruokolainen L, Ovaskainen O, Hanski I (2019) Scaling up the effects of inbreeding depression from individuals to metapopulations. Journal of Animal Ecology 88: 1202–1214.
  73. Opler P (1974) Studies on nearctic Euchloe. Part 7. Comparative life histories, hosts, and the morphology of immature stages. Journal of Research on the Lepidoptera 13: 1–20.
  74. Page AJ, Taylor B, Delaney AJ, Soares J, Seemann T, Keane JA, Harris SR (2016) SNP-sites: rapid efficient extraction of SNPs from multi-FASTA alignments. Microbial Genomics 2: e000056.
  75. Palash A, Paul S, Resha SK, Khan MK (2022) Body size and diet breadth drive local extinction risk in butterflies. Heliyon 8. doi:10.1016/j.heliyon.2022.e10290
  76. Park E, Hwang D-S, Lee J-S, Song J-I, Seo T-K, Won Y-J (2012) Estimation of divergence times in cnidarian evolution based on mitochondrial protein-coding genes and the fossil record. Molecular Phylogenetics and Evolution 62: 329–345.
  77. Patterson N, Price AL, Reich D (2006) Population structure and eigenanalysis. PLOS Genetics 2: e190.
  78. Peakall R, Smouse PE (2012) GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28: 2537–2539.
  79. Pérez-Pereira N, Wang J, Quesada H, Caballero A (2022) Prediction of the minimum effective size of a population viable in the long term. Biodivers Conserv 31: 2763–2780.
  80. Phillips JD, Gillis DJ, Hanner RH (2019) Incomplete estimates of genetic diversity within species: Implications for DNA barcoding. Ecology and Evolution 9: 2996–3010.
  81. Ponce M, Micol J (1992) PCR amplification of long DNA fragments. Nucleic Acids Research 20: 623–623.
  82. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155: 945–959.
  83. Puechmaille SJ (2016) The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Molecular Ecology Resources 16: 608–627.
  84. Qiagen (2016) CLC Genomics Workbench.
  85. R Core Team (2024) R: A Language and Environment for Statistical Computing.
  86. Rousselle M, Simion P, Tilak M-K, Figuet E, Nabholz B, Galtier N (2020) Is adaptation limited by mutation? A timescale-dependent effect of genetic diversity on the adaptive substitution rate in animals. PLOS Genetics 16: e1008668.
  87. Saarinen EV, Reilly PF, Austin JD (2016) Conservation genetics of an endangered grassland butterfly (Oarisma poweshiek) reveals historically high gene flow despite recent and rapid range loss. Insect Conserv Diversity 9: 517–528.
  88. Saccheri I, Kuussaari M, Kankare M, Vikman P, Fortelius W, Hanski I (1998) Inbreeding and extinction in a butterfly metapopulation. Nature 392: 491–494.
  89. Sánchez-García D, Wynhoff I, Kajzer-Bonk J, Sztencel-Jabłonka A, Nowicki P, Casacci LP, Witek M (2024) Temporal and spatial variation of morphological traits and genetic structure in Phengaris teleius myrmecophilous butterflies following habitat and climate changes three decades after reintroduction. Global Ecology and Conservation 54: e03104.
  90. Schultz AJ, Strickland K, Cristescu RH, Hanger J, de Villiers D, Frère CH (2022) Testing the effectiveness of genetic monitoring using genetic non-invasive sampling. Ecology and Evolution 12: e8459.
  91. Schultz CB, Haddad NM, Henry EH, Crone EE (2019) Movement and demography of at-risk butterflies: Building blocks for conservation. Annual Review of Entomology 64: 167–184.
  92. Schultz CB, Russell C, Wynn L (2008) Restoration, reintroduction, and captive propagation for at-risk butterflies: A review of British and American conservation efforts. Israel Journal of Ecology & Evolution 54: 41–61.
  93. Schwartz MK, McKelvey KS (2009) Why sampling scheme matters: the effect of sampling scheme on landscape genetic results. Conserv Genet 10: 441–452.
  94. Serrote CML, Reiniger LRS, Silva KB, Rabaiolli SM dos S, Stefanel CM (2020) Determining the Polymorphism Information Content of a molecular marker. Gene 726: 144175.
  95. Shepard JH (2000) Status of Five Butterflies and Skippers in British Columbia (No. Wildlife Working Report No. WR-101). British Columbia Ministry of Environment, Lands and Parks Wildlife Branch and Resources Inventory Branch, Victoria, BC.
  96. Sittenthaler M, Fischer I, Chovanec A, Koblmüller S, Macek O, Sattmann H, Szucsich N, Zangl L, Haring E (2023) DNA barcoding of exuviae for species identification of Central European damselflies and dragonflies (Insecta: Odonata). J Insect Conserv 27: 435–450.
  97. Stoffel MA, Esser M, Kardos M, Humble E, Nichols H, David P, Hoffman JI (2016) inbreedR: an R package for the analysis of inbreeding based on genetic markers. Methods Ecol Evol 7: 1331– 1339.
  98. Teixeira JC, Huber CD (2021) The inflated significance of neutral genetic diversity in conservation genetics. Proceedings of the National Academy of Sciences 118: e2015096118.
  99. Tibihika PD, Curto M, Dornstauder-Schrammel E, Winter S, Alemayehu E, Waidbacher H, Meimberg H (2019) Application of microsatellite genotyping by sequencing (SSR-GBS) to measure genetic diversity of the East African Oreochromis niloticus. Conserv Genet 20: 357–372.
  100. Tourvas N (2024) PopGenUtils: A collection of useful functions to deal with genetic data in
  101. R. van Oosterhout C, Zulstra WG, van Heuven MK, Brakefield PM (2000) Inbreeding depression and genetic load in laboratory metapopulations of the butterfly Bicyclus Anynana. Evolution 54: 218–225.
  102. Vartia S, Villanueva-Cañas JL, Finarelli J, Farrell ED, Collins PC, Hughes GM, Carlsson JEL, Gauthier DT, McGinnity P, Cross TF, et al. (2016) A novel method of microsatellite genotyping-by- sequencing using individual combinatorial barcoding. Royal Society Open Science 3: 150565.
  103. Waits JL, Leberg PL (2000) Biases associated with population estimation using molecular tagging. Animal Conservation 3: 191–199.
  104. Waits LP, Luikart G, Taberlet P (2001) Estimating the probability of identity among genotypes in natural populations: cautions and guidelines. Molecular Ecology 10: 249–256.
  105. Wang J (2009) A new method for estimating effective population sizes from a single sample of multilocus genotypes. Molecular Ecology 18: 2148–2164.
  106. Wang J (2018) Estimating genotyping errors from genotype and reconstructed pedigree data. Methods Ecol Evol 9: 109–120.
  107. Wang J (2022) A joint likelihood estimator of relatedness and allele frequencies from a small sample of individuals. Methods Ecol Evol 13: 2443–2462.
  108. Webster MT, Beaurepaire A, Neumann P, Stolle E (2023) Population Genomics for Insect Conservation. Annual Review of Animal Biosciences 11: 115–140.
  109. Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, et al. (2019) Welcome to the Tidyverse. Journal of Open Source Software 4: 1686.
  110. Witkowski Z, Kosior A, Adamski P, Plonka (1997) Extinction and reintroduction of Parnassius apollo in the Pieniny National Park (Polish Carpathians). Biologia, Bratislava 52: 199–208.
  111. Zhan X, Zheng X, Bruford MW, Wei F, Tao Y (2010) A new method for quantifying genotyping errors for noninvasive genetic studies. Conserv Genet 11: 1567–1571.
  112. Zhang J, Aunins AW, King TL, Cong Q, Shen J, Song L, Schuurman GW, Knutson RL, Grundel R, Hellmann J, et al. (2024) Range-wide population genomic structure of the Karner blue butterfly, Plebejus (Lycaeides) samuelis. Ecology and Evolution 14: e70044.

https://i.liadm.com/sync-container?duid=28e3293678dc–01k1hv9vygvqagd526m1kfdfv6&ds=did-0049&euns=1&s=CgA&version=v3.13.0&cd=.biorxiv.org&pv=3bf36413-505f-4389-8af6-4315d15efbab