Data and code for Satokangas et al. 2026 Satokangas I, Martin SH, Seifert B, Puukko T, Schultz R, Helanterä H, Kulmuni J. "Introgression and divergence in a young species group". Molecular Ecology. Contact: ina.satokangas@ed.ac.uk | j.k.kulmuni@uva.nl

OVERVIEW

This archive contains data and code supporting the manuscript. Materials are organised into two deposits:

Dryad (this archive) — VCF files and analysis-specific input/output files Zenodo — Bioinformatic scripts (linked from this Dryad record)

Raw whole-genome sequencing reads and aligned reads (BAM files) are deposited at the European Nucleotide Archive (ENA) under project accession PRJEB112102. Outgroup F. exsecta reads: SRA accession SRS2348679.

All analyses were run on the CSC Puhti HPC cluster (Finland) using conda environments deployed via tykky containerization (see 000_conda_envs.txt in the Zenodo scripts archive).

ALIGNED READS (BAM FILES)

Deduplicated, overlap-clipped BAM files (mapped reads) for all individuals are available at the European Nucleotide Archive (ENA) under project accession PRJEB112102, alongside the raw sequencing reads. They are required to re-run the SNP calling pipeline from Step 6 onwards, by phyml_gene_trees.sh (WhatsHap read-backed phasing), and by coverage_100kb_windows.sh (mosdepth per-sample depth).

VCF FILES (root of this archive)

These files are provided individually as they are large and used as starting points for multiple analyses. All files are standard VCF format, block-gzipped (bgzip); open/query with BCFtools or VCFtools (both free, open source). Sample names in the VCF header follow the individual ID naming convention described under ANALYSIS-SPECIFIC INPUT AND OUTPUT FILES below.

highcovsamples_filtered_noindels_rmdup.vcf.gz All-sites VCF (invariant + variant sites) for 25 high-coverage individuals (5 per species). Produced by 8_allsites_vcf.sh. Used as input for pixy (π, dxy, Fst).

DP8.93inds.AN10.noScaff0003.mac2.vcf.gz Main SNP VCF: 93 F. rufa group individuals, 1,890,044 biallelic SNPs, DP≥8, ≤10% missing data, MAC≥2. Produced by 7_VCF_filtering.sh. Used for NJ tree, NeighbourNet, fd statistics, Dsuite.

DP8.93inds.AN10.noScaff0003.mac2.Fexs.SNP.gtFix.vcf.gz Main SNP VCF + F. exsecta outgroup (94 samples). Produced by 10_Outgroup_genotyping.sh. Used for Dsuite/Fbranch and as input to fd_statistics.sh (.geno conversion).

master.whatsHapClean.whap.shapeit.allScafs.outgroup.gtfix.nomiss.vcf.gz Phased VCF (WhatsHap + ShapeIt4) for 73 non-admixed ingroup samples + F. exsecta (74 total), filtered to sites with <50% missing data (1,881,405 SNPs). Produced by phyml_gene_trees.sh. Deposited to allow re-use of the phasing step (takes ~2 days on HPC).

ANALYSIS-SPECIFIC INPUT AND OUTPUT FILES

These files are in Analysis_specific_input_and_output.zip. Subfolder structure mirrors the scripts that produced them (see Scripts/README.txt in the Zenodo archive). All files are plain text unless noted otherwise, and can be opened, e.g., with any text editor. Files with a .gz extension are gzip-compressed and can be decompressed with gunzip or read directly by most software (e.g. R read.table()).

Individual ID naming convention:

IDs follow the pattern -F, e.g. 102-Frufa (F. rufa), 104-Faqu (F. aquilonia), 108-Flug (F. lugubris), 10-Fpol (F. polyctena), 117-Fprat (F. pratensis). A trailing "H" (e.g. 44-FaquH) denotes an individual morphologically classified as hybrid but should be ignored as these are inaccurate. Some samples use site-code IDs instead of numbers (e.g. CAGa_1w, CBCH1_1w). The outgroup F. exsecta individual is labelled "Fexs" or "Outgroup" depending on the file. See the manuscript's Supporting Information (S1 Table) for species, sampling location, and country for every individual ID.

Species abbreviation key (used throughout this folder):

aqu, aqu_fi F. aquilonia (aqu_fi = Finnish individuals specifically) aqu_ss F. aquilonia, Switzerland/Scotland subset lug F. lugubris lug_fi_1/_2 F. lugubris, Finland, (split into 2 subgroups used for the "split" Dsuite/fd analyses, see below) pol F. polyctena pol_eswi/_wswi F. polyctena, eastern/western Switzerland subgroups prat, prat_fi_1/_2 F. pratensis (Finland), split into 2 subgroups rufa, rufa_fi_1/_2 F. rufa (Finland), split into 2 subgroups Outgroup, Fexs F. exsecta (outgroup species) xxx individual not used in that particular file

Dsuite_Fbranch/ — genome-wide introgression (f-branch statistic)

Software: Dsuite v0.5 (Malinsky et al. 2021; github.com/millanek/Dsuite)

Script: Dsuite_Fbranch.sh

Input files (tab-separated, 2 columns, no header: individual ID, species/group label using the abbreviation key above):

indset_simple_189milSNPs.csv Grouping for the "simple" run (one 5-individual set per species)

indset_split_189milSNPs.csv Grouping for the "split" run (each species divided into 2 within-species subgroups of 2-3 individuals)

treehyp_simple_189milSNPs.txt Newick tree hypothesis matching the "simple" grouping, supplied to Fbranch

treehyp_split_189milSNPs.txt Newick tree hypothesis matching the "split" grouping, supplied to Fbranch

Output files:

indset_simple_189milSNPs_tree.txt and indset_split_189milSNPs_tree.txt\ Dsuite Dtrios output tables (one row per taxon trio); direct input to Dsuite Fbranch.

Tab-separated columns: P1, P2, P3 the three ingroup taxa tested (P3 = putative introgression source) Dstatistic Patterson's D-statistic (ABBA-BABA asymmetry) Z-score D-statistic in standard-deviation units from 0 (block-jackknife); |Z| > 3 taken as significant p-value two-tailed p-value corresponding to the Z-score f4-ratio estimated genome-wide admixture proportion BBAA, ABBA, BABA allele-sharing site-pattern counts among P1, P2, P3 and the outgroup used to compute D and f4-ratio

formica_Fbranch_simple_189milSNPs.txt and formica_Fbranch_split_189milSNPs.txt\ Fbranch (fb) statistic matrices. Tab-separated; "branch" = internal/terminal branch of the tree hypothesis; "branch_descendants" = taxa descending from that branch; remaining columns = each taxon in the tree. Cell value = fb, the estimated proportion of introgression from the taxon named in the column header into the ancestor of the branch named in that row. "nan" = comparison not possible given the tree topology.

formica_Fbranch_simple_189milSNPs_zscore.txt Z-scores for each fb value

formica_Fbranch_split_189milSNPs_zscore.txt in the matching table above (same row/column layout); "nan" = not tested/not possible.

fd_statistics/ — windowed introgression (fd statistic)

Software: genomics_general scripts ABBABABAwindows.py and parseVCF.py

(S. H. Martin; github.com/simonhmartin/genomics_general)

Script: fd_statistics.sh

Input files:

abwind.simple.pops.txt Population grouping file (tab-separated, 2 columns, no header: individual ID, species label) for "simple" (5 individuals per species) trios.

abwind.split.pops.txt Population grouping file for "split" (2-3 individuals per within-species subgroup) trios, using the subgroup labels in the abbreviation key above.

ploidy.file.tab Tab-separated, 2 columns, no header: individual ID; ploidy (2 = diploid ingroup individual, 1 = haploid outgroup "Fexs"). Required by parseVCF.py for correct allele-frequency conversion (.geno file).

Output files (100kb and 50kb window sizes):

File names give the three taxa tested as P1.P2.P3 (ABBA-BABA/Dsuite convention; see manuscript Table 1 for the full rationale for each trio):

lug.aqufi.pol_ABBABABA_100kb_out.csv lug.aqufi.pol_ABBABABA_50kb_out.csv lug.aqufi.prat_ABBABABA_100kb_out.csv lug.aqufi.prat_ABBABABA_50kb_out.csv lug.pol.prat_ABBABABA_100kb_out.csv lug.pol.prat_ABBABABA_50kb_out.csv pol.rufa.lug_ABBABABA_100kb_out.csv pol.rufa.lug_ABBABABA_50kb_out.csv rufa.pol.aqufi_ABBABABA_100kb_out.csv rufa.pol.aqufi_ABBABABA_50kb_out.csv rufafi2.rufafi1.polwswi_ABBABABA_100kb_out.csv rufafi2.rufafi1.polwswi_ABBABABA_50kb_out.csv (lug.aqufi.prat and lug.pol.prat are control trios with no introgression detected genome-wide by f-branch; the remaining four trios are those where introgression was detected, see Dsuite_Fbranch/ above)

Comma-separated columns in each output file, one row per genomic window: scaffold reference genome scaffold name start, end window start/end coordinates (bp, 1-based) mid window midpoint coordinate (bp) sites number of biallelic sites in the window sitesUsed number of sites used after missing-data filtering ABBA, BABA frequency-weighted sums of the ABBA and BABA site patterns across sites in the window D windowed Patterson's D-statistic, (ABBA-BABA)/(ABBA+BABA) fd windowed estimator of introgression proportion (Martin et al. 2015) fdM modified fd (Malinsky et al., as implemented in genomics_general)

NJtree/ — whole-genome neighbour-joining species tree (Fig 1)

Software: vcf2phylip v2.8 (Ortiz 2019) to build the alignment; R package phangorn (Schliep 2010) to build the tree; FigTree v1.4.4 for visualisation.

Script: nj_tree.sh

Files:

• .phy Phylip-format alignment (73 non-admixed F. rufa group individuals + F. exsecta outgroup; main SNP VCF thinned to a minimum 1kb spacing between SNPs, individuals with >1% missing data excluded; 165,332 SNPs). Header line gives the number of taxa and alignment length ("74 165332"); each following line is one individual ID followed by its IUPAC-coded genotype sequence. Input to the R/phangorn NJ tree.

• .tree Output NJ tree with 100 bootstrap replicates, in NEXUS format (readable in FigTree, R ape/phangorn, etc.), rooted with F. exsecta ("Fexs"). Corresponds to manuscript Fig 1. Tip labels are individual IDs (see manuscript S1 Table for species and geographic origin of each ID).

Neighbournet/ — phylogenetic network (S1 Fig)

Software: distMat.py (S. H. Martin, genomics_general) to compute distances; SplitsTree v4.17.1 (Huson & Bryant 2006), NeighbourNet method (Bryant & Moulton 2004), to build and view the network.

Script: neighbournet.sh

Files:

DP8.93inds.AN10.noScaff0003.mac2.thin20kb_numonly.dist Pairwise genetic distance matrix for all 93 F. rufa group individuals (main SNP VCF thinned to a minimum 20kb spacing between SNPs; 9,816 SNPs). Plain-text, Phylip-style distance-matrix format: first line gives the number of individuals (93); each following line gives one individual's numeric ID followed by its genetic distances to all other individuals, in the same order. Import into SplitsTree (Distances -> Import) and apply NeighbourNet with default parameters to reproduce the network (S1 Fig). Individual IDs are the numeric part of the sample names used elsewhere in this archive (species and geographic origin given in manuscript Fig 1 and Supporting Information S1 Table).

Gene_trees/ — windowed gene trees and topology weighting (TWISST, S3 Fig)

Software: phyml_sliding_windows.py, which calls PhyML internally (S. H. Martin, genomics_general), for gene trees; TWISST (Martin & Van Belleghem 2017; github.com/simonhmartin/twisst) for topology weighting.

Scripts: phyml_gene_trees.sh (trees), twisst.sh (weights)

Files:

groupfile_genetrees.tab Tab-separated, 2 columns, no header: haplotype ID (individual ID with an "_A" or "_B" suffix for the two phased haplotypes; the haploid F. exsecta outgroup has no suffix); species/taxon label defining the 6 TWISST groups (prat, aqu, lug, rufa, pol, and the F. exsecta outgroup).

minsites20.minsnps20.HKY85.noOpt.phyml_bionj.trees.gz and minsites50.minsnps50.HKY85.noOpt.phyml_bionj.trees.gz

Gzipped Newick gene trees, one per genomic window, built with PhyML under the HKY85 substitution model (no parameter optimisation, "noOpt") from coding-sequence windows with a minimum of 20 or 50 variant sites respectively. Readable with any Newick-compatible parser (e.g. R ape, dendropy, ete3).

minsites20.minsnps20.HKY85.noOpt.phyml_bionj.weights.csv.gz minsites50.minsnps50.HKY85.noOpt.phyml_bionj.weights.csv.gz

Gzipped CSVs, one row per gene tree (same order as the matching .trees.gz file). One column per possible unrooted topology among the 5 ingroup taxa; values are TWISST topology weights (relative frequency, 0-1, of that topology within the window; weights in a row sum to 1).

minsites20.minsnps20.HKY85.noOpt.phyml_bionj.data.tsv minsites50.minsnps50.HKY85.noOpt.phyml_bionj.data.tsv

Tab-separated per-window metadata, row order matching the corresponding trees/weights files.

Columns: scaffold, start, end genomic coordinates (bp) of the window mid window midpoint coordinate (bp) sites number of variant sites used to build the tree lnL PhyML tree log-likelihood ("NA" throughout, as no model parameters were optimised, "noOpt")

PIXY/ — diversity, divergence, differentiation in genomic windows (Fig 2, 3)

Software: pixy v1.2.7 (Korunes & Samuk 2021)

Script: pixy.sh

Input file:

group_pixy_100kb_5inds_highcovsamples.tab Tab-separated, 2 columns, no header: individual ID; species label (aqu, lug, pol, prat, rufa) for the 5 high-coverage individuals per species (25 individuals total) used for pixy; the same individuals used for mosdepth coverage calculation (see Spearman_and_chi2/ below).

Output files (100kb and 50kb non-overlapping windows; units are per-site averages unless stated otherwise):

pixy_100kb_highcov_5inds_dxy.txt pixy_50kb_highcov_5inds_dxy.txt pop1, pop2 the two species compared chromosome scaffold name window_pos_1, window_pos_2 window start/end coordinates (bp) avg_dxy average pairwise nucleotide divergence between pop1 and pop2, per site no_sites sites with data in the window count_diffs pairwise sequence differences, summed over all between-population comparisons count_comparisons total number of pairwise comparisons count_missing comparisons with missing genotype data

pixy_100kb_highcov_5inds_fst.txt pixy_50kb_highcov_5inds_fst.txt pop1, pop2, chromosome, window_pos_1, window_pos_2 as above avg_wc_fst Weir & Cockerham (1984) Fst estimate for the window (negative values truncated to zero) no_snps number of SNPs used in the Fst calculation

pixy_100kb_highcov_5inds_pi.txt pixy_50kb_highcov_5inds_pi.txt pop species chromosome, window_pos_1, window_pos_2 as above avg_pi average within-species nucleotide diversity per site no_sites, count_diffs, count_comparisons, count_missing as above, but for within-population comparisons

Spearman_and_chi2/ — coverage filtering and genome-wide correlation tests

Software: mosdepth v0.3.3 (Pedersen & Quinlan 2018); R (base stats package: cor.test for Spearman's rank correlation, chisq.test for Pearson's chi-squared test)

Scripts: coverage_100kb_windows.sh (mosdepth), Spearman_and_chi2.R

Files:

mosdepth_highcovinds_100kb_regions_depth_coreregions.txt Tab-separated, no header row, one row per 100kb genomic window.

Column 1: scaffold; column 2: window start (bp, 0-based); column 3: window end (bp); columns 4-28: mean sequencing depth (X, from BAM files) for each of the 25 high-coverage individuals, in the same order as the 25 individuals listed in ind_cov_species.list below. Produced by coverage_100kb_windows.sh. Used to identify and exclude genomic windows where any species' coverage exceeded twice its mean (likely collapsed/duplicated regions) from the genome-wide correlation analyses (Spearman_and_chi2.R).

ind_cov_species.list Tab-separated, 2 columns. The first 3 rows ("chr chr", "start start", "end end") are column labels for the 3 coordinate columns of the depth table above; the remaining 25 rows give individual ID and species label (aqu, lug, pol, prat, rufa) for the same 25 high-coverage individuals used in the depth table and in the pixy analysis, in matching order. Note: the recombination rate file (3A_windows_20kb_loter_CDS_rho.tsv) used by Spearman_and_chi2.R is not re-deposited here; download from Nouhaud et al. 2022b figshare: https://figshare.com/articles/dataset/21511008

SCRIPTS (Zenodo)

See Scripts/README.txt in the Zenodo archive for the full annotated list.
Scripts are organised into two subfolders:
SNP_calling_pipeline_scripts/ — Steps 1–10: raw reads → filtered VCF Other_analyses_scripts/ — Downstream population genomic analyses