Extensive Interspecific Gene Flow Shaped Complex Evolutionary History and Underestimated Species Diversity in Rapidly Radiated Dolphins

Recently diverged taxa are often characterized by high rates of hybridization, which can complicate phylogenetic reconstruction. For this reason, the phylogenetic relationships and evolutionary history of dolphins are still not very well resolved; the question of whether the genera Tursiops and Stenella are monophyletic is especially controversial. Here, we performed re-sequencing of six dolphin genomes and combined them with eight previously published dolphin SRA datasets and six whole-genome datasets to investigate the phylogenetic relationships of dolphins and test the monophyly hypothesis of Tursiops and Stenella. Phylogenetic reconstruction with the maximum likelihood and Bayesian methods of concatenated loci, as well as with coalescence analyses of sliding window trees, produced a concordant and well-supported tree. Our studies support the non-monophyletic status of Tursiops and Stenella because the species referred these genera do not form exclusive monophyletic clades. This suggests that the current taxonomy of both genera might not reflect their evolutionary history and may underestimate their diversity. A four-taxon D-statistic (ABBA-BABA) test, five-taxon DFOIL test, and tree-based PhyloNet analyses all showed extensive gene flow across dolphin species, which could explain the instability in resolving phylogenetic relationship of oceanic dolphins with different and limited markers. This study could be a good case to demonstrate how genomic data can reveal complex speciation and phylogeny in rapidly radiating animal groups.


Introduction
Cetaceans, including whales, dolphins, and porpoises, are a group of secondarily adapted marine mammals with a history of transitioning from terrestrial (land) to fully aquatic habitats, representing one of the most fascinating evolutionary events of the planet (Thewissen et al. 2009;Uhen 2010). Approximately 56-53 million years ago, cetacean ancestors transformed their habitat from a completely terrestrial to a fully aquatic environment and subsequently radiated around the world with a series of outstanding morphological changes (e.g., degeneration of limbs, lengthening of the skull, thickening of the blubber, loss of hair, etc.) (Reidenberg 2007;Thewissen et al. 2007). After they transitioned to the sea, cetaceans further diverged into two highly distinct groups: Mysticeti (baleen whales, characterized by baleen plates and two blowholes) and Odontoceti (toothed whales, dolphins, and porpoises, which have a single blowhole) (Gatesy et al. 2013).
Establishing correct phylogenetic relationships is an elementary step in unraveling the evolutionary processes responsible for cetacean diversity. The secondarily aquatic nature of whales has inspired numerous phylogenetic studies and sparked debates about their origin from terrestrial mammals (Gatesy and O'Leary 2001;Thewissen et al. 2001;Árnason et al. 2004;O'Leary and Gatesy 2008;Uhen 2010). Several problems concerning the systematics of cetaceans have been well resolved, such as the monophyly of cetaceans, Mysticeti, and Odontoceti (Gatesy et al. 1999;Nikaido et al. 1999;Geisler and Sanders 2003), the sister relationship between cetaceans and Hippopotamidae (Lum et al. 2000;McGowen et al. 2009;Thewissen et al. Weijian Guo and Di Sun Both authors contributed equally. 2009; Zhou et al. 2011), and the non-monophyletic status of river dolphins (Geisler et al. 2011;Hassanin et al. 2012;McGowen et al. 2020). Despite widespread popular and scientific interest in some major cetacean lineages, many aspects of their phylogeny and evolutionary history remain unresolved. For example, the phylogenetic relationships within Balaenopteroidea and Ziphiidae have been studied using extensive DNA sequences, but several studies have produced contradictory results (Nikaido et al. 2006;Dalebout et al. 2008;McGowen et al. 2009;Steeman et al. 2009;Árnason et al. 2018). Beyond rorquals and beaked whales, the relationships within Delphinidae remains problematic.
Delphinidae is the largest cetacean family (~37 species), and resolving phylogenetic relationships among its members has been disheartening despite numerous recent attempts because of the rapidity of speciation, which facilitated incomplete lineage sorting (ILS) and hybridization between emerging species (Kingston et al. 2009;Amaral et al. 2012;Perrin et al. 2013). A particular conundrum surrounds the phylogenetic relationships among bottlenose-like dolphins (e.g., subfamily Delphininae), which also have experienced recent, rapid radiations (Steeman et al. 2009;Slater et al. 2010). Although some morphological (Shirakihara et al. 2003;Guidarelli et al. 2014Guidarelli et al. , 2018Jedensjö et al. 2017) and molecular (McGowen et al. 2009(McGowen et al. , 2020Chen et al. 2011;Moura et al. 2013Moura et al. , 2020 analyses have been conducted, the evolutionary relationships among Delphininae generally have remained unclear. In particular, the monophyly of the genera Tursiops and Stenella is still under debate. As early as 20 years ago, most taxonomic studies of dolphin evolutionary relationships relied on morphological characters, particularly those of the skull. For example, the studies of skull, teeth, and vertebrae showed that species of Tursiops and Stenella can be distinguished from those of other genera (Archer and Perrin 1999;Shirakihara et al. 2003). Ever since the comprehensive molecular analysis of Delphinidae by Leduc et al. (1999), which included 33 species and was based on full mitochondrial cytochrome b sequences, an increasing number of studies using genetic data including multiple nuclear loci (Kingston et al. 2009), short interspersed elements (SINE) (Chen et al. 2011), and even whole mitochondrial genomes (Hassanin et al. 2012;Moura et al. 2013) have not supported the monophyly of Tursiops or Stenella (i.e., the species of these two groups do not cluster with their congeners). These results have partly contrasted with other studies integrating nuclear and mitochondrial data (McGowen et al. 2009;McGowen 2011), multiple protein-coding genes (McGowen et al. 2020, and RADseq data (Moura et al. 2020), which have recovered Tursiops as a monophyletic group (Fig. 1).
Two evolutionary processes, hybridization and ILS, may trigger the uncertain phylogenetic relationships observed in some groups of mammals, such as wild pigs, bears, and baboons (Frantz et al. 2013;Kumar et al. 2017;Rogers et al. 2019). Due to the invariable number of chromosomes (2n = 44) and the common karyotic arrangement in most cetaceans, a recent review found that hybridization is a common phenomenon among oceanic cetaceans (Crossman et al. 2016). Incongruence of gene trees across the genome results from hybridization and ILS that can be detected by analyzing whole-genome sequences (Degnan and Rosenberg 2009). In this context, it is becoming increasingly important to use genome-scale data sets to resolve phylogenetic relationships among cetaceans.
The genomic era allows us to sample massive genomic data, which undoubtedly helps us to address the phylogenetic issues that cannot be well resolved using limited genetic markers (Delsuc et al. 2005). Moreover, genomic data could also provide a special opportunity to detect the footprints of hybridization and allow exhaustive analyses of how gene flow affects phylogenetic relationships (Payseur and Rieseberg 2016). Recently, Árnason et al. (2018) analyzed the evolutionary history of baleen whales using genomic data and found some gene flow across rorquals. Here, we present new genomic data from six oceanic dolphin species. Using some analytical methods specially developed to handle complex genomic data, the phylogenetic relationships of oceanic dolphins were investigated, providing new insights into whether the genera Tursiops and Stenella are monophyletic or non-monophyletic as well as signals of recent and ancestral introgression among dolphins.

DNA Isolation, Sequencing, and Pseudo-genome Generation
The verified specimens of six species (nine individuals) sequenced in the present study were preserved at Nanjing Normal University. All samples were collected from individuals that were incidentally killed by unknown reasons in the wild (and thus no ethical approval was necessary). All tissue samples were frozen at -20 °C before this study. DNA was isolated using a standard phenol-chloroform method (Sambrook et al. 1989). Sequencing libraries were prepared with insert sizes of 550 bp and sequenced using Illumina Nextseq 500 technology. The whole genomes were sequenced at 8 × to 25.3 × depth coverage (Online Resource 1). Quality control was performed using FastQC (www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/), and the retrieved raw paired reads were trimmed using Trimmomatic version 0.33 (Bolger et al. 2014) with parameters "SLIDINGWINDOW:4:15" to cut the read once the average quality within the window fell below a threshold and "LEADING:3" as well as "TRAILING:3", which removed low quality bases from the beginning and the end, respectively. As the common bottlenose dolphin is phylogenetically positioned in the genus Tursiops, and a possible mapping bias against its genome is likely to affect phylogenetic and gene flow analyses, the killer whale genome (Foote et al. 2015) with soft mask was used as a reference in order to avoid mapping biases that can affect phylogenetic analyses. Scaffolds in the killer whale genome shorter than 100 kb were excluded. Paired-end reads were mapped to the killer whale genome with BWA version 0.7.12-r1039 (Li and Durbin 2009) using the "mem" and "-M" option (The "mem" option stands for Maximal Exact Match and the "-M" option means mark shorter split hits as secondary), and duplicates were marked with Picard (https:// github. com/ broad insti tute/ picard). From the mapped reads, single-nucleotide variants (SNVs) and short insertion or deletions (Indels) were called by SAMtools ) with the default settings. Consensus sequences were created from VCF files using BCFtools (https:// samto ols. github. io/ bcfto ols/), choosing the allele with more reads (> = 5) in case of heterozygous sites. Indels were removed, and ambiguously called sites were masked as "N". At this point, the pseudo-genome of each individual was created.

Pseudo-genome Alignment
The generation of multiple sequence alignments is a crucial step for phylogenetic analyses. Our genomic sequence alignment (GSA) of all pseudo-genomes spanned 21 individuals from 12 cetacean species and an outgroup species  Moura et al. (2020). Abbreviations: nu (nuclear), mt (mitochondrial), SINE (short interspersed elements) (H. amphibius) as well as the genome of killer whale. First, the pairwise alignments were conducted by LAST version 963 package (Kiełbasa et al. 2011), with the killer whale genome as the reference genome. Firstly, the killer whale genome "database" was built using "lastdb" with the parameter "-cR11", which excluded lowercase letters from initial matches but kept the lowercase letters in the sequences. Secondly, each pseudo-genome was aligned to the reference genome using the "lastal" command with the parameter "-E0.05, -m100". Subsequently, the "maf-swap" command was used to change the order of the sequences in the MAFformat alignments to obtain the best pairwise aligned blocks. Lastly, the pairwise alignments were merged into multiple genome alignments by MULTIZ version 11.2 (Blanchette et al. 2004) using the killer whale genome as a reference. Only 1:1 aligning to the killer whale sequence was allowed, choosing the hits with the highest identity, to reduce the probability of finding paralogs. After the multiple-genome alignments were created, the repetitive regions (lowercase in the killer whale genome) were removed using in-house Python scripts.

Phylogenetic Tree Construction
In total, 2,562,378 syntenic blocks were concatenated using in-house Python scripts, and a FASTA-formatted alignment file containing ~1,014 Mb sequences was then generated in our prior analysis. First, non-overlapping genome fragments of 10 kb were extracted using custom Python scripts, creating 101,408 windows from the GSA. Second, for each window, 200 rapid bootstrap analyses were estimated by using RaxML version 8.2.12 (Stamatakis 2014) with the default parameters under the GTR + GAMMA model. In order to ensure that the segments were appropriate for phylogenomic analysis, all the window trees were checked using a custom script for abnormal clusters in which individuals from the same species showed non-monophyletic status, indicating a certain window provided wrong information. After the abnormal trees were filtered, 14,355 windows remained. The median number of parsimony informative sites across 23 individuals for these 10 kb windows was 366, which was relatively informative for further phylogeny reconstruction, and then a frequency distribution graph (Online Resource 2) was plotted by ggplot2 (Wickham 2016). After individual window-based trees were reconstructed, the average bootstrap support values of these trees were summarized using a custom Python script. Again, these 10 k windows were concatenated using in-house Python scripts to form a new matrix (hereafter referred to as NewGSA). Phylogenetic relationships were inferred using both maximum likelihood (ML) and Bayesian methods as implemented in RaxML and ExaBayes version 1.4 (Aberer et al. 2014) based on New-GSA with an unpartitioned scheme due to thousands of windows in NewGSA. Moreover, in our pilot analyses (with 100 bootstrap replicates), we found that partitioning schemes had almost no effect on the final phylogenetic results. For ML analyses, 500 bootstrap replicates were estimated by RaxML under the GTR + GAMMA model. For the Bayesian analyses, default parameters and the GTR + GAMMA model of evolution was used to run two independent unpartitioned analyses. Each run was conducted for 10 6 generations, and Metropolis-Coupling (numCoupledChains = 2) was used to speed up convergence. The initial 10% of runs were discarded as burn-in. Results were checked for convergence and effective sampling sizes of all parameters using the tools included with the ExaBayes package. Supporting values are expressed as posterior probabilities.
Standard concatenation approaches did not model discordance among gene trees beyond the differences in sequence evolution history (Gatesy and Baker 2005). Recent studies have suggested that concatenation methods could yield results that differ from those obtained by coalescent-based methods, especially in the presence of ILS (Degnan and Rosenberg 2006;Liu et al. 2015;Roch and Steel 2015). This possible limitation can theoretically be overcome with the gene-tree-based coalescent approach using software such as ASTRAL version 5.6.2 (Zhang et al. 2018). Here, the remained windows were used as the input to ASTRAL. Briefly, for each window, 200 rapid bootstrap analyses were estimated using RaxML with the default parameters under the GTR + GAMMA model. Furthermore, three tree subsets were generated, selecting windows with average bootstrapping values over 70% (8,941 windows), 75% (4,222 windows), and 80% (1,085 windows). ASTRAL was applied to reconstruct the species tree based on all window-based trees and the three tree subsets using the default parameters. For the species tree based on all windowbased trees, 200 multi-locus bootstrapping replicates were implemented using ASTRAL-III (option: -r 200) (Seo 2008). In addition, another coalescent-based approach was employed to estimate a species tree. SVDquartets Kubatko 2014, 2015), as implemented in PAUP* version 4.0a (Swofford 2003) was performed on the same 200 bootstrap replicates that we generated window-based tree based on ML analysis.

Divergence Time Calibration
Divergence dating analyses were conducted using the software MCMCTree version 4.9 h, part of the PAML package (Yang 2007), with the approximate likelihood method (dos Reis et al. 2012). Phylogenetic analyses of different genomic segments that were shaped by different evolutionary histories will yield trees with unique clade clusters and branch lengths, which can produce a range of time estimates in divergence dating analyses (Fontaine et al. 2015;Springer et al. 2019). To overcome difficulties arising from computational efficiency and complicated evolutionary histories, we only used 10 k window sequences that had the same topology as the main topology. Although these divergence times are recent-biased, the amount of error is expected to be relatively small considering the deep time scale in this analysis (Frantz et al. 2013). Meanwhile, we only used specieslevel datasets (retaining the highest individual coverage of particular species) to ensure a high level of computational efficiency. To do this, a species-level dataset was extracted by SeqKit version 0.12.0 (Shen et al. 2016) from the New-GSA. For our main analyses, correlated rate (clock = 3) was used, and the parameters of the prior birth-death process were fixed at birth rate λ = 1, death rate μ = 1, and sampling fraction ρ = 0.1, which led to a uniform distribution of node ages (Yang and Rannala 2006). The posterior distribution of parameters was estimated using Markov chain Monte Carlo (MCMC) sampling. MCMC runs were conducted twice for 10 8 iterations with a sampling frequency of 1000, with the first 10 4 steps of each run discarded as burn-in. The results were examined in Tracer version 1.7 (Rambaut et al. 2018) to evaluate whether parameters, node ages, and likelihood values had converged. We checked that the estimated sample size (ESS) for each parameter was above 200 (Nascimento et al. 2017).
Five fossil calibrations were used in our molecular dating analyses: 1) the divergence between Cetacea and Hippopotamidae was calibrated at 66-52.4 Ma based on the oldest reported crown artiodactyl, Himalayacetus subathuensis (Bajpai and Gingerich 1998;Benton et al. 2015); 2) the divergence between baleen and toothed whales was calibrated using the earliest record of a cetacean from 36.4 ± 1.5 Ma (Lambert et al. 2017); 3) the divergence between Phocoenidae (porpoises) and Monodontidae (narwhals) was calibrated at 19.50-7.50 Ma using the oldest Phocoenidae + Monodontidae record of Salumiphocoena stocktoni (Barnes 1985) and the oldest delphinidan record of Kentriodon pernix (Kellogg 1927); 4) the age of the root of Delphinidae (exclusive of Lagenorhynchus albirostris) was calibrated to 19.50-8.50 Ma based on the oldest stem Orcinus record of Eodelphinus kabatensis (Murakami et al. 2014) and the oldest delphinidan record of Kentriodon pernix (Kellogg 1927); 5) the age of the root of Delphininae (exclusive of Sotalia guianensis) was calibrated at 8.50-3.98 Ma using the oldest delphinine record of Etruridelphis giulii (Bianucci 2013) and the oldest stem Orcinus record of Eodelphinus kabatensis (Murakami et al. 2014).

Gene Flow Analyses
Seemingly conflicting or unstable phylogenetic reconstructions can be caused by gene flow among species ). Here, we tried two different types of gene flow analyses. First, we employed ABBA-BABA based analyses (Green et al. 2010;Durand et al. 2011), which are relatively simple, direct, and robust. Then, we applied the tree-based method, which tries to accommodate ILS and gene flow together to infer evolution history.
The ABBA-BABA tests measured the excess of shared polymorphisms (SNPs) of two closely related species concerning a third species (Green et al. 2010) and discriminated between gene flow and ILS. To examine possible gene flow among different species of Delphinidae, the program ANGSD (Korneliussen et al. 2014) was used for admixture analysis among all the individuals in Delphinidae using killer whale as the outgroup. All possible four-taxon topologies of the Delphinidae species including Tursiops aduncus, T. truncatus, Stenella coeruleoalba, Stenella attenuata, Delphinus delphis, Sousa chinensis, Grampus griseus, and Sa. obliquidens were involved in the gene flow analysis using D-statistics. In total, 369 combinations were used in this analysis. ANGSD allows the calculation of a Z-value based on jackknifing with 5 Mb blocks (Efron 1981) for the null hypothesis that the D-statistic is 0, which means that the existence of introgression can be rejected if the Z-value is below the significance level (an absolute value of the Z above 3 is often used as a critical value).
In addition, we analyzed the data using D FOIL -statistics  to detect the direction of introgression. This test required an asymmetric five-taxon tree with a specific topology (((P 1 ,P 2 ),(P 3 ,P 4 )), Outgroup), and (P 3 ,P 4 ) divergence occurring before the divergence of P 1 and P 2 . Therefore, not all combinations of oceanic dolphin taxa could be analyzed (though all individuals of the genera Tursiops and Stenella could be), and in total, 48 combinations were calculated. Again, the killer whale was used as the outgroup. For each combination, we used the 100 k bp non-overlap windows, which were extracted from NewGSA using a custom Python script, and SeqKit was used to extract combinations from these 100 k bp windows. Other parameters of D FOIL -statistics were left at default.
To investigate whether hybridization windows contribute to specific biological functions, after each window was identified, the corresponding sequence of the killer whale was retrieved and compared against all coding DNA sequences of killer whale using blastn, requiring an E value ≤ 10 −40 , identity ≥ 70%, and coverage ≥ 60% or better. Then, all genes from the hybridization windows were enriched to assigned GO terms (Gene Ontology Consortium 2004) for Biological Process (BP). Finally, we kept all genes that met a significance threshold of P < 0.05 after false discovery rate (FDR) correction (Benjamini and Hochberg 1995).
In the PhyloNet (Wen et al. 2018) analysis, we used the trees based on 10 k non-overlap windows. The highest individual coverage of particular species was retained from the input window trees because the phylogenetic positions of other individuals were unambiguous. Two methods were applied to infer phylogenetic networks, including maximum parsimony (MP) and maximum pseudo-likelihood (MPL). The MP and MPL models were run with 50 iterations, and the maximum number of reticulations was set to either one or two, yielding five optimal networks. All other parameters were left at default settings. Analyzing networks with more than two reticulations was too complex and not interpretable from the PhyloNet results.

Phylogenetic Relationships
A ML tree and a Bayesian tree based on the concatenated sequence alignments of ~143 Mbp revealed the following relationship with the highest bootstrap supports or posterior probabilities for all nodes (Fig. 2a): the respective monophyletic status of Cetacea and its two suborders (Mysticeti and Odontoceti), the differentiation between Delphinidae and other two families (Phocoenidae and Monodontidae) within superfamily Delphinoidea, the placement of O. orca as the earliest branching lineage in Delphinidae and the subsequent successive divergence of Sa. obliquidens, G. griseus, So. chinensis and D. delphis. In particular, both Tursiops and Stenella were found to be non-monophyletic, with T. aduncus most closely related to St. attenuata and T. truncatus most closely related to St. coeruleoalba.
In addition to the concatenation method above, which did not consider ILS, phylogenetic relationships were also inferred using the coalescent-based method. The ASTRAL tree based on a sliding window approach was congruent with the concatenated analyses (Fig. 2b). All nodes received the highest support values. Unsurprisingly, ASTRAL trees based on the three tree subsets produced identical topologies as the ASTRAL tree based on all window-based trees (Online Resource 3). SVDquartets analysis also supported the nonmonophyly of Tursiops and Stenella, although this result conflicted with the inferred topology from the RaxML, ExaBayes and ASTRAL analyses (Online Resource 4). Moreover, we summed up the topologies that supported non-monophyly or monophyly of Tursiops and Stenella. We noticed more trees supported the non-monophyletic status of Tursiops and Stenella than supported the monophyletic Fig. 2 Phylogeny of cetaceans. a The concatenation-based phylogenetic tree from NewGSA data set and 5 fossil calibrations; b Coalescentbased tree inferred from 10 k bp window-based trees by ASTRAL; c Distribution of average bootstrap support value for individual window. All nodes in a have the highest support or posterior probabilities. In a, clade letters are identical to those in Table 1. Red boxes in a at each node represent calibration points and blue boxes indicate divergence dates estimated without prior calibration constraints for that node. In a, the bounds of the boxes correspond to the 95% highest posterior density (HPD) of each node. All nodes in b have the highest support values status of both genera, which further corroborated the nonmonophyly of Tursiops and Stenella (Fig. 2c).

Divergence Time Calibration
To provide a tree-calibrated tree for the gene flow analyses, D FOIL , the ML tree, and corresponding sequence data, which can build the same topology as the ML tree, were used as input for the program MCMCTree in the software package PAML4 to estimate divergence time. All estimated divergence times for nodes with labels from A to M are presented in Fig. 2

Extensive Gene Flow among Different Dolphin Species
We investigated interspecies gene flow signals using some gene flow analyses. The D-statistics analyses, also called ABBA-BABA tests, found evidence of gene flow among species within the family Delphinidae (Fig. 3, Online Resource 5). Our ABBA-BABA tests showed that the ABBA-BABA comparisons (Z-value > 3) are significantly different from zero, suggesting that gene flow may be an important factor in generating the reticulations among these lineages. These analyses supported the view that D is significant in many of the groups (95%, 352 out of 369 groups; Online Resource 5), which revealed a complex network of ancestors among several dolphin lineages (Fig. 3) (Fig. 3). The fact that the D-statistics analyses found evidence of gene flow between most dolphin sister species may be related to their complex evolutionary history (Fig. 3).
To detect the direction of introgression between different genus pairs, we implemented the D FOIL analyses. In all 48 tests, gene flow can be detected in approximately 9 ~ 12% of individual permutations (Fig. 4), showing that gene flow among the four dolphins existed and was strong. Meanwhile, in all 48 tests, among the detected gene flow events, the number of ancestral gene flow scenarios (P 12 < = > P 3 and P 12 < = > P 4 ) without distinguished direction were higher than the number of directional gene flow scenarios (recent gene flow scenarios, e.g., P 1 = > P 3 ). Tursiops and Stenella exhibited widespread signals of ancient gene flow, likely due to their broad historical ranges throughout the ocean, overlapping with several congeneric species.
To test whether the introgressed windows in Delphinidae were related to fitness or adaptation, the 314 introgressed windows identified by D FOIL analyses were compared against all killer whale coding-protein sequences and were found to contain notable matches to 6,378 genes. Of these 6,378 genes, 167 genes (e.g., PBX1, IL-17F) were associated with "skeletal system development" (GO: 0001501, P = 0.00466 < 0.01), 129 genes (e.g., BARHL2, RBP4) were associated with "regulation of developmental growth" (GO: 0048638, P = 2.07e −05 < 0.01), and 106 genes (e.g., PTPRS, PAK1) were associated with the GO pathway "developmental growth involved in morphogenesis" (GO: 0060560, P = 1.74e −07 < 0.01) (Table 2, Online Resource 6). In addition, among 6,378 genes, 134 genes were related to "visual system development" (GO: 0150063, P = 0.000970 < 0.01) and 130 genes were related to "response to drug" (GO: 0042493, P = 0.00259 < 0.01) (Online Resource 6). These results might imply that differential introgression of genes gave rise to or contributed to similar morphological features and fitness in dolphins, as suggested for Neanderthal and Denisovan alleles in modern humans (Racimo et al. 2015). For example, PBX1 plays a crucial role in osteogenesis and matrix mineralization (Cheung et al. 2009). Although our analysis did not explore the specific role of these introgression genes, it did provide pertinent evolutionary hypotheses for future investigation. In addition to the character-based analysis, gene flow may preferably be studied by tree-based phylogenetic networks using both MP and MPL models. Although these network topologies differed from the concatenation and coalescentbased method, both pseudolikelihood and parsimony networks with a maximum of two allowed reticulation supported the non-monophyletic status of Tursiops and Stenella. In both parsimony networks, St. attenuata was connected by a hybrid edge to the lineage of T. truncatus (Fig. 5a-b). The inferred inheritance probabilities for the hybrid edge differed slightly between the networks with a maximum of one and two allowed reticulations. When we allowed for a The results of 48 tests, which were based on 100 k windows maximum of two reticulations in the pseudolikelihood networks, the connection between T. aduncus and the ancestor of T. truncatus and St. coeruleoalba was observed (Fig. 5d). Interestingly, T. aduncus had a reticulation edge to the stem lineage of T. truncatus and St. coeruleoalba, resembling the result of the D-statistics analyses, which meant that some loci in T. truncatus or St. coeruleoalba shared a most recent common ancestor with loci from T. aduncus.

Non-monophyly of Tursiops and Stenella
The phylogenetic relationships of the family Delphinidae have been notoriously difficult to disentangle due to hybridization and ILS arising from recent explosive speciation (Wiens et al. 2006;Whitfield and Lockhart 2007;McGowen et al. 2020). Although a number of morphological and genetic studies have been conducted in an attempt to resolve their evolutionary relationships (e.g., Perrin et al. 1987;Archer and Perrin 1999;Chen et al. 2011;Moura et al. 2013Moura et al. , 2020Jedensjö et al. 2017), many aspects of Delphinidae phylogeny remain unresolved, especially the monophyly of Tursiops and Stenella. Traditionally, morphologybased taxonomy has considered Tursiops and Stenella to be monophyletic groups, based on features of the skull and superficial morphological observations, respectively (Archer and Perrin 1999;Shirakihara et al. 2003). However, Leduc et al. (1999) published a phylogeny for the 33 thenrecognized Delphinidae species based on cytochrome b, which revealed that T. truncatus formed a sister group with a polytomy formed by St. coeruleoalba + Stenella clymene, Delphinus spp., Stenella frontalis, and T. aduncus. Their results challenged the monophyly of Tursiops and Stenella for the first time. Since then, some researchers also have suggested that Tursiops is a non-monophyletic group (e.g., Kingston et al. 2009;Zurano et al. 2019). For example, Moura et al. (2013) showed T. aduncus clustering with Delphinus capensis based on the whole mitochondrial genome, while McGowen (2011) found a sister relationship between T. aduncus and T. truncatus + St. coeruleoalba, based on 20 nuclear genes. However, in another study based on 3,191 protein-coding genes, McGowen et al. (2020) suggested that T. aduncus and T. truncatus form a monophyletic group that constitutes a sister group within Delphinidae. For Stenella, Leduc et al. (1999) suggested Stenella may be an artificial assemblage, as members of this genus are both morphologically and genetically very dissimilar. The genus Stenella has been rendered non-monophyletic in nearly all analyses since then (e.g., Nishida et al. 2007;Chen et al. 2011;Dornburg et al. 2012). However, due to the considerable phylogenetic uncertainty among the species within Stenella, a formal redescription has not been attempted (McGowen et al. 2020). For example, Steeman et al. (2009) showed St. coeruleoalba and St. frontalis to be more closely related to members of Delphinus than to their congeners, based on six mitochondrial and nine nuclear genes; moreover, Zurano et al. (2019) showed that St. coeruleoalba clustered with St. clymene, and that St. frontalis was sister to the clade including D. capensis + D. delphis and T. aduncus, using a whole mitochondrial genome. Differences in molecular markers, taxon sampling, and phylogenetic inference methods can interfere with resolving phylogenetic relationships and result in contrasting conclusions. For this reason, we tried to conduct phylogenetic analyses with large-scale genomic sequence data in this study, avoiding the impact of different markers, and provide a novel opportunity to address the problem of Tursiops and Stenella monophyly.
In the present study, phylogenetic reconstruction was performed based on concatenation and coalescent-based approaches to infer the relationships within Delphinidae. Consistency between these results based on different treebuilding methods should be viewed as better support for a hypothesis than bootstrap values from a single analysis (Suh 2016). Trees derived from the concatenation and coalescentbased method (ASTRAL) recovered a T. truncatus + St. coeruleoalba clade and a T. aduncus + St. attenuata clade, both of which were supported by the highest support value ( Fig. 2a-b). This non-monophyletic result coincided with previous results based on whole mitogenomes (Moura et al. 2013;Zurano et al. 2019). Furthermore, the window tree statistics showed approximately five times as many window trees supporting the non-monophyly of Tursiops as window trees that recovered Tursiops as a monophyletic lineage; similarly, window trees overwhelmingly provided support for the non-monophyly of Stenella (Fig. 2c). These results implied that the genomic fragments supporting Tursiops and Stenella as non-monophyletic groups are widely distributed throughout the whole genome. In addition, two phylogenetic network topologies supported the non-monophyletic status of Tursiops and Stenella (Fig. 5b, d). Taking all these evidences together, based on different phylogenetic analyses and the largest genomic data so far, our analyses provided novel strong support for the non-monophyly of Tursiops and Stenella and constitute the best current non-monophyletic hypothesis for the genera Tursiops and Stenella. Early molecular studies based on a few nuclear DNA or the whole mitochondrial genomes showed that, with increased taxon sampling, the genera Tursiops and Stenella are recovered as non-monophyletic groups although internal nodes often have < 70% bootstrap support (e.g., Kingston et al. 2009;Hassanin et al. 2012). In agreement with previous studies (Kingston et al. 2009;Hassanin et al. 2012), we resolved a non-monophyletic Tursiops and Stenella; in particular, we recovered strong support values for all internal nodes. Nevertheless, some phylogenomic studies have supported T. aduncus and T. truncatus as a monophyletic group within Delphininae (e.g., McGowen et al. 2020;Moura et al. 2020). Based on the Bayesian inference method, Moura et al. (2020) used RADseq data including nearly all delphinine species and recovered Tursiops as a monophyletic clade with 0.92 posterior probability. However, it should be noted that the Bayesian inference method consistently overestimates support and therefore should not be interpreted as probabilities that clades are correctly resolved (Simmons et al. 2004). In particular, compared with Moura et al. (2020), phylogenetic trees inferred from concatenation (ML and Bayesian) and coalescent-based approaches were concordant with each other in the present study, which supported non-monophyly of Tursiops. Based on the largest genetic assessment of the Tursiops species to date, we showed that the current taxonomic framework of Tursiops greatly underestimates its diversity. Previously, Charlton-Robb et al. (2011) found Tursiops to be non-monophyletic and named Tursiops australis sp. nov. with the common name of 'Burrunan Dolphin', based on morphological and genetic differentiation; Gray et al. (2018) also discovered a new 'aduncus' type lineage in the Arabian Sea. Considering that traditional species groupings of dolphins were based largely on the skull or superficial morphological observations (Perrin et al. 2013), it is not surprising that the taxonomy of Tursiops and Stenella has been surrounded by uncertainty and indicates gene flow in their evolutionary history (see the next section). In addition, we are aware that poor taxon sampling can result in an ambiguous result (Philippe et al. 2011); therefore, the current non-monophyletic status of the genera Tursiops and Stenella still needs further validation when more genome data for Tursiops and Stenella are available.

Widespread Gene Flow among Oceanic Dolphins
Hybridization might produce phylogenetic discordance in genomes, which has been especially observed in species undergoing rapid speciation such as big cats, butterflies, and cichlids (Figueiró et al. 2017;Malinsky et al. 2018;Edelman et al. 2019). Oceanic dolphins are iconic mammals with a complex evolutionary history due to a recent and rapid diversification process (Steeman et al. 2009;Zhou et al. 2011). A large number of phylogenetic discordances extensively reported in dolphins are commonly explained by rapid radiation associated with ILS or hybridization among lineages (e.g., Amaral et al. 2012;Perrin et al. 2013). However, only very few studies have detected ILS or hybridization in dolphins with molecular data (Kingston et al. 2009;Moura et al. 2020). For example, Kingston et al. (2009) applied the STRU CTU RE program (Pritchard et al. 2000;Falush et al. 2003) to 418 polymorphic markers and indicated some hybridization events between the two spotted dolphin species. However, the potential hybridization events outside these two species may be hidden because not all of their studied species were tested. Instead, it is necessary to take more combinations among dolphins into consideration to accurately explore hybridization events and their effects on dolphin evolution.
In this study, we applied different statistical methods, i.e., D-statistic, D FOIL and PhyloNet analyses, to illuminate the existence of gene flow in dolphins. The D-statistic found remarkable gene flow among all members of Delphinidae (Fig. 3). Although only one So. chinensis sample was included, the Z-values > 3 of most combinations including So. chinensis (74/75, Online Resource 5) were massive. Besides, only one individual of St. attenuata was sampled.
Of course, small sample size may cause biased results due to SNP (Hoggart et al. 2008) and thus, the current gene flow signals needed to be further confirmed in the future. Furthermore, considering the wide overlapping geographical distribution among different dolphin species (Jefferson and LeDuc 2018), as well as several studies showing successful mating between different dolphin species (Crossman et al. 2016;Gridley et al. 2018), it has been regarded that extensive gene flow could characterize all oceanic dolphins rather than just closely related species, and it may still be ongoing. Furthermore, ancestral gene flow is expected to leave signatures in genomes of their descendants, causing some phylogenetic conflict (Kumar et al. 2017). The instance of ancestral gene flow between T. aduncus and the ancestor of T. truncatus and St. coeruleoalba (Figs. 3 and 5d), might be the factor leading previous studies to infer an incorrect relationship between T. aduncus and T. truncatus in previous cases (e.g., McGowen et al. 2009;Amaral et al. 2012). Although it was inferred in the present study that gene flow might be responsible for part of the non-concordance within the dolphin evolutionary tree because the D values calculated by D-statistic and D FOIL were significantly different from zero (Figs. 3 and 4), phylogenetic discordance might also be attributed to ILS (Whitfield and Lockhart 2007). Therefore, future studies are required to further identify the role of these two processes. Based on our data, our study strongly suggested pervasive gene flow in oceanic dolphins, and this may help explain why previous phylogenetic studies with limited markers (e.g., Xiong et al. 2009;McGowen 2011;Amaral et al. 2012) were unable to reliably resolve the phylogenetic relationships of this highly rapidly radiated group of animals.

Conclusion
In this study, we produced a high-resolution phylogenetic reconstruction for dolphins with a special focus on addressing the monophyly hypothesis of Tursiops and Stenella. Based on the genomic sequences, we were able to obtain a well-resolved tree that revealed the power of multiple complete genomes from closely related species to comprehensively infer their evolutionary history. Based on the concordant phylogenetic relationships in the concatenation and coalescent-based analysis, the non-monophyletic status of both Tursiops and Stenella was strongly suggested, which implied that the current taxonomy of both genera might not reflect their evolutionary history and might underestimate their diversity. Extensive inter-specific gene flow among species detected with different methods could explain interspecific convergence among species from different genera and also instability in resolving phylogenetic relationships of oceanic dolphins with different and limited markers.

Declarations
Ethics Approval Not applicable.

Conflicts of Interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.