This dataset consists of (1) sequence data processing for all of the downstream analysis, (2-3) population genomic analyses, (4) linkage map constructions for the downstream analysis, (5) screening of the genomic loci associated with survival, (6) QTL mapping for habitat choice, (7) permutation test for barrier effect overlap with HMM States, (8) estimation of selection coefficients, and (9) individual-based simulations.

Description of the data and file structure

Detailed information on data collection and analytical methods is provided in the associated manuscripts.

1. Sequencing Data Processing

Overview

This section contains scripts for processing sequencing data and variant calling used in the manuscript. The first script (01_1_fastq_trimming_filtering_mapping.sh) trims and filters FASTQ files with fastp, aligns reads to a reference genome using BWA, and processes alignments with SAMtools. The second script (01_2_SNP_call_separate_chromosome.sh) performs variant calling using bcftools mpileup and bcftools call, generating gzipped VCF files per chromosome. The third script (01_3_Concat_filtering_Liftover.sh) concatenates VCF files, applies filtering, performs liftover to a new genome assembly, and prepares VCFs for downstream analyses such as ADMIXTURE and introgression.

• 01_1_fastq_trimming_filtering_mapping.sh: FASTQ Trimming, Filtering and Mapping
• 01_2_SNP_call_separate_chromosome.sh: SNP call
• 01_3_Concat_filtering_Liftover.sh: VCF Processing

Requirements

• fastp v 0.20.1 or higher
• bwa-mem v 0.7.17 or higher
• vcftools v 0.1.15 or higher
• GATK 4.2.0.0 or higher

Usage

Run 01_1_fastq_trimming_filtering_mapping.sh to:

• Trim and filter FASTQ files using fastp
• Align reads to the reference genome with BWA
• Sort and index BAM files using SAMtools. For paired-end libraries, Read 1 was trimmed to retain the terminal 100 bp, and only the first 50 bp of the resulting read were used to match the read length of the single-end 50 bp dataset.

Create a BAM path list file in your directory manually:

• Set path to ./path/to/bamlist.txt
• Run 01_2_SNP_call_separate_chromosome.sh to:
• Perform SNP calling per chromosome using bcftools mpileup and bcftools call
• Generate gzipped VCF files

Run 01_3_Concat_filtering_Liftover.sh to:

• Concatenate per-chromosome VCF files
• Apply filtering (remove indels, set missingness and quality thresholds)
• Perform liftover to a new genome assembly using GATK
• Prepare filtered VCF files for ADMIXTURE analysis, introgression analyses, and QTL mapping.

Notes

Modify scripts according to your cluster's job scheduler. For performing Hosoki et al., download fastq from the Sequence Read Archive (DRA accession numbers: PRJDB39618, PRJDB39649, PRJDB38001, PRJDB38013, PRJDB38020, PRJDB42376). Modify paths in scripts according to your directory structure.

2. Admixture and Introgression Analysis

Overview

This section contains scripts for population structure and introgression analyses. The first script (02_1_ADMIXTURE.sh) converts VCF files to PLINK format, performs LD pruning, and then runs ADMIXTURE analysis. The second script (02_2_introgress_run.R) is an R script that estimates hybrid index and heterozygosity using the introgress package and produces the triangle plots.

Scripts:

• 02_1_ADMIXTURE.sh: Admixture Analysis: VCF to PLINK, LD pruning, ADMIXTURE runs
• 02_2_Introgress_run.R: Introgression Analysis in R: Hybrid index, genomic clines, triangle plots

Requirements

• PLINK v1.9 or higher
• ADMIXTURE v1.3.0 or higher
• R (≥ 4.0) with packages: introgress, dplyr, ggplot2, ggdensity
• Python 3.x and vcf2introgress.py for input preparation:
• https://github.com/pimbongaerts/radseq/blob/master/vcf2introgress.py

Usage

Run 02_1_ADMIXTURE.sh to:

• Get path to the VCF from 01_3_Concat_filtering_Liftover.sh
• Convert VCF to PLINK format
• Perform LD pruning
• Run ADMIXTURE for K=1 to 30 (K=2 was used in Hosoki et al.)

Prepare introgress input:

• Get path to the VCF from 01_3_Concat_filtering_Liftover.sh
• Run the vcf2introgress.py to convert VCF into introgress-compatible files

Run 02_2_introgress_run.R in R to:

• Estimate hybrid index and heterozygosity
• Generate triangle plots by sampling year

Notes

Use the vcf files converted by the 01_3_Concat_filtering_Liftover.sh above. Modify paths in scripts according to your directory structure.

3. Define Species-Diagnostic Loci and Wright-Fisher_Simulation

Overview

This section contains scripts for extracting species-diagnostic loci, estimating allele frequencies, and simulating allele frequency changes under the Wright-Fisher model. The first script (03_1_Allele_frequency_estimation_at_species-diagnostic_loci.sh) identifies loci differentiated between parental species. From the loci definitions, the first script calculates allele frequencies by pooled groups. The second script (03_2_Wright-Fisher Simulation.R) runs Wright-Fisher simulations to compare observed vs. expected allele frequency changes.

Scripts:

• 03_1_Allele_frequency_estimation_at_species-diagnostic_loci.sh: Extract species-diagnostic loci and calculate allele frequencies.
• 03_2_Wright-Fisher_Simulation.R: Simulate allele frequency changes.

Requirements

• Python 3.x and vcf2introgress.py for input preparation
• https://github.com/pimbongaerts/radseq/blob/master/vcf2introgress.py
• vcftools v0.1.15 or higher
• PLINK v1.9 or higher
• R (≥ 4.0) with packages: dplyr, jmuOutlier

Usage

Prepare input data:

• Run 01_3_Concat_filtering_Liftover.sh to generate SNP data.
• Run 03_1_Allele_frequency_estimation_at_species-diagnostic_loci.sh to:
• Identify species-diagnostic loci using vcf2introgress.py
• Extract differentiated SNPs from hybrid VCF files
• Convert filtered VCF to PLINK format
• Calculate allele frequencies using PLINK

Run 03_2_Wright-Fisher Simulation.R to:

• Simulate allele frequency changes using Wright-Fisher model
• Compare observed vs. expected allele frequencies for species-diagnostic loci
• Perform permutation tests to assess observed changes vs Wright-Fisher changes

Population sizes tested: 100% (22485), 10% (2249), 1% (225) census size

Notes

Use the vcf files converted by the 01_3_Concat_filtering_Liftover.sh above. For Wright-Fisher Simulation, use the results from the 03_1_Allele_frequency_estimation_at_species-diagnostic_loci.sh. Modify paths in scripts according to your directory structure.

4. Linkage map construction

Overview

This section contains scripts for constructing linkage maps for identifying loci associated with survival in a tsunami-created habitat and QTL mapping of habitat choice below. The script (04_make_linkagemap.sh) builds linkage maps using Lep-MAP3.

Scripts:

04_make_linkagemap.sh: Construct linkage map using Lep-MAP3 (ParentCall2, Filtering2, SeparateChromosomes2, OrderMarkers2, genotype conversion).

Requirements

Lep-MAP3

Usage

Prepare input data:

Run 01_3_Concat_filtering_Liftover.sh to generate SNP data.

Run 04_make_linkagemap.sh to:

• Call parental genotypes using ParentCall2
• Filter markers using Filtering2
• Assign linkage groups using SeparateChromosomes2
• Order markers within linkage groups using OrderMarkers2
• Convert phased data to genotype format using map2genotypes.awk

Notes

Modify paths in scripts according to your directory structure. For performing Hosoki et al. inprep, download fastq from Sequence Read Archive (DRA accession numbers PRJDB39618, PRJDB39649, PRJDB38001, PRJDB38013, PRJDB38020, PRJDB42376)

5. Identifying loci associated with survival in a tsunami-created habitat

Overview

This section contains scripts for identifying loci associated with survival in a tsunami-created habitat based on allele frequency changes before release and after recapture. The first script (05_1_FWenc_AFD_fileconversion.R) converts Lep-MAP3 genotype output for downstream analyses. The second script (05_2_FWenc_AF_calculation.sh) calculates allele frequencies for experimental groups (before and after). Before running this script, construct a linkage map using 04_make_linkagemap.sh as described above.

Scripts:

• 05_1_FWenc_AFD_fileconversion.txt: Convert Lep-MAP3 genotype output to an input to downstream analysis and add metadata columns and parental genotypes.
• 05_2_FWenc_AF_calculation.sh: Calculate allele frequencies for FWenc experiment and compute allele frequency differences (AFD) between pre-release and recapture groups.

Requirements

• R (≥ 4.0) with packages: dplyr, tidyr
• vcftools v0.1.15 or higher
• PLINK v1.9 or higher
• Python 3.x and vcf2introgress.py for input preparation
• https://github.com/pimbongaerts/radseq/blob/master/vcf2introgress.py

Usage

Prepare linkage map and input data:

Run 04_1_make_linkagemap.sh to generate linkage map and genotype data.

Run 05_1_FWenc_AFD_fileconversion.R to:

Convert Lep-MAP3 genotype output to downstream analysis

Run 05_2_FWenc_AF_calculation.sh to:

• Identify parental informative markers using vcf2introgress.py
• Calculate allele frequencies for each experimental group using PLINK

Notes

Adjust sex chromosome allele frequencies manually before allele frequency difference calculation. Modify paths according to your directory structure. For performing Hosoki et al., download fastq from Sequence Read Archive (DRA accession number PRJDB39618, PRJDB39649, PRJDB38001, PRJDB38013, PRJDB38020, PRJDB42376).

6. QTL mapping of habitat choice

Overview

This script performs QTL mapping using the R/qtl package. It includes data import, genotype probability calculation, and interval mapping using Haley-Knott regression. Before running this script, construct a linkage map using 04_make_linkagemap.sh as described above.

Scripts:

06_QTL_mapping.R: QTL mapping

Requirements

R (≥ 4.0)

R packages:

• r/qtl
• ggplot2

Usage

Prepare linkage map and input data:

• Use 04_1_make_linkagemap.sh to generate linkage map and genotype data.
• Prepare phenotype-genotype combined CSV file for R/qtl.

Load data and check missing values:

Import phenotype-genotype data using read.cross()

Simulate genotypes and calculate probabilities:

• Use sim.geno() for genotype simulation
• Use calc.genoprob() for genotype probabilities

Run interval mapping:

• Perform Haley-Knott regression with scanone()
• Permutation test for significance threshold
• Run scanone() with n.perm = 1000
• Determine LOD threshold and identify significant QTLs

Notes

Modify paths in scripts according to your directory structure. For performing Hosoki et al., download fastq from Sequence Read Archive (DRA accession number PRJDB39618, PRJDB39649, PRJDB38001, PRJDB38013, PRJDB38020, PRJDB42376).

7. Permutation test for barrier effect overlap with HMM States

Overview

This script performs permutation tests to evaluate whether genomic regions associated with barrier effects significantly overlap with HMM state assignments across multiple generations. Before running, compute HMM states using HMM_log10FST+1_3norm.R from below:
https://github.com/marqueda/HMM-detection-of-genomic-islands/tree/master.

Scripts:

07_PermutationTest_for_Barrier_Effect_Overlap_with_HMM_States.R: test association between barrier loci and HMM outliers

Requirements

R (≥ 4.0)
HMM states using HMM_log10FST+1_3norm.R from below:

https://github.com/marqueda/HMM-detection-of-genomic-islands/tree/master.

R packages:

dplyr
ggplot2

Usage

• Run HMM
• Calculate HMM states from the second to third generation of the tsunami-created population.
• Calculate HMM states from two years to nine years after the tsunami
• Load input data
• Import position of barrier loci
• Run permutation tests
• For each generation and HMM state, shuffle barrier effect labels and compute overlap counts.
• Repeat for 100,000 iterations to build null distributions.
• Calculate empirical p-values for observed overlaps.

Note

Use the observed allele frequencies generated from the section ”3. Species-Diagnostic Loci and Wright-Fisher_Simulation”. Barrier loci data were from Kitano et al. 2009, Nature, “4. Linkage map construction and identifying loci associated with survival in a tsunami-created habitat”, and “5. QTL mapping of habitat choice”

8. Bayesian_etimation_submission

Overview

This repository contains Python scripts used for Bayesian estimation of selection coefficient in the manuscript. The main script (Bayesian_estimation_fitness.py) processes an input data file in CSV format and generates summary figures and a results table.

Contents

• 08_Bayesian_estimation_fitness.py: Main analysis and plotting script
• 08_genotype_gen_ligestage.csv: Input dataset used by the script

Requirements

Python 3.8 or higher
Required Python packages: pandas, numpy, matplotlib, arviz, pymc
(Install via: pip install pandas numpy matplotlib arviz pymc)

Usage

Run the main script from the command line as follows:

python 08_Bayesian_estimation_fitness.py

The script reads 08_genotype_gen_lifestage.csv () from the current directory and produces:

• Three PNG figures (Fig_Posterior_p_first.png, Fig_Posterior_sc.png, Fig_Posterior_tc.png)
• One results table (Result_Bayesian_estimation_fitness_Hosoki_et_al.csv)
• These output files are saved in the same directory as the script.

Notes

The 08_Bayesian_estimation_fitness.py
reads 08_genotype_gen_lifestage.csv from the current directory. This CSV file contains the following variables: stage (generation and life stage), CHR (chromosome identifier), locus (locus identifier), HOM.Gnip (number of G. nipponicus homozygotes), HET (number of heterozygotes), HOM.Gacu (number of G. aculeatus homozygotes), barrier_effect (reproductive barrier category assigned to the locus), and missing (proportion of missing genotype data, 0–1).

Fig_Posterior_XXX.png gives the plot_trace output of arviz for the posterior distribution of the estimated parameter XXX.

9. Individual-based simulation

Overview

The directory code_IBM contains the files used for individual-based simulations. It includes five source code files written in the Julia programming language, two text files, and two subdirectories. The Julia programs are used for running simulations, analyzing the results, and visualizing them. The text files serve as input files for the Julia programs. The subdirectory appendix_inputs_Figs contains the input files used to reproduce the simulations shown in the figures of the paper (Figs. S16, S18, and S19). The subdirectory appendix_outcomes_Figs contains the final output files from the simulations shown in the figures of the paper (Figs. S16, S18, and S19).

Contents

09_code_IBM.zip:

Package_Installation.jl: Installs all packages required to run the other programs.

Model.jl: Defines the core functions of the simulation model.

Simulation.jl: Contains functions for running simulations of the model under varied conditions.

Analysis.jl: Provides functions for analyzing and visualizing results from multiple simulations conducted under different conditions.

Run_Simulations.jl: Executes simulations using multiple threads.

parameter_defo.txt: Specifies the default values for all simulation parameters (see Table 1 in the paper).

parameter_list.txt: Lists alternative parameter values for the simulations.

Table.1.csv:

Parameters for the simulation program and their symbols in the main paper

Instructions for running the Julia programs

Running simulations with various parameter conditions

The Julia programs can be executed using the Julia executable (julia command), which is available at the official website of the Julia programming language (https://julialang.org/). Since the programs were developed using Julia version 1.7.2, this version of Julia is required to reproduce the simulations.

Installation of the required packages, simulations, and the subsequent analyses can be performed by executing the following command in the directory IBM using Julia version 1.7.2:

julia --project=. -p2 Run_Simulations.jl ./ dir_name

The argument -p2 specifies the number of threads used for parallel computation. The two command-line arguments (./ and dir_name) specify where to save the output files:

• The first argument (./) designates an existing directory in which a new subdirectory will be created to store the output files.
• The second argument (dir_name) specifies the name of this new subdirectory.

The programs read two input files: parameter_defo.txt and parameter_list.txt. The former defines the default parameter values (Table 1), and the latter lists alternative parameter values. Simulations are automatically run for all parameter combinations specified in parameter_list.txt. Each simulation produces raw output files, which are then automatically analyzed and integrated into a single summary file (Outcomes.txt). Additionally, the program generates figures visualizing the simulation results. Note: On Linux, the program to generate figures may cause the following error message: “Qt: Session management error: None of the authentication protocols specified are supported”. Although the figures will be generated even with this error, the error will be avoided by running the following command before calling the Julia program: unset SESSION_MANAGER

Simulations with any parameter conditions can be performed by editing parameter_list.txt. The subdirectory appendix_inputs_Figs includes seven example files of parameter_list.txt, with which simulations under the same conditions as those shown in figures of the paper can be conducted (Figs. S16, S18, S19). Note: Running all 13,680 simulations listed in these files requires a very long computation time.

Format of the final output file

Each row of the final output file (Outcomes.txt) represents the parameter settings and results of a single simulation. The first 15 columns record the parameter values used in that simulation. The remaining two columns, time_purging and obs_reproduced, summarize the simulation outcomes:

time_purging: The time until the genome-wide purging of G. nipponicus alleles (i.e., g_purging as defined in the Materials and Methods section of the paper).

obs_reproduced (true/false): Indicates whether the simulation reproduced the observed pattern (see Materials and Methods for details).

11. Gene Richness Window Analysis

Overview

This script extracts structural gene regions from a genomic annotation file (GFF3), generates sliding genomic windows at multiple scales (10kb and 100kb), and quantifies localized gene density using BEDtools. Downstream analysis integrates these structural density matrices with temporal allele frequency datasets to evaluate correlation via Spearman's rank coefficients across generations.

Scripts:

• 010_generichness_count.sh: Extracts features, generates sliding windows, profiles regional gene counts, and documentation for subsequent R integration.

Requirements

• Command-line Tools: BEDtools (v2.30+), BEDOPS (gff2bed), Samtools
• R (> 4.0) packages: dplyr (>1.1.0 for advanced join_by functions)

Usage

Prepare genomic features and reference index: Isolate explicit coordinate elements flagged as 'gene' from the core GFF3 file and convert to BED format via gff2bed. Generate a valid reference index file (.fai) using samtools faidx. Generate genomic windows and calculate gene richness: Use bedtools makewindows to build non-overlapping sliding window tracks (10kb and 100kb) across the genome. Use bedtools intersect with the -c option to quantify total gene presence counts inside each window. Map locus positions and calculate mean frequencies in R: Import the generated window-based gene count file and allele frequency data. Perform a concatenation using dplyr::join_by() to securely anchor each marker locus into window boundaries. Calculate the mean allele frequency (A1_freq) within each genomic window using group_by() and summarise(). Run correlation analysis: Evaluate the correlation between gene richness and allele frequency changes for each year using cor.test() with method = "spearman".

Note

Modify paths in scripts according to your directory structure.