Genomic landscape of introgression from the ghost lineage in a gobiid fish uncovers the generality of forces shaping hybrid genomes

Kato, Shuya1

Published Nov 06, 2023; Updated Dec 18, 2023 on Dryad. https://doi.org/10.5061/dryad.7wm37pw09

Data files

Nov 06, 2023 version files 4.12 GB

Dec 18, 2023 version files 4.12 GB

Abstract

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.

README: 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

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.

Detailed description of the file structure

"scripts.tar.gz"

01ddRADseq_genotyping Scripts used for genotyping double digest restriction-site associated DNA (ddRAD-seq) data
1. 01trimmomatic Read filtering with Trimmomatic
2. 02mapping Mapping with bwa-mem
3. 03sambamba Filtering of BAM files by sambamba (retaining only uniquely mapped reads)
4. 04SNPcall_by_mpileup Genotyping with samtools mpileup / bcftools call
 1. 01pop1_example Example of a script for joint calling when targeting one group
5. 05merge_filtering Merging and filtering of VCF files to create the final genotype datasets
 1. 01dataset1_2 Merging and filtering scripts for datasets 1 and 2
 1. 01filtering_dataset1 Filtering scripts for dataet 1
 2. 02filtering_dataset2 Filtering scripts for dataet 2
 2. 02dataset3 Merging and filtering scripts for dataset 3
 1. 01filtering Filtering scripts for dataet 3
 3. 03dataset4 Merging and filtering scripts for dataset 4
 1. 01filtering Filtering scripts for dataet 4
 4. 04dataset5 Merging and filtering scripts for dataset 5
 1. 01filtering Filtering scripts for dataet 5
02WGS_genotyping Scripts used for genotyping whole genome resequencing (WGS) data
1. 01nuclear_SNP Genotyping whole nuclear SNPs
 1. 01mapping Mapping with bwa-mem
 2. 02sambamba Filtering of BAM files by sambamba (retaining only uniquely mapped reads)
 3. 03MarkDuplicate Removing PCR duplicates in BAMfiles by MarkDuplicate (GATK4)
 4. 04SNPcall_by_GATK Genotyping with GATK4
 1. 01SJ Genotyping for the Sea of Japan (SJ) lineage
 1. 01HaplotypeCaller Script for GATK HaplotypeCaller
 2. 02GenomicsDBImport_GenotypeGVCFs Script for GATK GenomicsDBImport and GenotypeGVCFs
 2. 02ECS Genotyping for the East China Sea (ECS) lineage
 1. 01HaplotypeCaller Script for GATK HaplotypeCaller
 2. 02GenomicsDBImport_GenotypeGVCFs Script for GATK GenomicsDBImport and GenotypeGVCFs
 3. 03PO Genotyping for the Pacific Ocean (PO) lineage
 1. 01HaplotypeCaller Script for GATK HaplotypeCaller
 2. 02GenomicsDBImport_GenotypeGVCFs Script for GATK GenomicsDBImport and GenotypeGVCFs
 4. 04outgroup Genotyping for the outgroup individual
 1. 01HaplotypeCaller Script for GATK HaplotypeCaller
 2. 02GenomicsDBImport_GenotypeGVCFs Script for GATK GenomicsDBImport and GenotypeGVCFs
 5. 05merge_filtering Merging and filtering of VCF files to create the final genotype datasets
 1. 01Cannularis_only_dataset Merging and filtering scripts for C. annularis only dataset
 2. 02with_outgroup_dataset Merging and filtering scripts for with-outgroup dataset
 5. 05indelrealign Local realignment around indels for BAM files with GATK3
 6. 06phasing Read-aware-phasing with shape-it2
2. 02mitogenome_sequence Genotyping whole mitogenome sequence
 1. 01mapping Mapping with bwa-mem
 2. 02sambamba Filtering of BAM files by sambamba (retaining only uniquely mapped reads)
 3. 03MarkDuplicates Removing PCR duplicates in BAMfiles by MarkDuplicate (GATK4)
 4. 04SNPcall_by_GATK Genotyping with GATK4
 1. 01HaplotypeCaller Script for GATK HaplotypeCaller
 2. 02GenomicsDBImport_GenotypeGVCFs Script for GATK GenomicsDBImport and GenotypeGVCFs
 3. 03filtering_merge Script for merging and filtering VCF records
 4. 04consensus Script for obtaining consensus sequences from filtered VCF files
03repeats_and_gene_annotation Repeat annotation and gene annotation of the reference sequence
1. 01repeat_annotation Repeat annotation with RepeatModeller and RepeatMasker
2. 02gene_annotation Gene annotation
 1. 01mapping_merge Script for RNA-seq reads mapping and subsequent merging
 2. 02braker2 Gene annotation with braker2 pipeline
 1. 01braker2_RNAseq Gene prediction using RNA-seq data with braker2
 2. 02braker2_protein Gene prediction using protein sequences of related species with braker2
 3. 03TSRBRA Merging two gene prediction files by TSEBRA
 4. 04filtering Retaining only predicted genes whose expression are confirmed by RNA-seq data
 5. 05reformatting_by_agat Reformatting og gff file by agat
 6. 06blast Homology search against closely related protein sequences using blastp
04population_recombination_rate_estimation Estimation of genome wide population recombination rate with LDhelmet
1. 01calculate_theta Calculating watterson's theta
2. 02vcf2fasta Converting phased VCF to fasta
3. 03find_confs 1st step of LDhelmet recombination rate estimation
4. 04table_gen 2nd step of LDhelmet recombination rate estimation
5. 05pade 3rd step of LDhelmet recombination rate estimation
6. 06rjmcmc 4th step of LDhelmet recombination rate estimation
7. 07post_to_text Final step of LDhelmet recombination rate estimation
8. 08window Summarizing the estimated population recombination rates in non-overlapping sliding windows
 1. 01_10kb window size: 10kb
 2. 02_30kb window size: 30kb
 3. 03_50kb window size: 50kb
 4. 04_100kb window size: 100kb
05potentially_deleterious_SNPs Estimatiing potentially deleterious SNPs with PROVEAN
1. 01SNPEff Extracting non-synonymous mutations by SNP annotation with SNPEff
2. 02provean Assessing the potential deleteriousness of each non-synonymous mutations with PROVEAN
06popultion_genetic_analyses Population genetic analyses using ddRAD-seq dataset
1. 01ADMIXTURE Scripts for clustering analysis with ADMIXTURE
2. 02PCA Scripts for principal component analysis with plink
3. 03diverisity_indices Scripts for calculating nucleotide diversity and pairwise *F_ST*between populations
07phylogenetic_analysis Phylogenetic analysis with RAxML
1. 01ddRADseq_RAxML Scripts for phylogenetic analysis for ddRAD-seq dataset
2. 02WGS_RAxML Scripts for phylogenetic analysis for WGS nuclear SNPs dataset (with-outgroup dataset)
3. 03whole_mitogenome_RAxML Scripts for phylogenetic analysis for whole mitogenome sequences obtained from WGS reads
08hybrid_detection Allele sharing pattern analysis to detect hybridization between populations / lineages
1. 01ddRADseq_analysis Analysis for ddRAD-seq dataset
 1. 01Dsuite D-statistics calculation and *f-*branch analysis with D-suite
 2. 02HyDe HyDe analysis
2. 02WGS_analysis Analysis for WGS dataset
 1. 01Dsuite D-statistics calculation and *f-*branch analysis with D-suite
 2. 02HyDe HyDe analysis
09demographic_estimation Demographic analyses
1. 01PSMC PSMC analysis to infer population size history
 1. 01coverage_calculation Calculating of coverage for each BAM file
 2. 02PSMC Scripts to run PSMC
 3. 03PSMC_bootstrap Scripts to perform bootstrapping analysis of PSMC
2. 02demographic modeling Demographic modeling analysis to test the ghost introgression origin of the ECS lineage
 1. 01easySFS Generating site frequency spectrum (SFS) from VCF with easySFS
 2. 02fastsimcoal2 Demographic modeling with fastsimcoal2
 1. 01Examined_models 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 Models without recent population size change (16 models)
 2. 02with_recent_sice_change Models incorporating recent population size change (16 models)
 3. 03recent_gene_flow_without_recent_size_change Additional models to examine whether recent gene flow can affect the results or not (four models)
 2. 02scripts_to_run_fsc Scripts to run fastsimcoal2
10introgression_landscape_characterization Characterizing genomic landscape of introgression from the ghost lineage
1. 01validation_by_forward_simulation Evaluation of perfomance on estinating genomic landscape of introgression from the ghost lineage using forward simulation with SLiM
 1. 01_750kb_simulation Preliminary forward simulations of 750 kb chromosome to determine demographic parameters
 1. 01slim Scripts to run SLiM.Contains SLiM file (.slim) with adjusted demographic parameters used in the final
 2. 02vcftools Filtering of simulated VCFs
 3. 03insertion_to_mimic_ghost_introgression Script to mimic ghost introgression for simulated VCF, and to compare diversity indices with the observed data
 2. 02_10kb_simulation Forward simulations of 10 kb chromosome to evaluate performance of each statistic
 1. 01slim Scripts to run SLiM
 2. 02vcftools Filtering of simulated VCFs
 3. 03insertion_to_mimic_ghost_introgression Script to mimic ghost introgression for simulated VCF with different admixture rate
 4. 04sliding_window Sliding window analysis on the simulated VCF
 1. 01geno_generate Script to prepare input file (.geno)
 2. 02sliding window Scripts to run sliding window analysis
 3. 03_30kb_simulation Forward simulations of 30 kb chromosome to evaluate performance of each statistic
 1. 01slim Scripts to run SLiM
 2. 02vcftools Filtering of simulated VCFs
 3. 03insertion_to_mimic_ghost_introgression Script to mimic ghost introgression for simulated VCF with different admixture rate
 4. 04sliding_window Sliding window analysis on the simulated VCF
 1. 01geno_generate Script to prepare input file (.geno)
 2. 02sliding window Scripts to run sliding window analysis
 4. 04_50kb_simulation Forward simulations of 50 kb chromosome to evaluate performance of each statistic
 1. 01slim Scripts to run SLiM
 2. 02vcftools Filtering of simulated VCFs
 3. 03insertion_to_mimic_ghost_introgression Script to mimic ghost introgression for simulated VCF with different admixture rate
 4. 04sliding_window Sliding window analysis on the simulated VCF
 1. 01geno_generate Script to prepare input file (.geno)
 2. 02sliding window Scripts to run sliding window analysis
 5. 05_100kb_simulation Forward simulations of 100 kb chromosome to evaluate performance of each statistic
 1. 01slim Scripts to run SLiM
 2. 02vcftools Filtering of simulated VCFs
 3. 03insertion_to_mimic_ghost_introgression Script to mimic ghost introgression for simulated VCF with different admixture rate
 4. 04sliding_window Sliding window analysis on the simulated VCF
 1. 01geno_generate Script to prepare input file (.geno)
 2. 02sliding window Scripts to run sliding window analysis
2. 02sliding_window Sliding window analysis on the WGS SNPs dataset (with-outgroup dataset)
 1. geno_generate Script to prepare input file (.geno)
 2. 02sliding_window_10kb_example Example scripts to perform sliding window analysis and subsequent characterization of introgression landscape
 1. 01sliding_window Scripts to run sliding window analysis
 2. 02characterizing_landscape_by_gamma Scripts to characterize introgression landscape by sliding-window γ
 1. 01summarizing_all_features Summarizing the target features into sliding windows 1. 01repeat Proportion of repetitive sequence 2. 02CDS Proportion of coding sequence (CDS) 3. 03N-mt Proportion of CDS of the representative N-mt genes (OXPHOS and mitoribosomal genes) 4. 04OXPHOS Proportion of CDS of the oxidative phosphorylation (OXPHOS) genes 5. 05mitoribo Proportion of CDS of the mitochondrial ribosomal (mitoribosomal) genes 6. 06deleterious Counts of potential deleterious alleles per lineage 7. 07CDS_length The length of CDS per window 8. 08summarizing_all Script to compile all features (including γ and population recombination rates) into bed files
 2. 02random_permutation Random permutation to test whether the characteristics of each γ category differ from the genomic background 1. 01other_than_recomb Tests for characteristics other than population recombination rate ##1. 01permutation ##Scripts to perform random permutations 2. 02recomb Test for population recombination rate ##1. 01permutation ##Scripts to perform random permutations
 3. 03circle_permutation Circular permutation to test whether the characteristics of each γ category differ from the genomic background 1. 01other_than_recomb Tests for characteristics other than population recombination rate ##1. 01permutation ##Scripts to perform circular permutations 2. 02recomb Test for population recombination rate ##1. 01permutation ##Scripts to perform circular permutations
 4. 04controlled_permutation "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 1. 01control_CDS When controlling CDS ##1. 01permutation ##Scripts to perform controlled permutations 2. 01control_averec When controlling population recombination rate (averaged for three lineages) ##1. 01permutation ##Scripts to perform controlled permutations
 3. genomics_general Modified script for two scripts in genomics_general provided by Dr. Simon Marin. (Please see https://github.com/simonhmartin/genomics_general)

"genotyping_data.tar.gz"

Filtered genotyping dataset used in this study.

01RADseq 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 For population structure analysis (PCA, ADMIXTURE)
2. RAD_dataset2.vcf.gz For caluculation of diversity indices
3. RAD_dataset3.vcf.gz For phylogenetic analysis
4. RAD_dataset4.vcf.gz For allele sharing pattern analysis to detect hybridization (D-suite and HyDe)
5. RAD_dataset5.vcf.gz For demographic modelling
02WGS Genotyping datasets for WGS data.
1. 01nuclear_SNP 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 Whole genome nuclear SNPs only for 18 individuals of Chaenogobius annularis (without outgroup)
 2. WGS_with-outgroup_SNPs.vcf.gz 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 Whole mitogenome sequence multi-FASTA file for 18 individuals of C. annularis obtained from WGS data.

"annotation_data.tar.gz"

01repeat_annotation
1. agohaze_sspace_x1.fa.masked.gz Fasta file of reference genome (accession number: GCA_015082035.1) soft-masked for the repetitive sequence
2. agohaze_sspace_x1.fa.out Catalog of repetitive sequences in the reference genome output by RepeatMasker
02gene_annotation
1. agat_braker_FULL.gff3 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"

01without_recent_size_change
1. sorted_AIC_dist_wo_recent.csv 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 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 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 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.
02with_recent_size_change
1. sorted_AIC_dist_with_recent.csv 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 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 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 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.
03recent_gene_flow_without_recent_size_change
1. sorted_AIC_rec-ongmig_without_recent_population_size.csv 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.
04observed_MSFS
1. MSFS.obs 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"

README_description_of_record_XXkb.txt Text file describing each column of record_XXkb.bed.gz
record_10kb.bed.gz
record_30kb.bed.gz
record_50kb.bed.gz
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).