Specimens of S. oualaniensis were collected from the Eastern Indian Ocean (10 individuals) and South China Sea (three individuals) by using a tuck net ( Figure 1 ; Table S2 ). Five grams of muscle tissue from each individual was removed and preserved in 95% ethanol. All capture and sampling of squid individuals were undertaken in accordance with an Animal Research ethics permit (W20211202) granted by the East China Normal University Animal Ethics Committee.

Total genomic DNA was extracted from each specimen using a standard phenol–chloroform extraction method (Zhang and Hewitt, 1998). The two mitochondrial DNA segments cytochrome b (Cyt b) and cytochrome c oxidase subunit 1 (COI) genes have proven to be good or very good at recovering interspecific and intraspecific phylogenetic relationships or population structure (Miya and Nishida, 2000). Thus, Cyt b and COI were amplified using the following two pairs of primers through polymerase chain reactions (PCRs): CytbF (5′-GAACCTCGACAMCCACCAACATG-3′) and CytbR (5′-GCTCTTGTAGATTTACCACCTGCTCC-3′) for Cyt b (Liu et al., 2017); COIF (5′-GTATTTGGGCAGGACTTCTAGGGA-3′) and COIR (5′-TAATGGCAAATATGGCGTGATCATGG-3′) for COI (Li et al., 2019). Amplifications were conducted in 50-μl volumes containing 5 μl of template DNA, 25 μl of 2× EasyTaq PCR SuperMix, 200 nM each of the appropriate primers (CytbF/CytbR, and COIF/COIR), and 18 μl of ddH₂O. PCR was performed in the following cycle parameters: 1) for Cyt b, denaturation at 94°C for 4 min, followed by 35 cycles of 94°C for 30 s, 60°C for 30 s, 72°C for 60 s, and a final extension of 72°C for 10 min; 2) for COI, denaturation at 94°C for 4 min, followed by 35 cycles of 94°C for 30 s, 58°C for 30 s, 72°C for 60 s, and a final extension of 72°C for 10 min. Amplified products were visualized on a 1.5% agarose gel and then purified with the Gel Extraction Mini Kit (Watson BioTechnologies, Shanghai, China). Purified DNA products were sequenced in both directions using primer pairs of CytbF/CytbR and COIF/COIR on an ABI Prism 3730 automatic sequencer. All nucleotide sequences were deposited in GenBank (Cyt b, MZ735375–MZ735387; COI, MZ712955–MZ712967). To infer the evolutionary history of S. oualaniensis, these sequences were supplemented with Cyt b and COI sequences obtained from three complete mitochondrial genomes of S. oualaniensis (EU658923, EU660576, and EU660577) from the Central-East Pacific Ocean. Intraspecific phylogenetic relationships of the sibling species Sthenoteuthis pteropus were also reconstructed based on COI sequences including MW817786–MW817789, OM988016, OM988017, MG591279, MF980608, MF980594, MF980595, GU145073, MW255556, and AB270942 sampled from the Atlantic Ocean ( Figure S1 ). Furthermore, sequences from S. pteropus (MW255556) and an additional squid species, Dosidicus gigas (EU068697), were utilized as outgroups in the following phylogenetic tree reconstructions of S. oualaniensis. All haplotypes and sampling information of S. oualaniensis, S. pteropus, and D. gigas used in this study are summarized in Table S2 .

Cyt b and COI sequences amplified from the same individual were concatenated and aligned in Cluster X (Thompson et al., 1997). Modeltest 3.6 (Posada and Crandall, 1998) was used to select the appropriate model of evolution for the COI (TrN+I), Cyt b (TIM+I), and concatenated (TrN+I) datasets by the Akaike information criterion (AIC). The model parameters of three datasets were listed as follows: 1) COI, base frequencies with A (0.2993), C (0.2124), and G (0.1253); rate matrix with R [A–C] = 1.0000, R [A–G] = 90.4377, R [A–T] = 1.0000, R [C–G] = 1.0000, R [C–T] = 13.3904, and R [G–T] = 1.0000; proportion of invariable sites (I) = 0.7282. 2) Cyt b, base frequencies with A (0.4438), C (0.2013), and G (0.1220); rate matrix with R [A–C] = 1.0000, R [A–G] = 63.8296, R [A–T] = 0.0106, R [C–G] = 0.0106, R [C–T] = 21.7894, and R [G–T] = 1.0000; proportion of invariable sites (I) = 0.7300. 3) Concatenation, base frequencies with A (0.3539), C (0.2086), and G (0.1236); rate matrix with R [A–C] = 1.0000, R [A–G] = 109.4730, R [A–T] = 1.0000, R [C–G] = 1.0000, R [C–T] = 20.6656, and R [G–T] = 1.0000; proportion of invariable sites (I) = 0.7322.

Phylogenetic trees were inferred using the maximum likelihood (ML) (Felsenstein, 1981) methodology in PHYML 3.0 (Guindon et al., 2010). Nodal supports of the ML tree were estimated using bootstrap values (Felsenstein, 1985) with 10,000 pseudo-repetitions. A Bayesian inference (BI) (Huelsenbeck and Ronquist, 2001) tree was further recovered in MrBayes 3.2 (Ronquist et al., 2012) with an evolutionary model (TrN+I). The model parameters of BI tree were estimated in the program based on BIC criteria using the following settings: ngen = 25,000,000, samplefreq = 100, burnin = 62,500, Nchains = 4, and Nruns = 2. The convergence of independent runs was achieved when white noise was seen in the overlay plot of generation versus the log-likelihood probability for both runs with the potential scale reduction factor (PSRF) approaching 1, a lower average standard deviation of split frequencies (0.000871 < 0.01). The nodal stability of the BI tree was assessed using posterior probability (PP) with a 60% consensus rule. Moreover, median-joining (MJ) network (Bandelt et al., 1999) was used to clarify the intraspecific and interspecific evolutionary relationships of S. oualaniensis and S. pteropus based on concatenation data set and COI segment in NETWORK (version 10200; fluxusengineering, 2021), respectively. A single haplotype from each clade of S. oualaniensis was added to COI sequences of S. pteropus to overcome questions of insufficient haplotypes in network reconstruction.

Species trees for COI and Cyt b data sets were also estimated in BEAST v.1.8.2 package (Drummond and Rambaut, 2007). Unlinked models of sequence evolution of each dataset were implemented using BEAUTi (TrN+I for COI, a similar GTR+I to TIM+I for Cyt b). A strict clock model and the species tree prior to a Yule Process were used. Four independent analyses with the same parameters were performed using Markov chain Monte Carlo (MCMC) simulations for a length of chain 100 million and a sampling frequency of every 10,000 generations in BEAST. Convergences of these analyses were inspected using Tracer 1.5 (Rambaut and Drummond, 2009), and the first 10% of the trees were discarded as burn-in. A single maximum clade credibility tree was reconstructed using TreeAnnotator v1.8.2 (Drummond and Rambaut, 2007) and FigTree v. 1.4.4 (Rambaut, 2009).

In order to distinguish ambiguous internal phylogenetic relationships among four clades from the northern East Indian Ocean, southern East Indian Ocean, South China Sea, and Central-Southeast Pacific, a likelihood mapping analysis was performed in TREE-PUZZLE version 5.2 (Schmidt et al., 2002). Likelihood mapping triangle diagrams were used to compare the competing relationship between the four clusters and to assess the support values of internal branches (Strimmer and von Haeseler, 1996; Strimmer and von Haeseler, 1997).

Pairwise net genetic distances between and within different lineages were inferred using D_net = D_xy − 1/2 * (D_x + D_y) (Nei, 1987) in MEGA X (Kumar et al., 2018), where D_xy, D_x, and D_y are the average numbers of nucleotide difference based on the Kimura 2 parameters correction (Kimura, 1980) between or within clades x and y (Nei and Li, 1979; Nei, 1987).

The divergence time between each pair of the S. oualaniensis lineages was estimated in BEAST v.1.8.2 based on the concatenation dataset. The TrN+I evolution model was used with empirically determined base frequencies. For the tree prior parameter, the Speciation : Birth–Death process (Gernhard, 2008) was used and modeled with an uncorrelated lognormal relaxed clock (Drummond et al., 2006). The following calibration node was constrained using lognormal distribution priors: 28.43 million years ago (Ma), coalescent time of S. oualaniensis and S. pteropus based on four fossils calibrated points (Pardo-Gandarillas et al., 2018). Three independent MCMC runs of 100 million generations with a sampling frequency of every 10,000 generations were performed. The log files were combined in Tracer to summarize posterior divergence times and 95% highest posterior density limits. Only runs where the ESS values over 200 for all marginal parameters were used, after discarding 10% of the first trees as burn-in.

The Bayesian binary method (BBM) implemented in RASA 4 (Yu et al., 2020) was used to reconstruct the ancestral area. Biogeographical data are based on the present distribution of S. oualaniensis and S. pteropus ( Figure 1 ; Figure S1 ). The BBM analysis was run for 50 million cycles with 10 chains and a sampling frequency of every 1,000 cycles. The temperature was set at 0.1, and a fixed JC model was used. The maximum number of areas for each node was four. The information of each node was summarized and plotted as a pie chart.

Similar topologies were recovered using different tree-building methodologies under the concatanated data set ( Figures 2 and S2 ) and species tree reconstruction ( Figure 3 ). Three new clades (I, II, and III) of S. oualaniensis, with reliable nodal supports (BP/PP > 90/0.98), were consistently recovered by both of the phylogenetic methods (ML/BI). Clade I from the South China Sea is composed of two subclades diverging at ca. 0.27 (0.90–0.05) Ma (D_net = 0.21%) ( Table 1 ). Clade II and clade III distribute in the northern East Indian Ocean and southern East Indian Ocean, respectively. Based on the branching order, two basal clades (V and VI) were consistently observed with strong statistical support (BP/PP = 84/0.97, 100/1.0) ( Figures 2 and S2 ). Clade VI/Equatorial, from the central Equatorial Pacific, was recovered to be basal to all other clades of S. oualaniensis coalescing at 14.27 (22.83–5.95) Ma (D_net = 10.67%). Clade V/Eastern Typical from the northern East Pacific is the second most basal clade diverging from the remaining lineages (I–IV) at 3.40 (7.38–0.83) Ma (D_net = 5.16%). A monophyletic East Indian Ocean branch grouping the two East Indian Ocean clades (II and III, D_net = 1.99%) at 0.88 (1.9–0.25) Ma was also successfully recovered in ML and BI trees with low statistical support (BP/PP = 0.45/0.72) ( Figures 2 and S2 ). Furthermore, this ancestral East Indian Ocean branch and clade IV/Pacific Typical coalesced at ca. 1.06 (2.48–0.31) Ma with medium-to-strong support (BP/PP = 0.69/0.97). The divergence time of clade I and clades II–IV was inferred to be 1.3 (3.0–0.4) Ma. The species tree reconstruction provided a similar topology to the ML and BI trees: among the four clades (I–IV), clade I was closely related to clade IV ( Figure 3 ).

In order to further resolve the ambiguous phylogenetic relationships among clades I–IV, a quartet topology was reconstructed based on the concatenation dataset ( Figure 4 ). The ancestral branch, grouping southern and northern East Indian Ocean clades (II and III) in Figure 4A , was strongly supported with the highest posterior weights (88.8% quartets; Figure 4D ) relative to two other trees ( Figures 4B, C ). Thus, this ancestral East Indian Ocean branch was confirmed with reliable statistical support.

The six lineages (clades I–VI) of S. oualaniensis with different divergence levels were also recovered in the network ( Figure S3 ). However, relative to the deep divergence observed in S. oualaniensis, only one clade with two closely related haplotypes (SP1 and SP2) of S. pteropus was revealed ( Figure S4 ).

The BBM revealed eight events of dispersal and seven events of vicariance along the phylogeny of D. gigas and Sthenoteuthis ( Figure 5 ). Estimation of ancestral locality for the basal node (39) indicated that the recent common ancestor of D. gigas and Sthenoteuthis likely originated from the East Pacific (Southeast Pacific, PP = 0.47; northern East Pacific, PP = 0.31) at ca. 38 Ma. Node 38 showed that the Atlantic Ocean has a higher probability (PP = 0.57) to be the ancestral distribution of Sthenoteuthis relative to the northern East Pacific (PP = 0.20). Node 37 showed a mixed ancestral origin of S. oualaniensis in either the northern East Pacific (PP = 0.41) or central Equatorial Pacific (PP = 0.40). Node 36 indicated an ancestral distribution of clades I–V in the northern East Pacific (PP = 0.80). For the ancestral node 35 of clades I–IV, the South China Sea origin was indicated (PP = 0.58). Node 30 showed an ancestral origin of clades II–IV in the Central-Southeast Pacific (PP = 0.68). Node 29 suggested that the Indian Ocean ancestral population of clades II–III likely originated from the southern East Indian Ocean (PP = 0.54).

The deepest genetic divergence (D_net = 10.67%) based on our concatenated COI and Cyt b datasets was observed between clade VI (dwarf form) from the central Equatorial Pacific and all other clades (I–V) ( Figures 2 , S3 and Table 1 ), which is similar to the divergence found between dwarf and middle-sized forms in previous studies (Li et al., 2019; Xu et al., 2020). The second deepest genetic divergence (D_net = 5.16%) between clade V and clades I–IV is also similar to that of middle-sized forms between the East Pacific Ocean and South China Sea based on COI gene sequences (4.85%; Xu et al., 2020). These inter-clade divergences between clade VI and clades I–V or clade V and clades I–IV of S. oualaniensis are much higher than the mean level of mollusk intraspecific distance (<3%) (Yu et al., 2016). In general, a standard threshold greater than a 10-fold difference between interspecific and intraspecific genetic divergences has been used for species identification (Hebert et al., 2003; Hebert et al.,2004; Ward et al., 2005). Following this criterion, the dwarf form of clade VI found in the central Equatorial Pacific and the middle-sized form of clade V from the northern East Pacific may represent separate species or subspecies (Clarke, 1966; Li et al., 2019; Xu et al., 2020). Although no dwarf forms were presently sampled from the East Indian Ocean and the South China Sea, this form has been observed in the South China Sea (Li et al., 2019; Xu et al., 2020), suggesting ongoing transoceanic gene flow between the central Equatorial Pacific and the West Pacific.

Two morphological populations of S. oualaniensis, viz., dwarf form and middle-sized form, have been identified in the South China Sea (Liu et al., 2019; Wang et al., 2019). Based on Cyt b gene sequences, Liu et al. (2017) found two genetic lineages of S. oualaniensis (minor A and dominant B) inhabiting the South China, East China, and Philippine Seas. The genetic divergence (>1.67%) found between the two lineages (A and B) corresponds to the inter-clade divergence levels between I, II, III, and IV (1.41%–2.16%) of the middle-sized form found here ( Table S3 ). Although no published data are available for comparison, the dominant cluster B could be analogous to clade I from the South China Sea in this study inferred from its high population abundance and sampled probability. Inferred from COI and H3 gene, the dwarf form and the dominant middle-sized form observed from the South China Sea by Xu et al. (2020) correspond well to clade VI and clade I in this study, respectively. Thus, at least three lineages including middle-sized form A, middle-sized form B/clade I, and dwarf form/clade VI have been revealed in the South China Sea, which also provide direct evidence for the three populations that were identified based on dorsal photophores by Jiang et al. (2015) (i.e., the middle-sized, dwarf-sized, and X3-sized forms). Two subclades of clade I with reliable nodal supports were further revealed from our limited sampling in the South China Sea ( Figures 2 , S2, S3 ), which may be indicative of high diversity and frequent population/lineages differentiation of S. oualaniensis in the Indo-West Pacific.

A divergence level of 1.99% was observed between clade II from the northern East Indian Ocean and clade III from the southern East Indian Ocean, which is far lower than that (10.33%–10.81%) observed between the dwarf form (clade VI) and middle-sized form (clades I–V), suggesting that both of these East Indian Ocean clades likely belong to different lineages of the middle-sized form in the Indian Ocean. Although the dominant middle-sized form and the dwarf form were found north of 10°N and south of 7°N along the neighboring southwest coast of India (Chembian and Mathew, 2014), no dwarf form was observed in the present sampling from the East Indian Ocean. Further studies will be needed to further distinguish these populations including the three middle-sized forms in the Indo-Pacific and the giant form in the Arabian Sea (Nesis, 1993; Snÿder, 1998) through integrating molecular and morphological data on a larger spatial scale.

Compared to deeply divergent lineages observed in S. oualaniensis from the Indo-Pacific, a single lineage was observed in S. pteropus from the Atlantic Ocean ( Figure S4 ). Three possibilities can be invoked to explain the low divergence of S. pteropus: 1) recent range expansion and ongoing high gene flow transported by circulation likely have caused a homogeneous population structure of S. pteropus (He et al., 2010). 2) Population bottlenecks likely have erased other lineages of S. pteropus due to historical disaster events (Hewitt, 1996). 3) Insufficient sampling likely caused to artificially underestimate diversity of S. pteropus (Yan and Zhang, 2004). Further sampling on a larger spatial scale and more sample size will be helpful to uncover the population structure of S. pteropus.

Sthenoteuthis is thought to have originated in the Atlantic Ocean and subsequently migrated eastward into the Indian and Pacific Oceans (Zuyev et al., 1975; Zuyev et al., 2002). During the late Pliocene and early Pleistocene, low sea levels are hypothesized to have led to the genetic differentiation and speciation of Sthenoteuthis due to isolation between the Atlantic and Indo-Pacific Ocean.

Present molecular data, fossil records, simulated paleo-ocean circulations, and geologic clues partly support the abovementioned evolutionary origin of Sthenoteuthis. Biogeographical analysis based on a concatenated dataset suggested that the common ancestor of D. gigas and Sthenoteuthis likely originated from the East Pacific with a high probability of 0.78 (Central-Southeast Pacific, PP = 0.47; northern East Pacific, PP = 0.31) at ca. 38 Ma ( Figure 5 ), and the ancestral origin of Sthenoteuthis was further inferred to be from the Atlantic Ocean with a probability of 0.57. Thus, the common ancestor of D. gigas and Sthenoteuthis was inferred to disperse from the Southeast Pacific to the Atlantic Ocean, and the ancestral population of Sthenoteuthis subsequently dispersed back to the northern East Pacific from the Atlantic Ocean.

The ancestral origin of D. gigas and Sthenoteuthis in the Pacific Ocean was also consistent with fossil clues of Ommastrephidae. The mean coalescent time of Ommastrephidae has been inferred to have occurred during the late Cretaceous (87.1 Ma; Pardo-Gandarillas et al., 2018), which preceded the earliest known Ommastrephidae fossils from the coastal northeast Atlantic Ocean, which are from the Middle Eocene (47.8–37.8 Ma) (Neige et al., 2016). However, the inferred coalescent age of Ommastrephidae is consistent with the lower jaw fossil of Teuthida from the Upper Cretaceous strata in Hokkaido, Japan (Tanabe et al., 2015). A large amount of the earliest known fossils and biodiversity hotspots of marine species including mangrove, corals, mollusk, and coral reef fish assemblages were observed in the proto-Pacific and Tethys Oceans (Piccoli et al., 1987; Bellwood, 1996; Renema et al., 2008; Neige et al., 2016). The common ancestor of D. gigas and Sthenoteuthis likely dispersed from the East Pacific to the Atlantic Ocean through the open Central American Seaway at 48–28 Ma followed by genetic divergence and speciation due to vicariance between the Pacific Ocean and Atlantic Ocean (Scotese, 2004).

Surface sea currents before 5 Ma promoted promoted the transoceanic migration of marine organisms between the Pacific Ocean and Atlantic Ocean. The North Equatorial Countercurrent extended eastward from the Pacific Ocean to the Caribbean Sea through the southern Panama Seaway during the Oligocene and Miocene (Maier-Reimer et al., 1990; Nisancioglu et al., 2003; Haug et al., 2004; von der Heydt and Dijkstra, 2006; Schneider and Schmittner, 2006; Sepulchre et al., 2014) ( Figure 6A ). However, the surface water of the upper 75 m flowed westward from the Atlantic Ocean to the East Pacific through the northern part of the open Central American Seaway. These bidirectional circulations might have carried larvae of marine organisms between the Pacific Ocean and Atlantic Ocean. Some transisthmus sister species pairs of octopus, hake (Merluccius spp.), and octopuses (Muusoctopus) from the East Pacific and Atlantic Oceans also shared a recent common ancestor from the East Pacific Ocean or the Atlantic Ocean (Grant and Leslie, 2001; Ibáñez et al., 2016; Lima et al., 2020).

Lastly, the geologic history of the Central American Seaway may have also influenced speciation processes between the East Pacific and Atlantic oceans. The divergence times of D. gigas × Sthenoteuthis and S. oualaniensis × S. pteropus were traced to ca. 38 (48–28) Ma and 28 (40–17 Ma), respectively (Pardo-Gandarillas et al., 2018) ( Figure 2 ). The western arc segment of Central America began a counterclockwise rotation between 38 and 28 Ma followed by a clockwise rotation of the central arc segment between 28 and 25 Ma, and the Central American Seaway narrowed ca. 25 Ma due to tectonic collision between Panama and South America (Farris et al., 2011; Montes et al., 2012). Moreover, an abrupt sea level regression of over 150 m occurred 35–30 Ma (Haq et al., 1987), which would have caused the exposure of shallow seas (Schaaf, 1996) and restricted water mass exchange between the East Pacific and West Atlantic Ocean preceding the near closure of the Central American Seaway (Montes et al., 2012). Low extinction and peak speciation rates have been observed in marine animals including neritic mollusks during the early Oligocene (McGowran, 1990). D. gigas × ancestral Sthenoteuthis and S. oualaniensis × S. pteropus were likely separated from each other due to habitat fragmentation and reduced gene flow (Hewitt, 1996) between the Pacific Ocean and the Atlantic Ocean during Oligocene. However, ancestral Sthenoteuthis origin should further be confirmed using some earlier fossils from the Atlantic Ocean because fossil statoliths from an extinct, potentially ancestral species (Symplectoteuthis pedroensis) of genus Sthenoteuthis have been found from Pliocene deposits in a coastal East Pacific locality (California, USA) (Clarke and Fitch, 1979).

Different from the aforementioned eastward migration pathway of ancestral Sthenoteuthis from the Atlantic Ocean to the Pacific Ocean through the Indian Ocean (Zuyev et al., 1975; Zuyev et al., 2002), multiple westward dispersal events were indicated by present molecular data ( Figure 5 ): ancestral population of S. oualaniensis was inferred to firstly colonize the northern East Pacific from the Atlantic through the open Panama Seaway during late Oligocene (28–23 Ma), and the northern East Pacific further acted as dispersal origin of S. oualaniensis to central Equatorial Pacific and the South China Sea. Moreover, the ancestral Indian Ocean population was also inferred to have originated from Central-Southeast Pacific (clade IV) based on phylogenetic topologies and biogeographical analysis ( Figures 2 , 5 ). Westward migrations from the East Pacific to West Pacific and from the West Pacific to the East Indian Ocean were also observed in different marine species including fringehead blennies in the genus Neoclinus (Hongjamrassilp et al., 2020) and intertidal barnacle (Tetraclita squamosa) (Chan et al., 2022), respectively.

These inferred westward colonization routes correspond well to the paleo-circulation directions in the Indo-Pacific Ocean during the late Oligocene–Miocene (Kennett et al., 1985) ( Figure 6A ). Westward warm North Pacific Equatorial Current likely transported larvae of ancestral S. oualaniensis from the northern East Pacific to the Hawaii seamount or atolls. The population in the central Equatorial Pacific (clade VI) was firstly isolated from the remaining population (clades I–V) at ca. 14.3 (22.8–6.0) Ma. Similar population differentiation or speciation around the Hawaii island chain was also observed in other marine organisms, e.g., Panulirus marginatus (Pollock, 1992). These Hawaii seamounts/island chains could have served as a barrier to surface circulation and impeded larval gene flow between the eastern and western populations due to sea level regression during 15–6 Ma (Pollock, 1992; Woodruff, 2003).

Clade V from the northern East Pacific was subsequently separated from other clades (I–IV) at 3.4 (7.4–0.82) Ma ( Figures 2 ; Tables 1 ), which indicates a possible evolutionary influence of the East Pacific Barrier (Ekman, 1953; Briggs, 1974) and seaway closure in the Pacific (Kennett et al., 1985). Panama Gateway was closed at ca. 3.5–3 Ma (O’Dea et al., 2016), which interrupted the westward water mass from the Atlantic Ocean to the East Pacific. With moving northward of the Australian Plate, the Indonesian Gateway became narrower than before and impeded South Pacific Equatorial Circulation into the Indian Ocean at ca. 3 Ma (Cane and Molnar, 2001). These seaways’ closures strengthened gyral circulations in the Indo-Pacific (Supplementary Material). As an endemic lineage, larvae of clade V were likely retained in the region by a counterclockwise circulation along the coast of Central America and could not drift into the westward North Equatorial Current ( Figure 6B ).

Clade I in the South China Sea was inferred to derive from the ancestral population in the northern East Pacific ( Figure 5 ). The westward North Equatorial Current likely transported the ancestral population of clades I–V from the northern East Pacific to the West Pacific and South China Sea ( Figure 6A ). Furthermore, transoceanic larval exchange between the northern East Pacific and the South China Sea continued till the closure of the Panama and Indonesian seaways at ca. 3.4 Ma ( Figure 2 ). Eastward Equatorial Countercurrent and Equatorial Undercurrent had been strengthened since these seaways’ closure in the Pacific (Kennett et al., 1985), which would be helpful to carry larvae of S. oualaniensis from the South China Sea to central Equatorial Pacific and Southeast Pacific ( Figure 5 ). Marine species including crustaceans, fishes, and molluscan fauna in the South China Sea or Northwest Pacific could be transported to the central Equatorial Pacific and Southeast Pacific through the North Pacific Current, North Equatorial Current, and Equatorial Countercurrent or through southward coastal current along the east of the Philippines and the Equatorial Countercurrent (Randall et al., 1985; Scheltema, 1988; Ding and Bai, 1993; He et al., 2021). These alternative migration pathways can be distinguished through future sampling in a larger spatial scale from the Philippines Sea, central Equatorial Pacific, and Southwest Pacific. The South China Sea population was then isolated from the Central-Southeast Pacific population due to the glacial low sea level and semi-closed sea basin during the early Pleistocene (ca. 1.3 Ma) (Wang, 1990; Woodruff, 2003).

Based on the close evolutionary relationship between the monophyletic East Indian Ocean clades II–III and Pacific Typical clade (IV) ( Figures 2 , S2 ), an counterclockwise South Equatorial Circulation could have carried larvae of the ancestral clade IV from the Southeast Pacific to the West Pacific and Indian Oceans through the open Palaeo-Indonesian Gateway ( Figures 5 , 6A ). The Indonesian Seaway and Indonesian Throughflow (ITF) were closely linked to the topography and bathymetry in this area (Supplementary Material). Closure of the Indonesian Gateway may have occurred more than once due to repeated tectonic events and sea level fluctuations since Miocene (Srinivasan and Sinha, 2000; von der Heydt and Dijkstra, 2006). Even if the deep-water exchange between the Pacific and Indian Oceans was restricted, surface circulation likely continued through narrower and shallower Indonesian seaway. The main composition of the ITF may have also gradually shifted from the warm South Equatorial Circulation to the cool North Equatorial water resulting from the closure of the Indonesian Seaway since 5–3 Ma (Cane and Molnar, 2001).

The divergence time of 1.1 (2.5–0.3) Ma between the East Indian Ocean clades (II–III) and Pacific Typical clade (IV) suggests an evolutionary influence of the early Pleistocene sea level fluctuation and closure of the Indonesian Seaway. Gene flow promoted by surface circulation between the Central-Southeast Pacific and East Indian Ocean could not have been interrupted completely until 1.1 Ma because a low proportion of the South Pacific water mass still kept flowing into the Indian Ocean through the Halmahera Sea and Banda Sea (Gordon and Fine, 1996). Nevertheless, frequent glacial sea level fluctuations over −120 m during the Pleistocene likely exposed shallow Indonesian seas and impeded gene flow between the Pacific and Indian Oceans (Woodruff, 2003). Similar speciation and lineage divergence due to ocean basins vicariance between the Indian Ocean and the Pacific Ocean were also observed in some other sister species pairs (Gaither and Rocha, 2013) and wide-ranging marine species including coconut crab (Birgus Zatro) (Lavery et al., 1996), green tiger prawn (Penaeus semisulcatus) (Abdul Halim et al., 2021), barnacle (Tetraclita kuroshioensis) (Chan et al., 2022), gastropod (Nerita albicilla) (Crandall et al., 2008), butterflyfishes species complex (Chaetodon punctatofasciatus and Chaetodon rhombochaetodon) (Mcmillan and Palumbi, 1995), and tropical reef fish (Chlorurus sordidus) (Bay et al., 2004). Thus, the closure of the Indonesian Seaway could have acted as a physical barrier to restrict transoceanic gene flow and promote population differentiation or speciation of widespread species between the Indian Ocean and the Pacific Ocean since the Miocene (Rocha et al., 2007).

Larvae of S. oualaniensis from the Central-Southeast Pacific firstly colonized the southern East Indian Ocean and then dispersed into the northern East Indian Ocean ( Figure 5 ). Monsoon circulations and Equatorial Countercurrent in Bengal Bay (Schott et al., 2009) subsequently separated the ancestral East Indian Ocean population into two lineages (II and III) ( Figure 6B ). A similar pattern has been inferred for the north–south population differentiation of D. gigas and S. oualaniensis, which are hypothesized to have been influenced by Equatorial Countercurrents in the East Pacific (Staaf et al., 2010) due to strengthened gyral circulations in the Indian and Pacific Oceans after the closure of the Indonesian Gateway (Kennett et al., 1985). Thus, oceanographic circulation systems could be responsible for population differentiation or speciation by impeding the gene flow of wide-ranging marine species (Tsang et al., 2012).