Summer rainfall drives adaptation with gene flow in a widespread butterfly
Data files
Jun 24, 2026 version files 1.91 GB
-
colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.subset.sigSNPs.annotate.vcf
15.98 MB
-
colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.vcf
960.08 MB
-
colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gz
932.97 MB
-
README.md
8.47 KB
Abstract
Understanding how environmental variation interacts with gene flow to shape population genomic patterns is a central goal in evolutionary biology. We investigated how geographic and environmental differences impact genomic variation in the clouded sulfur butterfly (Colias philodice eriphyle) by conducting whole-genome resequencing across replicated transects consisting of paired high- and low-elevation sites on both sides of a major mountain range. Despite sampling across steep environmental gradients, we found no evidence of discrete population structure, indicating high connectivity across the region. Nonetheless, significant isolation by distance – strongest in eastern populations – revealed that geographic distance still imposes limits on gene flow, and genetic diversity was also elevated in the east. Genotype-environment association analyses identified more than 16,000 loci associated with elevation, precipitation, and solar radiation. Our redundancy analysis identified precipitation as the strongest predictor of adaptive genomic differentiation, and candidate genes included those linked to melanization and thermoregulation (e.g., TH and yellow). These results demonstrate that even in a largely panmictic population, environmental variation can maintain regional-scale signals of local adaptation. Because insects are declining globally and remain underrepresented in genomic monitoring, conducting whole-genome analyses in a widespread species provides valuable context for assessing how insects today persist across such diverse landscapes and their potential for withstanding future environmental change.
Dataset DOI: 10.5061/dryad.gmsbcc34r
Description of the data and file structure
Overview
This repository contains genomic data for Colias philodice eriphyle collected across replicated elevational transects in Colorado and Wyoming, USA (2022). Whole-genome resequencing was conducted on 165 individuals across 15 sites east and west of the North American continental divide.
Files
1. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gz Full filtered VCF file containing all 1,454,624 SNPs across 165 C. p. eriphyle individuals and 31 scaffolds. Raw FASTQ files were trimmed using Trimmomatic, mapped to the Colias crocea reference genome (PRJEB42949) using BWA-MEM, and processed using GATK v4.3.0. Includes only biallelic sites with <20% missing data per site, a minor allele frequency filter of 5%–95%, base quality score >30, and minimum mapping quality of 20. PCR duplicates were removed with SAMtools and Bam files were downsampled to 4× coverage. Three related individuals were removed using NGSRelate.
2. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.vcf VCF file used for the redundancy analysis (RDA), containing 156 individuals (9 individuals from site Wy excluded due to missing GPS coordinates). Genotypes were imputed using BEAGLE v4.1 prior to running the RDA.
3. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.subset.sigSNPs.annotate.vcf VCF file containing the 16,492 significant outlier SNPs identified by the RDA using a standard deviation cutoff of 3 along one or more RDA axes. SNPs are annotated with the nearest gene from the Colias crocea reference genome using BEDtools, retaining loci within 500 bp of annotated genes. Functional annotation was performed using SnpEff.
Reference Genome
Colias crocea reference genome: European Nucleotide Archive, accession PRJEB42949 (Ebdon et al. 2022).
Analysis Scripts
All scripts used for data analysis are available on GitHub at https://github.com/lilyd-csu/colias-gea/. This repository contains scripts for all downstream analyses conducted using the genomic data files described above, including population structure analyses, redundancy analysis (RDA), isolation by distance, gene ontology enrichment, and allele frequency regression analyses. Scripts are organized by analysis type and include documentation for reproducing the results reported in Durkee et al. 2026.
Files and variables
1. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gz Compressed VCF (Variant Call Format) file containing all 1,454,624 filtered SNPs across 165 C. p. eriphyle individuals and 31 scaffolds. Standard VCF fields:
- CHROM: scaffold/chromosome identifier (scaffolds NC_059537.1–NC_059566.1 = autosomal chromosomes; NC_059568.1 = Z sex chromosome)
- POS: position of the SNP on the scaffold (bp)
- ID: SNP identifier
- REF: reference allele
- ALT: alternate allele
- QUAL: Phred-scaled quality score
- FILTER: PASS = variant passed all filters; other values indicate the filter failed (FS60, LowQual, MQ40, MQRankSum-12.5, QD2, QUAL30, ReadPosRankSum-8, SOR3)
- INFO fields: AC (allele count), AF (allele frequency), AN (total alleles), DP (read depth), FS (strand bias), MQ (mapping quality), QD (quality by depth), SOR (strand odds ratio), NMISS (number of missing genotypes), among others
- FORMAT fields: GT (genotype), AD (allelic depth), DP (read depth), GQ (genotype quality), PL (phred-scaled genotype likelihoods), PGT (physical phasing haplotype), PID (physical phasing ID), PS (phasing set)
- Sample columns: one column per individual, named by sample ID
Missing genotypes are indicated as "./." in the GT field.
2. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.vcf VCF file used for the redundancy analysis (RDA), containing 156 individuals (9 individuals from site Wy excluded due to missing GPS coordinates). Genotypes were imputed using BEAGLE v4.1. Generated from file 1 using PLINK v1.90. Contains only:
- CHROM, POS, ID, REF, ALT: standard VCF fields as above
- INFO: PR flag (provisional reference allele)
- FORMAT: GT (genotype) only
3. colias.4x.merged_gatk.rm.relate.SNP.filtered_gatkVQSR2.PASS.8miss.recode.vcf.gl_impute4.1.rm.WY.subset.sigSNPs.annotate.vcf VCF file containing the 16,492 significant outlier SNPs identified by the RDA. Same base format as file 2, with additional SnpEff functional annotation in the INFO field:
- ANN: functional annotations including allele, annotation type, annotation impact (LOW/MODERATE/HIGH/MODIFIER), gene name, gene ID, feature type, transcript biotype, HGVS notation, and distance to nearest feature
- LOF: predicted loss-of-function effects (gene name, gene ID, number of transcripts, percent transcripts affected)
- NMD: predicted nonsense-mediated decay effects (same format as LOF)
Abbreviations:
- SNP: single nucleotide polymorphism
- VCF: variant call format
- bp: base pairs
- RDA: redundancy analysis
- MAF: minor allele frequency
- GATK: Genome Analysis Toolkit
- HGVS: Human Genome Variation Society (nomenclature for variant annotation)
Code/software
All analysis scripts are available at https://github.com/lilyd-csu/colias-gea/. The following software and packages were used:
Bioinformatic pipeline (run on RMACC Alpine computing cluster, Colorado State University/University of Colorado):
- Trimmomatic v0.39 — quality trimming of raw FASTQ files
- BWA v0.7.17 (BWA-MEM algorithm) — mapping trimmed reads to the Colias crocea reference genome (PRJEB42949)
- SAMtools v1.16.1 — removing PCR duplicates and downsampling to 4× coverage
- Picard v2.27.5 — adding read group identifiers
- NGSRelate (Korneliussen & Moltke 2015) — identifying related individuals
- GATK v4.3.0 — genotype calling (HaplotypeCaller, GenomicsDBImport, GenotypeGVCFs, VariantFiltration, SelectVariants)
- BCFtools v1.16 — marking missing values
- BEAGLE v4.1 — genotype imputation
- PLINK v1.90 — linkage disequilibrium pruning and allele frequency calculation
- ANGSD (Korneliussen et al. 2014) — genotype likelihood-based analyses and sensitivity analysis
- PCAngsd (Meisner & Albrechtsen 2018) — principal component analysis from genotype likelihoods
- BEDtools v2.29.1 — intersecting SNP positions with genome annotation
- SnpEff v5.2 — functional annotation of SNPs
- Pannzer2 — gene ontology annotation of the C. crocea reference genome
- GOWINDA v1.12 — gene set enrichment analysis
Statistical analyses:
- R v4.1.2 — all statistical analyses and data visualisation
- vegan — principal component analysis, Mantel test
- geosphere — geographic distance calculations (Haversine)
- elevatr — extraction of elevation data
- ggplot2 — data visualisation
- patchwork / ggpubr — figure assembly
- dplyr, tidyr, stringr, readr — data wrangling
Workflow: Raw FASTQ files were processed through the bioinformatic pipeline to produce the VCF files deposited in this repository. Downstream statistical analyses were conducted in R using the VCF files and associated environmental data. Full workflow details are provided in the methods section of Durkee et al. 2026.
Access information
Other publicly accessible locations of the data:
- Raw sequence data: NCBI SRA, accession PRJNA1433049
Data derived from the following sources:
- Colias crocea reference genome: European Nucleotide Archive, accession PRJEB42949 (Ebdon et al. 2022)
- Climate variables (precipitation, temperature, solar radiation, wind speed): WorldClim v2, 30 arc-second resolution, 1970–2000 (Fick & Hijmans 2017); https://www.worldclim.org
- Elevation data: extracted using the elevatr R package (Hollister et al. 2023) via the AWS Terrain Tiles API
- Gene Ontology annotations: assigned using Pannzer2 (Törönen et al. 2018); https://ekhidna2.biocenter.helsinki.fi/sanspanz/
