DOI: 10.5061/dryad.7wm37pw09
title: Genomic landscape of introgression from the ghost lineage in a gobiid fish uncovers the generality of forces shaping hybrid genomes

description: Extinct lineages can leave legacies in the genomes of extant lineages
 through ancient introgressive hybridization. The patterns of genomic
 survival of these extinct lineages provide insight into the role of
 extinct lineages in current biodiversity. However, our understanding of
 the genomic landscape of introgression from extinct lineages remains
 limited due to challenges associated with locating the traces of unsampled
 “ghost” extinct lineages without ancient genomes. Herein, we conducted
 population genomic analyses on the East China Sea (ECS) lineage of
 Chaenogobius annularis, which was suspected to have originated from ghost
 introgression, with the aim of elucidating its genomic origins and
 characterizing its landscape of introgression. By combining
 phylogeographic analysis and demographic modeling, we demonstrated that
 the ECS lineage originated from ancient hybridization with an extinct
 ghost lineage. Forward simulations based on the estimated demography
 indicated that the statistic γ of the HyDe analysis can be used to
 distinguish the differences in local introgression rates in our data.
 Consistent with introgression between extant organisms, we found reduced
 introgression from extinct lineage in regions with low-recombination rates
 and with functional importance, thereby suggesting a role of linked
 selection that has eliminated the extinct lineage in shaping the hybrid
 genome. Moreover, we identified enrichment of repetitive elements in
 regions associated with ghost introgression, which was hitherto
 little-known but was also observed in the reanalysis of published data on
 introgression between extant organisms. Overall, our findings underscore
 the unexpected similarities in the characteristics of introgression
 landscapes across different taxa, even in cases of ghost introgression.
description: We performed population genomics analysis based on RAD sequence
 data (279 individuals from 21 populations, and one outgroup individual)
 and whole genome resequencing data (18 individuals from three populations
 and one outgroup individual) for <em>Chaenogobius
 annularis</em>, a coastal goby species inhabiting the Japanese
 archipelago. Our analyses included repeat annotation, gene annotation,
 population recombination rate estimation, estimation of potential
 deleterious mutations, population genetic analysis, test for
 hybridization, PSMC analysis, fastsimcoal2-based demographic modeling,
 SLiM-based forward simulation, sliding window analysis, and
 permutation-based characterization of introgression landscapes.
description: # Genomic landscape of introgression from the ghost lineage in a gobiid
 fish uncovers the generality of forces shaping hybrid genomes
 [https://doi.org/10.5061/dryad.7wm37pw09](https://doi.org/10.5061/dryad.7wm37pw09) ## **Brief description of the data and file structure** ### scripts.tar.gz #### **Note:** **"scripts.tar.gz" can be obtained from the Zenodo link (https://doi.org/10.5281/zenodo.10048869) tied to this Dryad page** (see "Related works" section in the upper right corner). **To get the script only, please visit this Zenodo link**. #### For your convenience, we have changed the "scripts.tar.gz" file to be available directly from Data files in this page as well (added on 12/18/2023). The scripts used in this study (bash, python, R). These scripts are categorized into the following 10 contents, which are hierarchized within each directory. 1\. ddRAD\-seq genotyping 2\. WGS \(whole genome resequencing\) genotyping 3\. repeats and gene annotation 4\. population recombination rate estimation 5\. potentially deleterious SNPs 6\. popultion genetic analyses 7\. phylogenetic analysis 8\. hybrid detection 9\. demographic estimation 10\. introgression landscape characterization The detailed hierarchical structure is given below in the section "*Detailed description of the file structure"*. ### genotyping\_data.tar.gz The compressed files of directories containing genotyping data generated in this study. Five VCF files from RAD-seq, two VCF files from whole genome resequencing data, and one FASTA file of whole mitogenomic sequences. The scripts used to analyze the demographic modeling are stored in "/01ddRADseq\_genotyping/" or "02WGS\_genotyping" in the scripts.tar.gz. ### annotation\_data.tar.gz The compressed files of directories containing repeat annotation data and gene annotation data in this study. The scripts used to analyze the demographic modeling are stored in "/03repeats\_and\_gene\_annotation/" in the scripts.tar.gz. ### demographic\_modeling.tar.gz The compressed file of a directory containing the results of the demographic modeling (distribution of AIC for each model, and maximum likelihood parameters for the best model) and the input site frequency spectrum. The scripts used to analyze the demographic modeling are stored in "/09demographic\_estimation/02demographic\_modeling/" in the scripts.tar.gz. ### slidingwindow\_results.tar.gz The compressed file of a directory containing the results of the sliding window anlysis (bed files summarizing the statistic *γ* in the HyDe analysis and some other features). Please see "README\_description\_of\_record\_XXkb.txt" in the slidingwindow\_results.tar.gz for column name descriptions. The scripts used for this analysis are stored in "/10introgression\_landscape\_characterization/02sliding\_window/" in the scripts.tar.gz. <br> ## **Detailed description of the file structure** ### "scripts.tar.gz" 1. **01ddRADseq\_genotyping** <br> Scripts used for genotyping double digest restriction-site associated DNA (ddRAD-seq) data 1. 01trimmomatic <br> Read filtering with Trimmomatic 2. 02mapping <br> Mapping with bwa-mem 3. 03sambamba <br> Filtering of BAM files by sambamba (retaining only uniquely mapped reads) 4. 04SNPcall\_by\_mpileup <br> Genotyping with samtools mpileup / bcftools call 1. 01pop1\_example <br> Example of a script for joint calling when targeting one group 5. 05merge\_filtering <br> Merging and filtering of VCF files to create the final genotype datasets 1. 01dataset1\_2 <br> Merging and filtering scripts for datasets 1 and 2 1. 01filtering\_dataset1 <br> Filtering scripts for dataet 1 2. 02filtering\_dataset2 <br> Filtering scripts for dataet 2 2. 02dataset3 <br> Merging and filtering scripts for dataset 3 1. 01filtering <br> Filtering scripts for dataet 3 3. 03dataset4 <br> Merging and filtering scripts for dataset 4 1. 01filtering <br> Filtering scripts for dataet 4 4. 04dataset5 <br> Merging and filtering scripts for dataset 5 1. 01filtering <br> Filtering scripts for dataet 5 2. **02WGS\_genotyping** <br> Scripts used for genotyping whole genome resequencing (WGS) data 1. 01nuclear\_SNP <br> Genotyping whole nuclear SNPs 1. 01mapping <br> Mapping with bwa-mem 2. 02sambamba <br> Filtering of BAM files by sambamba (retaining only uniquely mapped reads) 3. 03MarkDuplicate <br> Removing PCR duplicates in BAMfiles by MarkDuplicate (GATK4) 4. 04SNPcall\_by\_GATK <br> Genotyping with GATK4 1. 01SJ <br> Genotyping for the Sea of Japan (SJ) lineage 1. 01HaplotypeCaller <br> Script for GATK HaplotypeCaller 2. 02GenomicsDBImport\_GenotypeGVCFs <br> Script for GATK GenomicsDBImport and GenotypeGVCFs 2. 02ECS <br> Genotyping for the East China Sea (ECS) lineage 1. 01HaplotypeCaller <br> Script for GATK HaplotypeCaller 2. 02GenomicsDBImport\_GenotypeGVCFs <br> Script for GATK GenomicsDBImport and GenotypeGVCFs 3. 03PO <br> Genotyping for the Pacific Ocean (PO) lineage 1. 01HaplotypeCaller <br> Script for GATK HaplotypeCaller 2. 02GenomicsDBImport\_GenotypeGVCFs <br> Script for GATK GenomicsDBImport and GenotypeGVCFs 4. 04outgroup <br> Genotyping for the outgroup individual 1. 01HaplotypeCaller <br> Script for GATK HaplotypeCaller 2. 02GenomicsDBImport\_GenotypeGVCFs <br> Script for GATK GenomicsDBImport and GenotypeGVCFs 5. 05merge\_filtering <br> Merging and filtering of VCF files to create the final genotype datasets 1. 01Cannularis\_only\_dataset <br> Merging and filtering scripts for *C. annularis* only dataset 2. 02with\_outgroup\_dataset <br> Merging and filtering scripts for with-outgroup dataset 5. 05indelrealign <br> Local realignment around indels for BAM files with GATK3 6. 06phasing <br> Read-aware-phasing with shape-it2 2. 02mitogenome\_sequence <br> Genotyping whole mitogenome sequence 1. 01mapping <br> Mapping with bwa-mem 2. 02sambamba <br> Filtering of BAM files by sambamba (retaining only uniquely mapped reads) 3. 03MarkDuplicates <br> Removing PCR duplicates in BAMfiles by MarkDuplicate (GATK4) 4. 04SNPcall\_by\_GATK <br> Genotyping with GATK4 1. 01HaplotypeCaller <br> Script for GATK HaplotypeCaller 2. 02GenomicsDBImport\_GenotypeGVCFs <br> Script for GATK GenomicsDBImport and GenotypeGVCFs 3. 03filtering\_merge <br> Script for merging and filtering VCF records 4. 04consensus <br> Script for obtaining consensus sequences from filtered VCF files 3. **03repeats\_and\_gene\_annotation** <br> Repeat annotation and gene annotation of the reference sequence 1. 01repeat\_annotation <br> Repeat annotation with RepeatModeller and RepeatMasker 2. 02gene\_annotation <br> Gene annotation 1. 01mapping\_merge <br> Script for RNA-seq reads mapping and subsequent merging 2. 02braker2 <br> Gene annotation with braker2 pipeline 1. 01braker2\_RNAseq <br> Gene prediction using RNA-seq data with braker2 2. 02braker2\_protein <br> Gene prediction using protein sequences of related species with braker2 3. 03TSRBRA <br> Merging two gene prediction files by TSEBRA 4. 04filtering <br> Retaining only predicted genes whose expression are confirmed by RNA-seq data 5. 05reformatting\_by\_agat <br> Reformatting og gff file by agat 6. 06blast <br> Homology search against closely related protein sequences using blastp 4. **04population\_recombination\_rate\_estimation** <br> Estimation of genome wide population recombination rate with LDhelmet 1. 01calculate\_theta <br> Calculating watterson's theta 2. 02vcf2fasta <br> Converting phased VCF to fasta 3. 03find\_confs <br> 1st step of LDhelmet recombination rate estimation 4. 04table\_gen <br> 2nd step of LDhelmet recombination rate estimation 5. 05pade <br> 3rd step of LDhelmet recombination rate estimation 6. 06rjmcmc <br> 4th step of LDhelmet recombination rate estimation 7. 07post\_to\_text <br> Final step of LDhelmet recombination rate estimation 8. 08window <br> Summarizing the estimated population recombination rates in non-overlapping sliding windows 1. 01\_10kb <br> window size: 10kb 2. 02\_30kb <br> window size: 30kb 3. 03\_50kb <br> window size: 50kb 4. 04\_100kb <br> window size: 100kb 5. **05potentially\_deleterious\_SNPs** <br> Estimatiing potentially deleterious SNPs with PROVEAN 1. 01SNPEff <br> Extracting non-synonymous mutations by SNP annotation with SNPEff 2. 02provean <br> Assessing the potential deleteriousness of each non-synonymous mutations with PROVEAN 6. **06popultion\_genetic\_analyses** <br> Population genetic analyses using ddRAD-seq dataset 1. 01ADMIXTURE <br> Scripts for clustering analysis with ADMIXTURE 2. 02PCA <br> Scripts for principal component analysis with plink 3. 03diverisity\_indices <br> Scripts for calculating nucleotide diversity and pairwise \*F<sub>ST </sub>\*between populations 7. **07phylogenetic\_analysis** <br> Phylogenetic analysis with RAxML 1. 01ddRADseq\_RAxML <br> Scripts for phylogenetic analysis for ddRAD-seq dataset 2. 02WGS\_RAxML <br> Scripts for phylogenetic analysis for WGS nuclear SNPs dataset (with-outgroup dataset) 3. 03whole\_mitogenome\_RAxML <br> Scripts for phylogenetic analysis for whole mitogenome sequences obtained from WGS reads 8. **08hybrid\_detection** <br> Allele sharing pattern analysis to detect hybridization between populations / lineages 1. 01ddRADseq\_analysis <br> Analysis for ddRAD-seq dataset 1. 01Dsuite <br> *D*-statistics calculation and \*f-\*branch analysis with D-suite 2. 02HyDe <br> HyDe analysis 2. 02WGS\_analysis <br> Analysis for WGS dataset 1. 01Dsuite <br> *D*-statistics calculation and \*f-\*branch analysis with D-suite 2. 02HyDe <br> HyDe analysis 9. **09demographic\_estimation** <br> Demographic analyses 1. 01PSMC <br> PSMC analysis to infer population size history 1. 01coverage\_calculation <br> Calculating of coverage for each BAM file 2. 02PSMC <br> Scripts to run PSMC 3. 03PSMC\_bootstrap <br> Scripts to perform bootstrapping analysis of PSMC 2. 02demographic modeling <br> Demographic modeling analysis to test the ghost introgression origin of the ECS lineage 1. 01easySFS <br> Generating site frequency spectrum (SFS) from VCF with easySFS 2. 02fastsimcoal2 <br> Demographic modeling with fastsimcoal2 1. 01Examined\_models <br> All .est and .tpl files describing the demography and parameter rules of the examined models All models are classified into three categories based on the different elements incorporated into the model. The description of the bottom group of directories is omitted because only est and tpl are included in the directories of the corresponding models 1. 01without\_recent\_size\_change <br> Models without recent population size change (16 models) 2. 02with\_recent\_sice\_change <br> Models incorporating recent population size change (16 models) 3. 03recent\_gene\_flow\_without\_recent\_size\_change <br> Additional models to examine whether recent gene flow can affect the results or not (four models) 2. 02scripts\_to\_run\_fsc <br> Scripts to run fastsimcoal2 10. **10introgression\_landscape\_characterization** <br> Characterizing genomic landscape of introgression from the ghost lineage 1. 01validation\_by\_forward\_simulation <br> Evaluation of perfomance on estinating genomic landscape of introgression from the ghost lineage using forward simulation with SLiM 1. 01\_750kb\_simulation <br> Preliminary forward simulations of 750 kb chromosome to determine demographic parameters 1. 01slim <br> Scripts to run SLiM.Contains SLiM file (.slim) with adjusted demographic parameters used in the final 2. 02vcftools <br> Filtering of simulated VCFs 3. 03insertion\_to\_mimic\_ghost\_introgression <br> Script to mimic ghost introgression for simulated VCF, and to compare diversity indices with the observed data 2. 02\_10kb\_simulation <br> Forward simulations of 10 kb chromosome to evaluate performance of each statistic 1. 01slim <br> Scripts to run SLiM 2. 02vcftools <br> Filtering of simulated VCFs 3. 03insertion\_to\_mimic\_ghost\_introgression <br> Script to mimic ghost introgression for simulated VCF with different admixture rate 4. 04sliding\_window <br> Sliding window analysis on the simulated VCF 1. 01geno\_generate <br> Script to prepare input file (.geno) 2. 02sliding window <br> Scripts to run sliding window analysis 3. 03\_30kb\_simulation <br> Forward simulations of 30 kb chromosome to evaluate performance of each statistic 1. 01slim <br> Scripts to run SLiM 2. 02vcftools <br> Filtering of simulated VCFs 3. 03insertion\_to\_mimic\_ghost\_introgression <br> Script to mimic ghost introgression for simulated VCF with different admixture rate 4. 04sliding\_window <br> Sliding window analysis on the simulated VCF 1. 01geno\_generate <br> Script to prepare input file (.geno) 2. 02sliding window <br> Scripts to run sliding window analysis 4. 04\_50kb\_simulation <br> Forward simulations of 50 kb chromosome to evaluate performance of each statistic 1. 01slim <br> Scripts to run SLiM 2. 02vcftools <br> Filtering of simulated VCFs 3. 03insertion\_to\_mimic\_ghost\_introgression <br> Script to mimic ghost introgression for simulated VCF with different admixture rate 4. 04sliding\_window <br> Sliding window analysis on the simulated VCF 1. 01geno\_generate <br> Script to prepare input file (.geno) 2. 02sliding window <br> Scripts to run sliding window analysis 5. 05\_100kb\_simulation <br> Forward simulations of 100 kb chromosome to evaluate performance of each statistic 1. 01slim <br> Scripts to run SLiM 2. 02vcftools <br> Filtering of simulated VCFs 3. 03insertion\_to\_mimic\_ghost\_introgression <br> Script to mimic ghost introgression for simulated VCF with different admixture rate 4. 04sliding\_window <br> Sliding window analysis on the simulated VCF 1. 01geno\_generate <br> Script to prepare input file (.geno) 2. 02sliding window <br> Scripts to run sliding window analysis 2. 02sliding\_window <br> Sliding window analysis on the WGS SNPs dataset (with-outgroup dataset) 1. geno\_generate <br> Script to prepare input file (.geno) 2. 02sliding\_window\_10kb\_example <br> Example scripts to perform sliding window analysis and subsequent characterization of introgression landscape 1. 01sliding\_window <br> Scripts to run sliding window analysis 2. 02characterizing\_landscape\_by\_gamma <br> Scripts to characterize introgression landscape by sliding-window *γ* 1. 01summarizing\_all\_features <br> Summarizing the target features into sliding windows <br> 1\. 01repeat <br> Proportion of repetitive sequence <br> 2\. 02CDS <br> Proportion of coding sequence (CDS) <br> 3\. 03N\-mt <br> Proportion of CDS of the representative N-mt genes (OXPHOS and mitoribosomal genes) <br> 4\. 04OXPHOS <br> Proportion of CDS of the oxidative phosphorylation (OXPHOS) genes <br> 5\. 05mitoribo <br> Proportion of CDS of the mitochondrial ribosomal (mitoribosomal) genes <br> 6\. 06deleterious <br> Counts of potential deleterious alleles per lineage <br> 7\. 07CDS\_length <br> The length of CDS per window <br> 8\. 08summarizing\_all <br> Script to compile all features (including *γ* and population recombination rates) into bed files 2. 02random\_permutation <br> Random permutation to test whether the characteristics of each *γ* category differ from the genomic background <br> 1\. 01other\_than\_recomb <br> Tests for characteristics other than population recombination rate <br> ##1. 01permutation <br> ##Scripts to perform random permutations <br> 2\. 02recomb <br> Test for population recombination rate <br> ##1. 01permutation <br> ##Scripts to perform random permutations 3. 03circle\_permutation <br> Circular permutation to test whether the characteristics of each *γ* category differ from the genomic background <br> 1\. 01other\_than\_recomb <br> Tests for characteristics other than population recombination rate <br> ##1. 01permutation <br> ##Scripts to perform circular permutations <br> 2\. 02recomb <br> Test for population recombination rate <br> ##1. 01permutation <br> ##Scripts to perform circular permutations 4. 04controlled\_permutation <br> "Controlled permutation" to test whether the proportion of repetitive sequence in each *γ* category differs from the genomic background after controlling the proportion of CDS or the population recombination rate <br> 1\. 01control\_CDS <br> When controlling CDS <br> ##1. 01permutation <br> ##Scripts to perform controlled permutations <br> 2\. 01control\_averec <br> When controlling population recombination rate (averaged for three lineages) <br> ##1. 01permutation <br> ##Scripts to perform controlled permutations 3. genomics\_general <br> Modified script for two scripts in genomics\_general provided by Dr. Simon Marin. (Please see [https://github.com/simonhmartin/genomics_general](https://github.com/simonhmartin/genomics_general)) ### "genotyping\_data.tar.gz" Filtered genotyping dataset used in this study. 1. 01RADseq <br> Genotyping datasets for ddRAD-seq data (VCF files with SNPs / SNPs and invariant sites). Different merging and filtering were applied according to the purpose of the analysis (please see Supplementary Notes and Figure S1 in the paper). 1. RAD\_dataset1.vcf.gz <br> For population structure analysis (PCA, ADMIXTURE) 2. RAD\_dataset2.vcf.gz <br> For caluculation of diversity indices 3. RAD\_dataset3.vcf.gz <br> For phylogenetic analysis 4. RAD\_dataset4.vcf.gz <br> For allele sharing pattern analysis to detect hybridization (*D*-suite and HyDe) 5. RAD\_dataset5.vcf.gz <br> For demographic modelling 2. 02WGS <br> Genotyping datasets for WGS data. 1. 01nuclear\_SNP <br> Genotyping datasets for whole genome nuclear SNPs (VCF files with SNPs). The two datasets differ only in the presence or absence of outgroups. 1. WGS\_C.annularis\_only\_SNPs.vcf.gz <br> Whole genome nuclear SNPs only for 18 individuals of *Chaenogobius annularis* (without outgroup) 2. WGS\_with-outgroup\_SNPs.vcf.gz <br> Whole genome nuclear SNPs for 18 individuals of *C. annularis* and one individual of *C. gulosus* (with outgroup) 2. 02mitogenome\_sequence 1. whole\_mitogenome\_C.annularis.fasta <br> Whole mitogenome sequence multi-FASTA file for 18 individuals of *C. annularis* obtained from WGS data. ### "annotation\_data.tar.gz" 1. 01repeat\_annotation 1. agohaze\_sspace\_x1.fa.masked.gz <br> Fasta file of reference genome (accession number: GCA\_015082035.1) soft-masked for the repetitive sequence 2. agohaze\_sspace\_x1.fa.out <br> Catalog of repetitive sequences in the reference genome output by RepeatMasker 2. 02gene\_annotation 1. agat\_braker\_FULL.gff3 <br> Gene annotation files after reformatting by agat. It consists only of predicted genes whose expression were confirmed by RNA-seq data ("conservative annotation" in the paper) ### "demographic\_modeling.tar.gz" 1. 01without\_recent\_size\_change 1. sorted\_AIC\_dist\_wo\_recent.csv <br> The AIC distributions for each model without recent population size change, with the model name in the first column and the AIC in the second column. Input data for Figure 4a in the paper. 2. best\_parameter\_set\_in\_best\_models 1. model7a 1. model7a.bestlhoods <br> Best parameters in model7a. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/01without\_recent\_size\_change/model7a" in "scripts.tar.gz" and Figure S13 in the paper. 2. model7b 1. model7b.bestlhoods <br> Best parameters in model7b. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/01without\_recent\_size\_change/model7b" in "scripts.tar.gz" and Figure S13 in the paper. 3. model7c 1. model7c.bestlhoods <br> Best parameters in model7c. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/01without\_recent\_size\_change/model7c" in "scripts.tar.gz" and Figure S13 in the paper. 2. 02with\_recent\_size\_change 1. sorted\_AIC\_dist\_with\_recent.csv <br> The AIC distributions for each model with recent population size change, with the model name in the first column and the AIC in the second column.Input data for Figure 4d in the paper. 1. best\_parameter\_set\_in\_best\_models 1. model7a 1. model7a.bestlhoods <br> Best parameters in model7a. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/02with\_recent\_sice\_change/model7a" in "scripts.tar.gz" and Figure S14 in the paper. 2. model7b 1. model7b.bestlhoods <br> Best parameters in model7b. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/02with\_recent\_sice\_change/model7b" in "scripts.tar.gz" and Figure S14 in the paper. 3. model7c 1. model7c.bestlhoods <br> Best parameters in model7c. Each parameter is listed in "09demographic\_estimation/02demographic\_modeling/02fastsimcoal2/01 Examined\_models/02with\_recent\_sice\_change/model7c" in "scripts.tar.gz" and Figure S14 in the paper. 3. 03recent\_gene\_flow\_without\_recent\_size\_change 1. sorted\_AIC\_rec-ongmig\_without\_recent\_population\_size.csv <br> The AIC distributions for additional models incorporating recent gene flow (without recent population size change). The model name in the first column and the AIC in the second column. Input data for Figure S15 in the paper. 4. 04observed\_MSFS 1. MSFS.obs <br> Multi-dimensions site frequency spectrum (MSFS) used as an input file for the demographic modeling. It was generated from RAD\_dataset5.vcf.gz using easySFS ### "slidingwindow\_results.tar.gz" 1. README\_description\_of\_record\_XXkb.txt <br> Text file describing each column of record\_XXkb.bed.gz 2. record\_10kb.bed.gz 3. record\_30kb.bed.gz 4. record\_50kb.bed.gz 5. record\_100kb.bed.gz ## **Sharing/Access information** Sequencial data used in this study can be available from DDBJ (accession numbers: DRR174909, DRR175781–DRR175796, DRR175830–DRR175860, DRR175876–DRR175955, DRR489922–DRR490073 for ddRAD-seq, DRR489903–DRR489921 for whole genome resequencing, and DRR490074–DRR490084 for RNA-seq).