Whole genomes from the extinct Xerces Blue butterfly can help identify declining insect species
eLife assessment
This important study illustrates the value of museum samples for understanding past genetic variability in the genomes of populations and species, including those that no longer exist. The authors present genomic sequencing data for the extinct Xerces Blue butterfly and report convincing evidence of declining population sizes and increases in inbreeding beginning 75,000 years ago, which strongly contrasts to the patterns observed in similar data from its closest relative, the extant Silvery Blue butterfly. Such long-term population health indicators may be used to highlight still extant but especially vulnerable-to-extinction insect species – irrespective of their current census population size abundance.
https://doi.org/10.7554/eLife.87928.3.sa0Important: Findings that have theoretical or practical implications beyond a single subfield
- Landmark
- Fundamental
- Important
- Valuable
- Useful
Convincing: Appropriate and validated methodology in line with current state-of-the-art
- Exceptional
- Compelling
- Convincing
- Solid
- Incomplete
- Inadequate
During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments
Abstract
The Xerces Blue (Glaucopsyche xerces) is considered to be the first butterfly to become extinct in historical times. It was notable for its chalky lavender wings with conspicuous white spots on the ventral wings. The last individuals were collected in their restricted habitat, in the dunes near the Presidio military base in San Francisco, in 1941. We sequenced the genomes of four 80- to 100-year-old Xerces Blue, and seven historical and one modern specimens of its closest relative, the Silvery Blue (Glaucopsyche lygdamus). We compared these to a novel annotated genome of the Green-Underside Blue (Glaucopsyche alexis). Phylogenetic relationships inferred from complete mitochondrial genomes indicate that Xerces Blue was a distinct species that diverged from the Silvery Blue lineage at least 850,000 years ago. Using nuclear genomes, both species experienced population growth during the Eemian interglacial period, but the Xerces Blue decreased to a very low effective population size subsequently, a trend opposite to that observed in the Silvery Blue. Runs of homozygosity and deleterious load in the former were significantly greater than in the later, suggesting a higher incidence of inbreeding. These signals of population decline observed in Xerces Blue could be used to identify and monitor other insects threatened by human activities, whose extinction patterns are still not well known.
Introduction
The Xerces Blue butterfly (Glaucopsyche xerces) (Boisduval, 1852) was native to the coastal sand dunes of San Francisco in association with the common Deerwood (Acmispon glaber), which was the preferred food source for larval stage (Tilden, 1956). It was notable for its iridescent blue colouration on the dorsal (upper) wing surface, and conspicuous, variable white spots on the ventral surface (Downey, 1956). With the growth of San Francisco and the destruction of sand dune habitats, the Xerces Blue became restricted to a few sites in what is now Golden Gate National Recreation Area. The last specimens were reportedly collected by entomologist W. Harry Lange on 23 March 1941 (Downey, 1956). It is considered the first butterfly to have been driven to global extinction by human activities (Downey, 1956).
The Xerces Blue and the closely related Silvery Blue (Glaucopsyche lygdamus) were recently proposed to be distinct species based on mtDNA data from a single Xerces Blue specimen (Grewe et al., 2021). However, two nuclear genes analysed (ribosomal 28S and histone H3) were invariable and genome-wide data were unavailable for the Xerces Blue, hampered by the inherent difficulties of retrieving genome-wide data from historical insect specimens (Thomsen et al., 2009; Staats et al., 2013) and the absence of a suitable reference genome. The genus Glaucopsyche consists of 18 extant species distributed across the temperate regions of the northern hemisphere. To provide a relevant reference, we generated an annotated genome from the Palearctic Green-Underside Blue butterfly Glaucopsyche alexis (Hinojosa Galisteo et al., 2021). Using DNA extracted from five Xerces Blue and seven Silvery Blue (G. lygdamus) historical specimens from the vicinity of San Francisco, and also from a modern Silvery Blue male from Canada, we generated whole-genome resequencing data for both species and investigated their relationships and historical population genetics.
Results
Historic and modern butterfly genomes
We extracted DNA from 12 historical specimens (5 G. xerces, 7 G. lygdamus) (Table 1). One Xerces Blue sample did not yield detectable DNA in two independent extractions. For each of the successful extracts, we prepared a single library which was shotgun sequenced on the HiseqX Illumina platform. We mapped 124,101,622 and 184,084,237 unique DNA reads of Xerces Blue and Silvery Blue, respectively, against the G. alexis reference genome (Table 2). The DNA reads exhibited typical ancient DNA features, such as short mean read length (ranging from 47.55 to 67.41 bases on average), depending on the specimen and post-mortem deamination patterns at the 5′ and 3′ ends (Supplementary file 1A). As listed in the original museum records, we found one Silvery Blue and two Xerces Blue females. Inter-individual comparisons suggested no close kinship link among the studied individuals.
The historical genomes covered 49.3% (Xerces Blue) and 55.2% (Silvery Blue) of the G. alexis reference genome, largely because repetitive chromosomal regions cannot be confidently assessed with short, ancient DNA sequence reads (Supplementary file 1B). To estimate the mappable fraction of the reference G. alexis genome, we randomly fragmented it to 50–70 nucleotides and mapped the generated fragments back to the complete genome. An average of 57.8% of the G. alexis genome was covered with these read lengths. We suggest that reduced coverage from the historical specimens may be due to genomic divergence of G. xerces and G. lygdamus from the G. alexis reference. The annotation of genes located in those unrecoverable regions provided a putative list of 14 nuclear genes with diverse functions obtained from BLAST, that should be further explored to understand the uniqueness of the extinct species (Table 3).
Phylogenetic relationships
Maximum likelihood phylogenetic inference using whole mitochondrial genomes showed that the Xerces Blue specimens form a monophyletic clade, as do the Silvery Blue specimens (Figure 1a). We inferred a time-calibrated Bayesian phylogenetic tree from protein-coding genes analysis and 12 related butterflies in Polyommatinae subfamily (Supplementary file 1C), revealing high support for the sister group relationship (posterior probability = 1). We found the specie Shijimiaeoides (Sinia) divina inside the Glaucopsyche clade, in agreement with previous phylogenetic studies (Lukhtanov and Gagarina, 2022). Because there are no known fossils to calibrate the time since divergence, we first used a molecular clock that spanned the range of rates frequently used for arthropod mitochondrial genes (1.5–2.3% divergence/Ma). Our dated analysis yielded an origin of this subgroup of Polyommatinae at 12.4 Ma (8.82–16.27 Ma 95% HPD [highest posterior density] interval) and divergence of the Xerces Blue from the Silvery Blue at 900,000 years ago (0.61–1.19 Ma 95% HPD interval, Figure 1b). A second estimate based on larger-scale fossil-based calibrations (Espeland et al., 2018) fixed the origin of the subgroup to ca. 33 Ma (Chazot et al., 2019), inferred the subsequent divergence of the Xerces Blue and Silvery Blue to 2.40 Ma (1.95–2.73 Ma 95% HPD interval, Figure 1b). The recent speciation of Xerces and Silvery Blue is not obviously due to infection with the Wolbachia, as no evidence of infection of the sampled specimens with this alpha-proteobacterium is detected in the raw read data.
Principal component analysis (PCA) using PCAngsd (Meisner and Albrechtsen, 2018) and nuclear genome polymorphisms for the three Glaucopsyche species supports the relationships among them; the historical specimens are equally distant to the Green-Underside Blue (G. alexis) in the first principal components (PC), explaining 52.81% of the variance (Figure 2). The second PC separates the Xerces Blue from the Silvery Blue specimens.
Demographic history and diversity
We used the pairwise sequentially Markovian coalescent (PSMC) algorithm (Li and Durbin, 2011) to evaluate the demographic histories of both butterfly species, first exploring the two specimens with highest coverage (L05 and L13) (Figure 3). We found an increase in effective population size in both species that is roughly coincident with the interglacial Marine Isotopic Stage 7 (approximately from 240,000 to 190,000 years ago; Batchelor et al., 2019). After this timepoint the trends differ. We estimated a continuous decrease in Xerces Blue population size in parallel to the Wisconsin Glacial Episode, which started about 75,000 years ago. However, both the modern and the historical Silvery Blue do not appear to have been negatively affected by this event (Figure 3—figure supplement 1), suggesting different adaptive strategies to cope with cooling temperatures and/or food plant availability.
Second, we generated PSMC curves from the remaining lower-coverage individuals and down-sampled data from specimen L05 to 50% and 75% of the total coverage to explore the effects of coverage on estimation of heterozygous sites. Although there was a reduction in the effective population size estimates, as expected, the temporal trajectories in lower-coverage individuals were similar to their respective, higher-coverage Xerces Blue and Silvery Blue references (Figure 3—figure supplement 1).
We subsequently explored the heterozygosity of each individual and found that Xerces Blue had 22% less heterozygosity on average than the Silvery Blue historical samples, a difference that is statistically significant (T-test; p = 0.0072) (Figure 4, Supplementary file 1D). We searched for runs of homozygosity (RoH) that can indicate the existence of inbreeding in a dwindling population. The total fraction of the genome presenting RoH, although limited, is much higher in Xerces Blue (up to 6% of the genome) than in Silvery Blue, especially in short RoH of size between 100 and 500 kb (Figure 4—figure supplement 1), consistent with background inbreeding. The limited presence of long RoH discards consanguinity as a common scenario in Xerces Blue.
We identified amino acid-changing alleles that may be suggestive of a deleterious genetic load associated with long-term low population numbers in the Xerces Blue. The average Ka/Ks ratio is higher in Xerces Blue than in Silvery Blue; the former also carries a higher fraction of nonsense and functionally high-to-moderate effect variants in homozygosity and RoH with an increased concentration of high-to-moderate effect variants (Figure 5), as predicted with a functional prediction toolbox, SnpEff (Cingolani et al., 2012).
Discussion
We have used a modern reference genome and ancient DNA genome sequence data from museum specimens to explore the relationships and historical population genetic history of an extinct butterfly, the Xerces Blue; to our knowledge, this is the first ancient genome ever generated from an extinct insect. Based upon a near-complete mtDNA genome from a Xerces Blue specimen (Grewe et al., 2021) proposed that the Xerces Blue and the Silvery Blue were distinct species. We confirm this finding using full mitochondrial genomes and extensive nuclear genomic data from multiple specimens. Given the lack of evidence for Wolbachia infection, the recent speciation of Xerces Blue and Silvery Blue seems unrelated to cytoplasmic incompatibility cause by this endosymbiont (Telschow et al., 2005; Sucháčková Bartoňová et al., 2021); a detailed analysis of genomic architectures could help identify barriers to introgression between these species.
Our analyses indicate that the Xerces Blue had experienced a severe demographic decline for tens of thousands of years, likely associated with changing climatic factors. Thus, the destruction of the Xerces Blue habitat by humans was likely the final blow in the extinction process. We provide evidence for low population size in Xerces Blue, correlated with low genetic variation, a higher proportion of RoH and increased frequency of deleterious, amino acid-changing alleles (Szpiech et al., 2013; Spielman et al., 2004; Palkopoulou et al., 2015). However, there was no genetic evidence of recent inbreeding.
Inbreeding genetic signals in the form of long chromosomal sections with no variation sometimes occur in critically endangered species (van der Valk et al., 2019; Díez-Del-Molino et al., 2018) and in extinct species such as the last Mammoths from Wrangel Island (Rogers and Slatkin, 2017) or the Altai Neanderthal (Prüfer et al., 2014). The PSMC shows a continuous low effective population size for Xerces Blue; demographic declines are also seen in some extinct species, including Wrangel Mammoths (Palkopoulou et al., 2015) but not in others such as the Woolly Rhino that showed a pre-extinction demographic stability and relatively low inbreeding signals (Lord et al., 2020). In many endangered species there is little concordance between genome diversity, population sizes, and conservation status (Díez-Del-Molino et al., 2018); this decoupling was also observed in the genomes of the extinct passenger pigeon that despite being one of the world’s most numerous vertebrates, showed a surprisingly low genetic diversity (Murray et al., 2017). Despite being notoriously abundant, insects, and in particular butterflies, are very sensitive to climate fluctuations; therefore, we suggest that insects with observations of demographic traits indicative of long-term low effective population size such as those found in Xerces Blue should be considered to be especially vulnerable to extinction events.
Our study further demonstrates the value of museum insect specimens for estimating temporal changes in genetic diversity at a population scale (Lalueza-Fox, 2022). Despite being notoriously abundant, insects, and in particular butterflies, are very sensitive to climate fluctuations. We suggest that insects with genetic observations of long-term low effective population size such as those found in Xerces Blue should be considered to be especially vulnerable to extinction events. However, being the insect numbers usually very high, it is likely that their genomic signals of extinction could be different to those described in vertebrates in many cases. Therefore, this is a subject that should be further explored with genomic data from other declining insects.
Methods
Historical butterfly specimens
The Xerces Blue specimens analysed belong to the Barnes collection deposited at the Smithsonian National Museum of Natural History. Two of them were collected on 26 April 1923. The Silvery Blue specimens were mostly collected between 1927 and 1948, in Haywood City, Santa Cruz, Oakland, San José, Fairfax, and Marin County (these locations surround San Francisco Bay) (Table 1).
DNA extraction and sequencing of Xerces Blue and Silvery Blue specimens
All DNA extraction and initial library preparation steps (prior to amplification) were performed in a dedicated clean lab, physically isolated from the laboratory used for post-polymerase chain reaction (PCR) analyses. Strict protocols were followed to minimise the amount of human DNA in the ancient DNA laboratory, including the wearing a full body suit, sleeves, shoe covers, clean shoes, facemask, hair net, and double gloving, as well as frequent bleach cleaning of benches and instruments. DNA extraction was performed from 12 abdominal samples of historical Xerces Blue and Silvery Blue, as well as a modern Silvery Blue specimen from Canada.
For the extraction procedure, 1 ml of digestion buffer (final concentrations: 3 mM CaCl2, % SDS (sodium dodecyl sulfate), 40 mM DTT (dithiothreitol), 0.25 mg/ml proteinase K, 100 mM Tris buffer pH 8.0, and 100 mM NaCl) was added to each crushed butterfly residue, including an extraction blank, and incubated at 37°C overnight (24 hr) on rotation (750–900 rpm). Next, DNA extraction was continued following the method proposed by Dabney et al., 2013. Remaining butterfly sample was then pelleted by centrifugation in a bench-top centrifuge for 2 min at maximum speed (16,100 × g). The supernatant was added to 10 ml of binding buffer (final concentrations: 5 M guanidine hydrochloride, 40% (vol/vol) isopropanol, 0.05% Tween-20, and 90 mM sodium acetate (pH 5.2)) and purified on a High Pure Extender column (Roche). DNA extracts were eluted with 45 μl of low EDTA (Tris-ethylene-diamine-tetraacetic acid) TE buffer (pH 8.0) and quantified using a Qubit instrument.
Following extraction, the DNA extract was converted into Illumina sequencing libraries following the BEST protocol (Carøe, 2018). Each library was amplified by PCR using two uniquely barcoded primers, prior to being purified with a 1.5x AMPure clean (Beckman Coulter) and eluted in 25 μl of low EDTA TE buffer (pH 8.0). One Xerces Blue sample did not yield detectable DNA in two independent extractions. For each of the successful extracts, we prepared a single library which was shotgun sequenced on the HiseqX Illumina platform.
G. alexis genome sequencing and annotation
G. alexis was chosen as a congeneric reference to compare the demographic histories of both the Xerces Blue and the Silvery Blue. We generated a G. alexis reference genome from a male specimen collected in Alcalá de la Selva in Teruel (Spain). Its genome has a sequence length of 619,543,730 bp on 24 chromosomes – including the Z sex chromosome – and the mitochondrial genome. The genome sequence is biologically complete (BUSCO Lepidoptera completeness 97.1%) (Manni et al., 2021). The G. alexis genome was sequenced at the Sanger Institute as part of the Darwin Tree of Life Project following the extraction, sequencing, and assembly protocols developed for Lepidoptera (Hinojosa Galisteo et al., 2021).
Xerces Blue and Silvery Blue mapping and variant calling
The ancient DNA reads were clipped using AdapterRemoval2 (Schubert et al., 2016), and only reads longer than 25 bp were kept. Filtered reads were mapped against the G. alexis assembly with Burrows-Wheeler Aligner (BWA) (Li and Durbin, 2009) backtrack algorithm, with parameters optimised for the analysis of aDNA (-l 2, -n 0.01, -o 2). After mapping, duplicated reads were removed using picard MarkDuplicates. Mapped reads with mapping quality below 30 were removed using samtools. Finally, to avoid problems in the next steps derived from spurious callings due to aDNA at reads’ ends, we trimmed 2 nt from each read end using BamUtil trimbam. Basic mapping statistics were generated using Qualimap2 (Okonechnikov et al., 2016; Supplementary file 1). We used bedtools (Quinlan and Hall, 2010) to assess genome coverage across the reference using windows of 1 mbp for the nuclear fraction of the genome, as well as depth of coverage, read length, and edit distance distribution. Authenticity of the sequences was assessed by characterising aDNA damage patterns with pmdtools (Skoglund et al., 2014) and MapDamage2 (Jónsson et al., 2013).
We used snpAD (Prüfer, 2018), a program for genotype calling in ancient specimens. The mapped sequences were transformed from bam-format into snpAD-format files, priors for base composition estimated, and genotypes were called using standard settings. The variant call formats (VCFs) were combined and concatenated with CombineVariants and GatherVcfs from GATK (McKenna et al., 2010) and filtered with vcftools (Danecek et al., 2011) to keep only sites within the mappable fraction of the genome previously obtained with minimum read depth of 2, max read depth of 30, genotype quality >30, maximum missingness of 0.6, minor allele frequency of 5%, and excluding indels and multiallelic sites. Since RVcoll10-B005 is a modern individual, we proceed to map it against the G. alexis reference genome with slightly altered parameters. As with the historical samples, pair was collapsed using AdapterRemoval2. BWA mem with default parameters was used as the mapping algorithm. As with the other samples, reads were filtered with Samtools (min. quality of 30) and duplicates were removed by coordinate with picard.
Genotype likelihoods were obtained with ANGSD (Korneliussen et al., 2014) using the GATK model with the following parameters for all the samples: -uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 -trim 10 -C 50 -baq 1 -minInd 5 -skipTriallelic 1 -GL 2 -minMapQ 30.
Sex determination
The sex of the specimens was determined by differential coverage of the Lepidopteran Z chromosome (females are the heterogametic sex in the Lepidoptera and show reduced coverage on the Z chromosome) (Supplementary file 1).
Mitochondrial phylogenetic tree and divergence dating
Haploid variants were called using bcftools (Danecek et al., 2021) with a ploidy of 1, filtering low-quality indels and variants, after which a consensus sequence was exported. We downloaded 14 complete mitochondrial genomes for Polyommatinae from NCBI (Supplementary file 1).
All mitochondrial genomes were annotated with MitoFinder (Allio et al., 2020) using Shijimiaeoides divina as the reference. The 11 protein-coding genes were aligned with the codon-aware aligner MACSE (Ranwez et al., 2018) and the ribosomal rRNAs were aligned with MAFFT l-ins-i (Katoh and Standley, 2013). We first investigated phylogenetic relationships among five G. xerces and eight G. lygdamus individuals, with G. alexis as the outgroup. We used IQ-TREE2 (Minh et al., 2020) to select the best fitting nucleotide substitution model for each partition and merge similar partitions (Kalyaanamoorthy et al., 2017), built a maximum likelihood tree and assessed support with 1000 ultrafast bootstrap replicates (Hoang et al., 2018).
To infer a time-calibrated phylogenetic hypothesis, we selected one individual of Xerces Blue (L003) and Silvery Blue (RVcoll10-B005) and analysed with 13 other Polyommatinae species. We used BEAST2 (Ranwez et al., 2018) with the bModelTest (Bouckaert and Drummond, 2017) package to perform phylogenetic site model averaging for each of the merged partitions. Because there is no accepted molecular clock rate for butterflies and no fossils to apply in this part of the phylogeny, we used two strategies to apply time constraints to the analysis. First, we used two published molecular clock rates for the mitochondrial COX1 gene (1.5% divergence/Ma) estimated for various invertebrates (Quek et al., 2004), and the ‘standard’ insect mitochondrial clock (2.3% divergence/Ma) (Van Zandt Brower, 1994). We applied a strict clock with a normal prior set up to span 1.5–2.3% with the 95% HPD interval (mean = 1.9%, sigma = 0.00119). Second, we borrowed the age of the most recent common ancestor of our sampled taxa from fossil-calibrated analyses across butterflies (Chazot et al., 2019; Wiemers et al., 2020). We fixed the root age to 33 Ma and allowed the remaining node ages to be estimated using a strict clock. Analyses were run twice from different starting seeds for 10 million MCMC generations and trees were sampled every 1000 generations. Runs were checked for convergence with Tracer and all effective sample size values were >200. Runs were combined with the BEAST2 package LogCombiner (Drummond and Rambaut, 2007), after removing the first 10% of topologies as burn-in, and a maximum credibility tree was generated with TreeAnnotator (Drummond and Rambaut, 2007). Phylogenetic analyses were performed on the National Life Science Supercomputing Center – Computerome 2.0 (https://www.computerome.dk/).
Xerces Blue and Silvery Blue population histories
We used the PSMC model (Li and Durbin, 2011) to explore the demographic history of both butterfly species. We obtained a consensus fastq sequence of the mappable fraction of the genome for each autosomal chromosome (total of 22 chromosomes of G. alexis assembly). Only positions with a depth of coverage above 4× and below 15× were kept. Posteriorly, a PSMC was built using the following parameters: -N25 -t15 -r5 -p ‘28*2+3+5". We used 1 year for the generation time and a mutation rate of 1.9 × 10−9, estimated in Heliconius melpomene (Martin et al., 2016). Considering that calling consensus sequences from low-coverage samples (<10×) can underestimate heterozygous sites (Keightley et al., 2015), and given the different coverage between samples, we corrected by false negative rate the samples with coverage lower than the coverage of L005 (for Xerces Blue) and L013 (for Silvery Blue), as recommended by the developers of the software, so that all samples are comparable with each other. However, since in our dataset we do not reach a coverage >20×, we acknowledge that we are not capturing the whole diversity and thus our PSMC might infer lower historical effective population sizes.
Population stratification and average genome heterozygosity
PCA was performed using PCAngsd (Meisner and Albrechtsen, 2018) after obtaining genotype likelihoods with ANGSD including all individuals. To assess global levels of heterozygosity, the unfolded site frequency spectrum (SFS) was calculated for each sample separately using ANGSD (Korneliussen et al., 2014) and realSFS with the following quality filter parameters: -uniqueOnly 1 - remove_bads 1 -only_proper_pairs 1 -trim 10 -C 50 -baq 1 -minMapQ 30 -minQ 30 -setMaxDepth 200 - doCounts 1 -GL 2 -doSaf 1.
Runs of homozygosity
RoH were called based on the density of heterozygous sites in the genome using the implemented hidden Markov model in bcftools (Danecek et al., 2021) roh with the following parameters: -G30 --skip-indels --AF-dflt 0.4 --rec-rate 1e−9 from the mappable fraction of the genome with the filtered VCF file. We kept the RoH with a phred score >85. We divided the RoH into different size bins: very short RoH (<100 kb), short RoH (100–500 kb), intermediate RoH (500 kb to 1 Mb), and long (1–5 or >5 Mb). Short RoH reflect LD patterns, intermediate size RoH describe background inbreeding due to genetic drift, and long RoH appear in the case of recent inbreeding (Ceballos et al., 2018).
Deleterious load
We used the G. alexis annotations to create a SNPeff database that we used to annotate our callings. Using SNPeff (Cingolani et al., 2012) again and the set of variants discovered by angsd, we predicted the putative effect of those variant in the analysed individuals (Supplementary file 1). In addition to wide genome mutations, we specifically focussed on mutations present in homozygosis, heterozygosis, and the previously annotated RoH.
Unrecoverable regions
To further explore how the genomic divergence can influence our genome reconstruction success, we undertook a similar approach as the genome of the Christmas Island rat (Lin et al., 2022), and explored the chromosomal regions in the G. alexis reference that were significantly depleted of Xerces DNA reads. We used bedtools (Quinlan and Hall, 2010) and some in-home bash scripting to calculate the mean coverage per gene of the G. alexis genome for Xerces Blue sequencing DNA reads. We first used bedtools’ algorithms bamtobed and genomecov to estimate the genome-wide per-site coverage of the reference genome in these two species. Then, we extracted the coordinates of all protein-coding genes from the annotation file (gff file) and used the intersect to estimate the average coverage of each protein-coding gene. We performed a functional analysis of all genes uncovered in G. alexis, excluding those that are present in G. lygdamus with more than 5× coverage (as we were looking for evolutionary novelties in the Xerces Blue lineage alone) using profile- IntersProScan (Jones et al., 2014) and sequence similarity-based (blasp) searches (Gish and States, 1993; Supplementary file 1).
Colouration genes variability
To find possible amino acid-changing variants that could explain phenotypical differences between G. lygdamus and G. xerces, we have identified and explored three well-known genes associated to colour patterns in butterflies: optix, cortex, and Wnt genes (Zhang and Reed, 2016; Zhang et al., 2017; Mazo-Vargas et al., 2017; Fenner et al., 2020; Banerjee et al., 2021). First, we located those genes in our annotation with BLAST and their homologs in other butterfly species, setting an E-value lower than 0.001 and an Identity value above 60% (Table 3). Then, the coordinates were called using GATK UnifiedGenotyper. Variants were filtered for indels and minimum Genotype Quality of 30 using. Variants were kept regardless of their coverage. A variant is considered as fixed in one species if it is covered in at least two individuals of each species, it is in homozygous state, and when one of the species present all their genotypes calls as homozygous for the alternative allele while in the other are homozygous of the reference allele. No fixated mutations were identified in the regions covered at the same time by G. lygdamus and G. xerces sequences.
Wolbachia screening
Wolbachia are endosymbiotic alpha-proteobacteria that are present in about 70% of butterfly species and induce diverse reproductive alterations, including genetic barriers when two different strains infect the same population or when two populations – one infected and one uninfected – meet (Telschow et al., 2005). As potential evidence for a reproductive barrier promoting the separation of Xerces Blue and Silvery Blue, we searched for Wolbachia DNA reads in our specimens, taking advantage of the high coverage and the shotgun approach. First, we collapsed unique reads from the butterfly-free sequences with BBmap (Bushnell, 2014) and removed from the dataset low complexity sequences using Prinseq (Schmieder and Edwards, 2011). Afterwards, we used Kraken2 (Wood et al., 2019) to assign reads against the standard plus human Kraken2 database (bacteria, archaea, fungi, protozoa, and viral). The historical specimens did not display enough reads assigned to Wolbachia for us to suspect of the presence of the bacteria in those samples (Table 4).
Data availability
The accession numbers for the Xerces Blue and Silvery Blue genomes reported in this study are in the European Nucleotide Archive (ENA): PRJEB47122. Data on G. alexis are available in INSDC under BioProject PRJEB43798 and genome assembly accessions GCA_905404095.1 (primary haplotype) and GCA_905404225.1 (secondary, alternate haplotype).
-
European Nucleotide ArchiveID PRJEB47122. Whole-genomes from the extinct Xerces Blue butterfly reveal low diversity and long-term population decline.
-
NCBI BioProjectID PRJEB43798. ilGlaAlex (green-underside blue).
-
NCBI GenBankID GCA_905404095.1. Genome assembly ilGlaAlex1.1.
-
NCBI GenBankID GCA_905404225.1. Genome assembly ilGlaAlex1.1 alternate haplotype.
References
-
MitoFinder: efficient automated large-scale extraction of mitogenomic data in target enrichment phylogenomicsMolecular Ecology Resources 20:892–905.https://doi.org/10.1111/1755-0998.13160
-
Basic local alignment search toolJournal of Molecular Biology 215:403–410.https://doi.org/10.1016/S0022-2836(05)80360-2
-
The configuration of Northern Hemisphere ice sheets through the QuaternaryNature Communications 10:3713.https://doi.org/10.1038/s41467-019-11601-2
-
bModelTest: Bayesian phylogenetic site model averaging and model comparisonBMC Evolutionary Biology 17:42.https://doi.org/10.1186/s12862-017-0890-6
-
BEAST 2.5: An advanced software platform for Bayesian evolutionary analysisPLOS Computational Biology 15:e1006650.https://doi.org/10.1371/journal.pcbi.1006650
-
ConferenceBBMap: A fast, accurate, splice-aware aligner9th Annual Genomics of Energy & Environment Meeting.
-
Single-tube library preparation for degraded DNAMethods in Ecology and Evolution 9:410–419.https://doi.org/10.1111/2041-210X.12871
-
Runs of homozygosity: windows into population history and trait architectureNature Reviews. Genetics 19:220–234.https://doi.org/10.1038/nrg.2017.109
-
The variant call format and VCFtoolsBioinformatics 27:2156–2158.https://doi.org/10.1093/bioinformatics/btr330
-
Twelve years of SAMtools and BCFtoolsGigaScience 10:giab008.https://doi.org/10.1093/gigascience/giab008
-
Quantifying temporal genomic erosion in endangered speciesTrends in Ecology & Evolution 33:176–185.https://doi.org/10.1016/j.tree.2017.12.002
-
Analysis of Variation in a recently extinct polymorphic lycaenid butterfly, Glaucopsyche XercesBull. So. Calif. Acad. Sci 55:153–170.
-
BEAST: Bayesian evolutionary analysis by sampling treesBMC Evolutionary Biology 7:214.https://doi.org/10.1186/1471-2148-7-214
-
De novo genome assemblies of butterfliesGigaScience 10:giab041.https://doi.org/10.1093/gigascience/giab041
-
A comprehensive and dated phylogenomic analysis of butterfliesCurrent Biology 28:770–778.https://doi.org/10.1016/j.cub.2018.01.061
-
Wnt genes in wing pattern development of coliadinae butterfliesFrontiers in Ecology and Evolution 8:197.https://doi.org/10.3389/fevo.2020.00197
-
UFBoot2: improving the ultrafast bootstrap approximationMolecular Biology and Evolution 35:518–522.https://doi.org/10.1093/molbev/msx281
-
InterProScan 5: genome-scale protein function classificationBioinformatics 30:1236–1240.https://doi.org/10.1093/bioinformatics/btu031
-
MAFFT multiple sequence alignment software version 7: improvements in performance and usabilityMolecular Biology and Evolution 30:772–780.https://doi.org/10.1093/molbev/mst010
-
Estimation of the spontaneous mutation rate in Heliconius melpomeneMolecular Biology and Evolution 32:239–243.https://doi.org/10.1093/molbev/msu302
-
ANGSD: Analysis of Next Generation Sequencing DataBMC Bioinformatics 15:356.https://doi.org/10.1186/s12859-014-0356-4
-
The Sequence Alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
Probing the genomic limits of de-extinction in the Christmas Island ratCurrent Biology 32:1650–1656.https://doi.org/10.1016/j.cub.2022.02.027
-
BUSCO: assessing genomic data quality and beyondCurrent Protocols 1:e323.https://doi.org/10.1002/cpz1.323
-
IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic EraMolecular Biology and Evolution 37:1530–1534.https://doi.org/10.1093/molbev/msaa015
-
SNPAD: an ancient DNA genotype callerBioinformatics 34:4165–4171.https://doi.org/10.1093/bioinformatics/bty507
-
Codiversification in an ant-plant mutualism: stem texture and the evolution of host use in Crematogaster (Formicidae: Myrmicinae) inhabitants of Macaranga (Euphorbiaceae)Evolution; International Journal of Organic Evolution 58:554–570.https://doi.org/10.1111/j.0014-3820.2004.tb01678.x
-
Posterior summarization in bayesian phylogenetics using tracer 1.7Systematic Biology 67:901–904.https://doi.org/10.1093/sysbio/syy032
-
MACSE v2: Toolkit for the alignment of coding sequences accounting for frameshifts and stop codonsMolecular Biology and Evolution 35:2582–2584.https://doi.org/10.1093/molbev/msy159
-
SoftwareA language and environment for statistical computingR Foundation for Statistical Computing, Vienna, Austria.
-
Long runs of homozygosity are enriched for deleterious variationAmerican Journal of Human Genetics 93:90–102.https://doi.org/10.1016/j.ajhg.2013.05.003
-
The effect of Wolbachia versus genetic incompatibilities on reinforcement and speciationEvolution; International Journal of Organic Evolution 59:1607–1619.
-
Phylogeny of Heliconius butterflies inferred from mitochondrial DNA sequences (Lepidoptera: Nymphalidae)Molecular Phylogenetics and Evolution 3:159–174.https://doi.org/10.1006/mpev.1994.1018
-
Improved metagenomic analysis with Kraken 2Genome Biology 20:257.https://doi.org/10.1186/s13059-019-1891-0
Article and author information
Author details
Funding
Ministerio de Ciencia e Innovación (PID2021-124590NB-100)
- Carles Lalueza-Fox
European Research Council (864203)
- Tomàs Marquès
Ministerio de Ciencia e Innovación (BFU2017-86471-P)
- Tomàs Marquès
Ministerio de Ciencia e Innovación (CEX2018-000792-M)
- Tomàs Marquès
Generalitat de Catalunya (GRC 2017-SGR-880)
- Tomàs Marquès
Wellcome Trust
https://doi.org/10.35802/206194- Marcela Uliano-Silva
Wellcome Trust
https://doi.org/10.35802/218328- Mark Blaxter
Ministerio de Ciencia e Innovación (PID2019-107078GB-I00)
- Roger Vila
Generalitat de Catalunya (GRC 2017-SGR-991)
- Roger Vila
Howard Hughes Medical Institute (Howard Hughes International Early Career)
- Tomàs Marquès
The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.
Acknowledgements
CL-F is supported by a PID2021-124590NB-100 grant (MCIU/AEI/FEDER, UE) of Spain; TM-B is supported by funding from the European Research Council (ERC) (grant agreement No. 864203), BFU2017-86471-P (MINECO/FEDER, UE), 'Unidad de Excelencia María de Maeztu', funded by the AEI (CEX2018-000792-M), Howard Hughes International Early Career, and Generalitat de Catalunya, GRC 2017-SGR-880; RV is supported by grant PID2019-107078GB-I00, funded by MCIN/AEI/10.13039/501100011033, and by GRC 2017-SGR-991 (Generalitat de Catalunya). We are grateful to the SCIENCE Faculty at University of Copenhagen for free access to Computerome 2.0. This research was funded in whole, or in part, by the Wellcome Trust Grants 206194 and 218328 (MU, MB). For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
- Version of Record published:
- Version of Record updated:
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.87928. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, de-Dios, Fontsere, Renom et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,247
- views
-
- 159
- downloads
-
- 8
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Evolutionary Biology
Spatial patterns in genetic diversity are shaped by individuals dispersing from their parents and larger-scale population movements. It has long been appreciated that these patterns of movement shape the underlying genealogies along the genome leading to geographic patterns of isolation-by-distance in contemporary population genetic data. However, extracting the enormous amount of information contained in genealogies along recombining sequences has, until recently, not been computationally feasible. Here, we capitalize on important recent advances in genome-wide gene-genealogy reconstruction and develop methods to use thousands of trees to estimate per-generation dispersal rates and to locate the genetic ancestors of a sample back through time. We take a likelihood approach in continuous space using a simple approximate model (branching Brownian motion) as our prior distribution of spatial genealogies. After testing our method with simulations we apply it to Arabidopsis thaliana. We estimate a dispersal rate of roughly 60 km2/generation, slightly higher across latitude than across longitude, potentially reflecting a northward post-glacial expansion. Locating ancestors allows us to visualize major geographic movements, alternative geographic histories, and admixture. Our method highlights the huge amount of information about past dispersal events and population movements contained in genome-wide genealogies.
-
- Evolutionary Biology
Understanding the genomic basis of natural variation in plant pest resistance is an important goal in plant science, but it usually requires large and labor-intensive phenotyping experiments. Here, we explored the possibility that non-target reads from plant DNA sequencing can serve as phenotyping proxies for addressing such questions. We used data from a whole-genome and -epigenome sequencing study of 207 natural lines of field pennycress (Thlaspi arvense) that were grown in a common environment and spontaneously colonized by aphids, mildew, and other microbes. We found that the numbers of non-target reads assigned to the pest species differed between populations, had significant SNP-based heritability, and were associated with climate of origin and baseline glucosinolate contents. Specifically, pennycress lines from cold and thermally fluctuating habitats, presumably less favorable to aphids, showed higher aphid DNA load, i.e., decreased aphid resistance. Genome-wide association analyses identified genetic variants at known defense genes but also novel genomic regions associated with variation in aphid and mildew DNA load. Moreover, we found several differentially methylated regions associated with pathogen loads, in particular differential methylation at transposons and hypomethylation in the promoter of a gene involved in stomatal closure, likely induced by pathogens. Our study provides first insights into the defense mechanisms of Thlaspi arvense, a rising crop and model species, and demonstrates that non-target whole-genome sequencing reads, usually discarded, can be leveraged to estimate intensities of plant biotic interactions. With rapidly increasing numbers of large sequencing datasets worldwide, this approach should have broad application in fundamental and applied research.