Whole genomes from the extinct Xerces Blue butterfly can help identify declining insect species

  1. Toni de Dios Martinez
  2. Claudia Fontsere
  3. Pere Renom
  4. Josefin Stiller
  5. Laia Llovera
  6. Marcela Uliano-Silva
  7. Alejandro Sánchez-Gracia
  8. Charlotte Wright
  9. Esther Lizano
  10. Berta Caballero
  11. Arcadi Navarro
  12. Sergi Civit
  13. Robert K Robbins
  14. Mark Blaxter
  15. Tomàs Marquès  Is a corresponding author
  16. Roger Vila  Is a corresponding author
  17. Carles Lalueza-Fox  Is a corresponding author
  1. Institute of Evolutionary Biology, Spain
  2. Institute of Genomics, University of Tartu, Estonia
  3. Section for Evolutionary Genomics, The Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Denmark
  4. Centre for Biodiversity Genomics, University of Copenhagen, Denmark
  5. Wellcome Sanger Institute, United Kingdom
  6. Departament of Genetics, Microbiology and Statistics-Institut de Recerca de la Biodiversitat (IRBio), Universitat de Barcelona, Spain
  7. Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Spain
  8. Museu de Ciències Naturals de Barcelona, Spain
  9. Catalan Institution of Research and Advanced Studies (ICREA), Spain
  10. Department of Entomology, National Museum of Natural History, Smithsonian Institution, United States
  11. CNAG-CRG, Centre for Genomic Regulation, Barcelona Institute of Science and Technology (BIST), Spain

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.sa0

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.

Table 1
List of historical specimens analysed in this study.
Genome #SpeciesSubspp.StateLocalityDateCollection
USNMENT101413G. xercesCaliforniaSan FranciscoNABarnes
USNMENT101402G. xercesCaliforniaSan Francisco16/4/1923Barnes
USNMENT101441G. xercesCaliforniaSan FranciscoNABarnes
USNMENT101406G. xercesCaliforniaSan FranciscoNABarnes
USNMENT101434G. xercesCaliforniaSan Francisco16/4/1923Barnes
USNMENT00181297G. lygdamusincognitusCaliforniaMarin CountryNABarnes
USNMENT00181298G. lygdamusincognitusCaliforniaFairfax27/5/1932WMD Field
USNMENT00181299G. lygdamusincognitusCaliforniaOakland14/4/1948Graham Heid
USNMENT00181300G. lygdamusincognitusCaliforniaSan Jose27/3/1964Opler
USNMENT00181301G. lygdamusincognitusCaliforniaHaywood City1/5/1931WMD Field
USNMENT00181302G. lygdamusincognitusCaliforniaSanta Cruz1/4/1932JW Tilden/Field
USNMENT00181303G. lygdamusincognitusCaliforniaSanta Cruz8/4/1927GW Rawson
Table 2
Mapping statistics of the analysed historical specimens.

Mapping statistics of the four historical G. xerces (L003, L005, L007, and L009) and the seven historical G. lydagmus (L002, L004, L006, L008, L011, L012, and L013) specimens mapped against the G. alexis reference genome. Average depth is displayed for the covered regions of each individual.

Sample identifierGenerated readsQ25 unique mapped readsBreadth of coverage (%)Average depth covered regions
L002300,294,24823,337,75137.275.105
L003405,198,06032,547,82036.866.78
L004357,165,43828,722,18538.776.55
L005776,312,37856,459,03745.712.42
L006359,520,16828,498,72040.076.18
L007348,916,87026,758,35634.796.21
L008508,120,15632,107,19242.087.422
L009322,955,38439,312,61740.68.02
L011236,886,53424,165,28238.65.40
L012328,359,66918,683,73833.374.29
L013385,635,64452,612,93747.212.3

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).

Table 3
Coordinates of the analysed colouration genes.

Genomic coordinates in G. alexis reference genomes of different wing colouration genes described in other butterfly species.

ChromosomeStartEndGene
FR990043.15,387,7065,403,599Wnt1
FR990043.15,417,9025,423,677Wnt6
FR990043.15,519,3535,539,737Wnt10b
FR990043.15,553,6665,554,753Wnt10a
FR990043.126,972,85626,974,487WntA
FR990046.12,343,4672,357,667Wnt7b
FR990046.16,255,2756,271,623Wnt5b
FR990046.119,475,63619,486,554Wnt9
FR990050.116,200,97816,212,495Wnt11
FR990054.120,633,40020,655,261Cortex
FR990059.120,254,46020,255,275Optix

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.

Phylogenetic placement of the Xerces Blue.

(a) Maximum likelihood tree from whole mitochondrial genomes of Xerces Blue, Silvery Blue, and Green-Underside Blue. Node labels are bootstrap support values. (b) Time-calibrated phylogeny from Bayesian inference using mitochondrial protein-coding genes of Xerces Blue and related butterflies. Node values show median age estimates from dating analysis with a molecular clock (above nodes) or from fixing the age of the root (below nodes). Bars are 95% HPD intervals for node ages. All posterior probabilities were 1, except for one node annotated in black.

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.

Plotting of PC1 and PC2 of the principal component analysis (PCA).

The PCA was generated with nuclear DNA data (N = 6,682,591 SNPs (single nucleotide polymorphisms)) from 11 historical butterfly specimens (4 G. xerces and 7 G. lygdamus), a modern G. lygdamus from Canada (RVcoll10-B005) and a modern G. alexis reference genome. The PCA shows a clear separation of both historical species and the reference in the first PC (explaining 52.81% of the variance), and separation of G. xerces and G. lygdamus by the second PC (explaining 6.09% of the variance), supporting they are separated lineages.

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.

Figure 3 with 1 supplement see all
Pairwise sequentially Markovian coalescent (PSMC) plot of one Xerces Blue (Glaucopsyche xerces) (L05) specimen and one Silvery Blue specimen (Glaucopsyche lygdamus).

The two historical samples are those with higher average coverage. Individual PSMC plots were bootstrapped 100 times each (lighter lines). One year of generation time and a mutation rate of µ = 1.9 × 10−9 were used. The peak of the Marine Isotopic Stage 7 interglacial is marked in yellow.

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.

Figure 4 with 1 supplement see all
Runs of homozygosity (RoH) in the genomes of Xerces Blue and Silvery Blue (modern and historical).

(a) Percentage of the autosomal genome in RoH by size bins: very short RoH (<100 kb), short RoH (100–500 kb), intermediate RoH (500 kb to 1 Mb), and long (1–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 very recent inbreeding due to consanguinity. Error bars show the standard deviation. (b) Distribution of RoH in the autosomal genome of a Xerces specimen, L05. (c) Distribution of RoH in the autosomal genome of a Silvery specimen L13.

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).

Functional effect prediction on the fixed amino acid-changing alleles observed in Xerces Blue and Silvery Blue.

(a) Wide genome Ka/Ks ratio comparison. (b) High-to-moderate effect variant comparison in homozygous sites. (c) High-to-moderate effect variant comparison in heterozygous sites. (d) Presence of high-to-moderate variants in regions of the genome in runs of homozygosity (RoH). Error bars show the standard deviation.

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

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Biological sample (Glaucopsyche xerces; female)L003This paperSAMEA114094142See Materials and methods
Biological sample (G. xerces; male)L005This paperSAMEA114094143See Materials and methods
Biological sample (G. xerces; male)L007This paperSAMEA114094144See Materials and methods
Biological sample (G. xerces; female)L009This paperSAMEA114094145See Materials and methods
Biological sample (Glaucopsyche lygdamus; male)L002This paperSAMEA114094134See Materials and methods
Biological sample (G. lygdamus; male)L004This paperSAMEA114094135See Materials and methods
Biological sample (G. lygdamus; male)L006This paperSAMEA114094136See Materials and methods
Biological sample (G. lygdamus; male)L008This paperSAMEA114094137See Materials and methods
Biological sample (G. lygdamus; male)L011This paperSAMEA114094138See Materials and methods
Biological sample (G. lygdamus; female)L012This paperSAMEA114094139See Materials and methods
Biological sample (G. lygdamus; male)L013This paperSAMEA114094140See Materials and methods
Biological sample (G. lygdamus; male)RVcoll10-B005This paperSAMEA114094141See Materials and methods
Biological sample (Glaucopsyche alexis; male)G. alexisHinojosa Galisteo et al., 2021ilGlaAlex1.1; GCA_905404095.1
Biological sample (Aricia agestis)A. agestisHayward et al., 2023LR990279.1
Biological sample (Aricia artaxerxes)A. artaxerxesEbdon et al., 2022OW569311.1
Biological sample (Celastrina argiolus)C. argiolusHayward et al., 2021LR994603.1
Biological sample (Cyaniris semiargus; male)C. semiargusLohse et al., 2023LR994570.1
Biological sample (G. alexis; male)G. alexisHinojosa Galisteo et al., 2021FR990065.1
Biological sample (G. xerces)G. xercesGrewe et al., 2021MW677564.1
Biological sample (Lysandra bellargus; female)L. bellargusLohse et al., 2022HG995365.1
Biological sample (Lysandra coridon; male)L. coridonVila et al., 2023HG992145.1
Biological sample (Plebejus argus)P. argusZhou et al., 2020MN974526.1
Biological sample (Plebejus melissa)P. melissaEllis et al., 2021DWQ001000057.1
Biological sample (Plebejus anna)P. annaEllis et al., 2021DWTA01000073.1
Biological sample (Polyommatus icarus; male)P. icarushttps://www.darwintreeoflife.org/OW569343.1
Biological sample (Shijimiaeoides divina)S. divinaJeong et al., 2017NC_029763.1
Biological sample (Zizina emelina)Z. emelinaLiu et al., 2020MN013031.1
Software, algorithmBUSCOManni et al., 2021v.5.1.2
Software, algorithmAdapterRemovalSchubert et al., 2016v.2.2.2
Software, algorithmBWA – backtrackLi and Durbin, 2009v.0.7.1
Software, algorithmBWA – memLi, 2013v.0.7.1
Software, algorithmQualimap2Okonechnikov et al., 2016v.2.2.2
Software, algorithmpmdtoolsSkoglund et al., 2014v.0.50
Software, algorithmMapDamage2Jónsson et al., 2013v.2.7.12
Software, algorithmBedtoolsQuinlan and Hall, 2010v.2.27.1
Software, algorithmsnpADPrüfer, 2018v.0.3.2
Software, algorithmGATKMcKenna et al., 2010v.3.5–3.7
Software, algorithmvcftoolsDanecek et al., 2011v.0.1.12b–0.1.14b
Software, algorithmangsdKorneliussen et al., 2014v.0.916
Software, algorithmbcftoolsDanecek et al., 2021v.1.9
Software, algorithmMitofinderAllio et al., 2020v.1.4
Software, algorithmMACSERanwez et al., 2018v.2.05
Software, algorithmMAFFTKatoh and Standley, 2013v.7.490
Software, algorithmIQ-TREE2Minh et al., 2020v.2.1.3
Software, algorithmModelFinderKalyaanamoorthy et al., 2017Available in IQ-TREE2
Software, algorithmUFBoot2Hoang et al., 2018Available in IQ-TREE2
Software, algorithmBEAST2Bouckaert et al., 2019v.2.6.3
Software, algorithmbModelTestBouckaert and Drummond, 2017v.1.2.1
Software, algorithmTracerRambaut et al., 2018v.1.7.2
Software, algorithmPSMCLi and Durbin, 2011v.0.6.5
Software, algorithmPCAngsdMeisner and Albrechtsen, 2018v.20180209
Software, algorithmBcftools-rohNarasimhan et al., 2016v.1.9
Software, algorithmSNPeffCingolani et al., 2012v.4.3
Software, algorithmPicardBroad Institute, 2015v.2.0.1
Software, algorithmSamtoolsLi et al., 2009v.1.6
Software, algorithmBamUtilJun et al., 2015v.1.0.13
Software, algorithmBedtoolsQuinlan and Hall, 2010v.2.27.1
Software, algorithmBLASTAltschul et al., 1990v.2.2.2
Software, algorithmBBMapBushnell, 2014v.38.18
Software, algorithmPrinseqSchmieder and Edwards, 2011v.0.20.4
Software, algorithmKraken2Wood et al., 2019v.2.1.1
Software, algorithmRR Core Team, 2019v.3.6.3–4.1.0
Software, algorithmGgplot2Wickham, 2016v.3.0.0

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).

Table 4
Wolbachia DNA reads assigned using Kraken2.
SpecimenWolbachia genus readsWolbachia spp. reads
L0021905
L0031313
L0042135
L0053118
L0062429
L0071522
L00841421
L0092366
L0111849
L0121689
L01352324

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).

The following data sets were generated

References

    1. Boisduval JA
    (1852)
    Lépidoptères de la Californie
    Ann. Soc. Ent. Fr 21:275–324.
  1. Conference
    1. Bushnell B
    (2014)
    BBMap: A fast, accurate, splice-aware aligner
    9th Annual Genomics of Energy & Environment Meeting.
    1. Downey JC
    (1956)
    Analysis of Variation in a recently extinct polymorphic lycaenid butterfly, Glaucopsyche Xerces
    Bull. So. Calif. Acad. Sci 55:153–170.
  2. Software
    1. R Core Team
    (2019) A language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    1. Telschow A
    2. Hammerstein P
    3. Werren JH
    (2005)
    The effect of Wolbachia versus genetic incompatibilities on reinforcement and speciation
    Evolution; International Journal of Organic Evolution 59:1607–1619.
    1. Tilden JW
    (1956)
    San Francisco’s vanishing butterflies
    Lepid. News 10:133–1445.

Article and author information

Author details

  1. Toni de Dios Martinez

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Institute of Genomics, University of Tartu, Tartu, Estonia
    Contribution
    Formal analysis, Investigation
    Contributed equally with
    Claudia Fontsere and Pere Renom
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9260-8846
  2. Claudia Fontsere

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Section for Evolutionary Genomics, The Globe Institute, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Methodology
    Contributed equally with
    Toni de Dios Martinez and Pere Renom
    Competing interests
    No competing interests declared
  3. Pere Renom

    Institute of Evolutionary Biology, Barcelona, Spain
    Contribution
    Conceptualization, Investigation
    Contributed equally with
    Toni de Dios Martinez and Claudia Fontsere
    Competing interests
    No competing interests declared
  4. Josefin Stiller

    Centre for Biodiversity Genomics, University of Copenhagen, Copenhagen, Denmark
    Contribution
    Formal analysis, Investigation, Visualization, Writing – original draft
    Competing interests
    No competing interests declared
  5. Laia Llovera

    Institute of Evolutionary Biology, Barcelona, Spain
    Contribution
    Resources, Formal analysis
    Competing interests
    No competing interests declared
  6. Marcela Uliano-Silva

    Wellcome Sanger Institute, Saffron Walden, United Kingdom
    Contribution
    Conceptualization, Formal analysis
    Competing interests
    No competing interests declared
  7. Alejandro Sánchez-Gracia

    Departament of Genetics, Microbiology and Statistics-Institut de Recerca de la Biodiversitat (IRBio), Universitat de Barcelona, Barcelona, Spain
    Contribution
    Methodology
    Competing interests
    No competing interests declared
  8. Charlotte Wright

    Wellcome Sanger Institute, Saffron Walden, United Kingdom
    Contribution
    Formal analysis, Investigation
    Competing interests
    No competing interests declared
  9. Esther Lizano

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain
    Contribution
    Conceptualization, Formal analysis
    Competing interests
    No competing interests declared
  10. Berta Caballero

    Museu de Ciències Naturals de Barcelona, Barcelona, Spain
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  11. Arcadi Navarro

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Catalan Institution of Research and Advanced Studies (ICREA), Barcelona, Spain
    Contribution
    Resources, Supervision, Project administration
    Competing interests
    No competing interests declared
  12. Sergi Civit

    Departament of Genetics, Microbiology and Statistics-Institut de Recerca de la Biodiversitat (IRBio), Universitat de Barcelona, Barcelona, Spain
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  13. Robert K Robbins

    Department of Entomology, National Museum of Natural History, Smithsonian Institution, Washington, United States
    Contribution
    Conceptualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  14. Mark Blaxter

    Wellcome Sanger Institute, Saffron Walden, United Kingdom
    Contribution
    Conceptualization, Supervision
    Competing interests
    No competing interests declared
  15. Tomàs Marquès

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain
    3. Catalan Institution of Research and Advanced Studies (ICREA), Barcelona, Spain
    4. CNAG-CRG, Centre for Genomic Regulation, Barcelona Institute of Science and Technology (BIST), Barcelona, Spain
    Contribution
    Resources, Supervision, Investigation, Project administration
    For correspondence
    tomas.marques@upf.edu
    Competing interests
    No competing interests declared
  16. Roger Vila

    Institute of Evolutionary Biology, Barcelona, Spain
    Contribution
    Conceptualization, Supervision, Writing – original draft
    For correspondence
    roger.vila@csic.es
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2447-4388
  17. Carles Lalueza-Fox

    1. Institute of Evolutionary Biology, Barcelona, Spain
    2. Museu de Ciències Naturals de Barcelona, Barcelona, Spain
    Contribution
    Conceptualization, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    carles.lalueza.fox@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1730-5914

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

  1. Preprint posted:
  2. Sent for peer review:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Version of Record published:
  6. 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 Martinez, 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,246
    views
  • 159
    downloads
  • 8
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Toni de Dios Martinez
  2. Claudia Fontsere
  3. Pere Renom
  4. Josefin Stiller
  5. Laia Llovera
  6. Marcela Uliano-Silva
  7. Alejandro Sánchez-Gracia
  8. Charlotte Wright
  9. Esther Lizano
  10. Berta Caballero
  11. Arcadi Navarro
  12. Sergi Civit
  13. Robert K Robbins
  14. Mark Blaxter
  15. Tomàs Marquès
  16. Roger Vila
  17. Carles Lalueza-Fox
(2024)
Whole genomes from the extinct Xerces Blue butterfly can help identify declining insect species
eLife 12:RP87928.
https://doi.org/10.7554/eLife.87928.3

Share this article

https://doi.org/10.7554/eLife.87928

Further reading

    1. Evolutionary Biology
    Matthew Osmond, Graham Coop
    Research Article Updated

    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.

    1. Evolutionary Biology
    Dario Galanti, Jun Hee Jung ... Oliver Bossdorf
    Research Article

    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.