Homoploid hybrid speciation in a marine pelagic fish

Homoploid hybrid speciation (HHS) is an enigmatic evolutionary process where new species arise through hybridization of divergent lineages without changes in chromosome number. Although increasingly documented in various taxa and ecosystems, convincing cases of HHS in marine fishes have been lacking. This study presents a possible case of HHS in a pelagic marine fish based on comprehensive genomic, morphological, and ecological analyses. Population genomics, species tree estimation, and tests of introgression and admixture identified three sympatric clusters in Megalaspis cordyla in the western Pacific and the admixed nature of one cluster between the others. Moreover, model-based demographic inference favored a hybrid speciation scenario over introgression for the origin of the admixed cluster. While contemporary gene flow suggested partial reproductive isolation, examination of occurrence data and ecologically relevant morphological characters suggested ecological differences between the clusters, potentially contributing to the reproductive isolation and niche partitioning in sympatry. The clusters are also morphologically distinguishable and thus can be taxonomically recognized as separate species. The hybrid cluster is restricted to the coasts of Taiwan and Japan, where all three clusters coexist. The parental clusters are additionally found in lower latitudes, where they display non-overlapping distributions. Given the geographical distributions, estimated times of the species formation, and patterns of historical demographic changes, we propose that the Pleistocene glacial cycles were the primary driver of HHS in this system. We also develop an ecogeographic model of HHS in marine coastal ecosystems, including a novel hypothesis to explain the initial stages of HHS.

Methods
We analyzed a total of 160 specimens of M. cordyla collected from various localities in the western Pacific. 

We obtained partial sequences of the mitochondrial cytochrome b gene (CYTB: 873 bp) from all specimens. Additionally, we obtained partial sequences of the following seven loci from a subset of the specimens: cytochrome c oxidase subunits 1 (COI: 612 bp) and 2 (COII: 657 bp), NADH dehydrogenase subunits 1 (ND1: 828 bp), 2 (ND2: 531 bp), 4 (ND4: 828 bp), and 5 (ND5: 801 bp), and the control region (CR: 760 bp). These sequences were used for divergent time estimation, haplotype network construction, and Neighbor-Joining tree estimation.

We genotyped 115 specimens using MIG-seq, a reduced representation sequencing method targeting SNPs in the inter-simple sequence repeat region. We processed the raw MIG-seq reads by the single-end mode of fastp v0.20.1, merging the processed reads 1 and 2 into a single FASTQ file, and mapped them to the Caranx ignobilis reference genome (GenBank JAFHLA000000000.1) using BWA-mem v0.7.17 with default settings. Unmapped or low-quality reads were removed using Samtools v1.12. The resulting BAM files were processed using Stacks v2.5.4, bcftools v1.8, VCFtools v0.1.16, and R package ‘SNPfiltR’ v1.0.0 to perpare suitable datasets for downstream analyses.

We examined 37 linear measurements and 17 countable characters in 37 specimens. 

The ecological data is based on the catch data of M. cordyla in Taiwan during 2000–2001, originally reported in Su et al. (2020:Scientific Reports 10: 16829).

More details of data collection, processing, and analyses are described in the Supplementary Material (Supplement_text.docx).