Evolution of a Planktonic Foraminifer during Environmental Changes in the Tropical Oceans

Ecological adaptation to environmental changes is a strong driver of evolution, enabling speciation of pelagic plankton in the open ocean without the presence of effective physical barriers to gene flow. The tropical ocean environment, which plays an important role in shaping marine biodiversity, has drastically and frequently changed since the Pliocene. Nevertheless, the evolutionary history of tropical pelagic plankton has been poorly understood, as phylogeographic investigations are still in the developing state and paleontological approaches are insufficient to obtain a sequential record from the deep-sea sediments. The planktonic foraminifer Pulleniatina obliquiloculata is widely distributed in the tropical area throughout the world’s oceans, and its phylogeography is well established. It is thus one of the best candidates to examine how past environmental changes may have shifted the spatial distribution and affected the diversification of tropical pelagic plankton. Such an examination requires the divergence history of the planktonic foraminifers, yet the gene marker (partial small subunit (SSU) rDNA) previously used for phylogeographic studies was not powerful enough to achieve a high accuracy in estimating the divergence times. The present study focuses on improving the precision of divergence time estimates for the splits between sibling species (genetic types) of planktonic foraminifers by increasing the number of genes as well as the number of nucleotide bases used for molecular clock estimates. We have amplified the entire coding regions of two ribosomal RNA genes (SSU rDNA and large subunit (LSU) rDNA) of three genetic types of P. obliquiloculata and two closely related species for the first time and applied them to the Bayesian relaxed clock method. The comparison of the credible intervals of the four datasets consisting either of sequences of the partial SSU rDNA, the complete SSU rDNA, LSU rDNA, or a combination of both genes (SSU+LSU) clearly demonstrated that the two-gene dataset improved the accuracy of divergence time estimates. The P. obliquiloculata lineage diverged twice, first at the end of the Pliocene (3.1 Ma) and again in the middle Pleistocene (1.4 Ma). Both timings coincided with the environmental changes, which indirectly involved geographic separation of populations. The habitat of P. obliquiloculata was expanded toward the higher latitudinal zones during the stable warm periods and subsequently placed on the steep environmental gradients following the global cooling. Different environmental conditions in the stable warm tropics and unstable higher latitudes may have triggered ecological divergence among the populations, leading to adaptive differentiation and eventually speciation. A comprehensive analysis of divergence time estimates combined with phylogeography enabled us to reveal the evolutionary history of the pelagic plankton and to find the potential paleoenvironmental events, which could have changed their biogeography and ecology.

Nevertheless, the evolutionary history of tropical pelagic plankton has been poorly understood, as phylogeographic investigations are still in the developing state and paleontological approaches are insufficient to obtain a sequential record from the deep-sea sediments. The planktonic foraminifer Pulleniatina obliquiloculata is widely distributed in the tropical area throughout the world's oceans, and its phylogeography is well established. It is thus one of the best candidates to examine how past environmental changes may have shifted the spatial distribution and affected the diversification of tropical pelagic plankton. Such an examination requires the divergence history of the planktonic foraminifers, yet the gene marker (partial small subunit (SSU) rDNA) previously used for phylogeographic studies was not powerful enough to achieve a high accuracy in estimating the divergence times. The present study focuses on improving the precision of divergence time estimates for the splits between sibling species (genetic types) of planktonic foraminifers by increasing the number of genes as well as the number of nucleotide bases used for molecular clock estimates. We have amplified the entire coding regions of two ribosomal RNA genes (SSU rDNA and large subunit (LSU) rDNA) of three genetic types of P. obliquiloculata and two closely related species for the first time and applied them to the Bayesian relaxed clock method. The comparison of the credible intervals of the four datasets consisting either of sequences of the partial SSU rDNA, the complete SSU rDNA, LSU rDNA, or a combination of both genes (SSU +LSU) clearly demonstrated that the two-gene dataset improved the accuracy of divergence time estimates. The P. obliquiloculata lineage diverged twice, first at the end of the Pliocene (3.1 Ma) and again in the middle Pleistocene (1.4 Ma). Both timings coincided with the environmental changes, which indirectly involved geographic separation of populations. The habitat of P. obliquiloculata was expanded toward the higher latitudinal zones during the stable warm periods and subsequently placed on the steep environmental gradients following the global cooling. Different environmental conditions in the stable warm tropics and

Introduction
Understanding the evolution of pelagic plankton during the drastic environmental changes of the past few million years is crucial for estimating the effects of global changes on future biodiversity. The pelagic plankton, particularly unicellular plankton, are known to disperse passively and widely because of the little motility during their entire life-cycles and the absence of geographic barriers in the open ocean (e.g. [1,2]). However, the phylogeographic studies have unexpectedly revealed the genetic differentiations and the limited dispersal among the geographically distant populations in the pelagic unicellular plankton [3,4,5]. Moreover, a recent study has demonstrated that speciation of pelagic protist occurred after differential adaptation of two sister species to trophic niches separated vertically in the water column [6]. Theoretically, pelagic plankton should be able to evolve rapidly in response to the environmental changes; however, their adaptation to a specific environment is still little understood [7]. Ecological differences between environments cause divergent selection that may lead to changes in various traits of individuals [8]. If these traits control reproductive isolation, adaptive divergence in those traits may reduce gene flow [8]. In the open ocean, therefore, the adaption to distinct environments can constrain gene flow without the split of populations by physical barriers. A combined analysis of divergence time estimate and phylogeography can help to find the tectonic and climatic events that likely may have changed the geographic distribution of pelagic plankton and their ecological characters. Thus, in the present study, we have attempted to determine whether any physical (geographic) separation of populations or rather ecological divergence related to the environmental changes is more effective in impeding the gene flow of pelagic plankton.
The tropical oceans are a biodiversity hotspot for both coastal and pelagic organisms [2,[9][10][11]. In particular, the tropical coasts of the Indo-Pacific Oceans have been a marine hotspot since the Miocene [10]. These warm environments expanded in latitudinal direction from the middle Miocene to Pliocene [12,13]. However, this warm state cooled around~4 million years ago (Ma) toward the Pleistocene owing to the geographic separations between the tropical Pacific and Atlantic Oceans [14] and the tropical Indian and Pacific Oceans [15], as well as to the atmospheric-oceanic environmental changes [16][17][18]. The modern geographic and climatic conditions have been established after the Pleistocene, when the large glacial/interglacial cycles started [16,17]. The dispersal and biogeographic patterns of the tropical coastal organisms changed, leading to speciation during the Plio-Pleistocene environmental developments, such as geographic connections and separations between the ocean basins [19,20]. However, little is known about the evolution of the tropical pelagic plankton that could be associated with the biogeographic and ecological changes in the past few million years.
Among pelagic plankton, the biodiversity of planktonic foraminifers has been well studied. They have been examined using the ribosomal gene marker [21] and are known to have one of the best fossil records since their appearance in the Jurassic [22]. The extant morphospecies, which are classified according to morphological differences of their calcareous shells, have appeared from the Miocene or Pliocene [23]. These morphospecies contain multiple genetically isolated species (i.e. genetic types; compiled in [21]), and almost all the existing genetic types are monophyletic within the morphospecies (e.g. [24,25]). Therefore, the divergences in these genetic types must have occurred after the appearance of the morphospecies lineage (after the Miocene). The morphospecies Pulleniatina obliquiloculata is distributed from the tropical to subtropical areas in the world's oceans since its first appearance at the end of the Miocene (~5.8 Ma) [23] and it is abundant in the tropical Indo-Pacific Oceans [26,27]. This species is a powerful candidate for examining the relationship between evolution and environmental changes in the tropical pelagic oceans. However, the sole use of fossils is insufficient in tracing the evolutionary history, as most of the tropical ocean basin is deeper than the carbonate compensation depth wherein the carbonates of foraminiferal shells are dissolved. Recently, the molecular phylogeographic studies demonstrated that the morphospecies P. obliquiloculata has three genetic types (Types I, IIa, and IIb) with distinctive distribution in the world's oceans [5,28] (Fig 1). Estimating the divergence times of these genetic types enables us now to understand the biodiversity responses to the environmental changes in the tropical oceans during the Plio-Pleistocene.
The molecular clock concept based on the substitution rates of nucleotide and amino acid sequences is useful for estimating the divergence time between the lineages. The divergence times of the genetic types of planktonic foraminifers were initially estimated using a strict molecular clock, in which the substitution rate was assumed to be the same among the lineages [31]. However, since the nucleotide substitution rates vary greatly among the planktonic foraminiferal lineages [32], this assumption introduced bias in the estimation of divergence time. Subsequently, the use of algorithms based on maximum likelihood (ML) and Bayesian frameworks improved the estimation of divergence for lineages with variable evolutionary rates (e.g. [33]). Recent studies have employed relaxed clocks to estimate the timing of genetic divergence in planktonic foraminifers [24,34,35], though the credible intervals (CI) in these estimations are still too large. Such large intervals (>3 million years) overlap with multiple geological Geographic distribution of the three genetic types of P. obliquiloculata. Pie chart shows the frequencies of three genetic types at each sampling site based on the previous studies [5,28]. The small and large circles indicate the locations where only one or more specimens were collected. The asterisks show the locations, where the samples were obtained for the present study. The genetic types of P. obliquiloculata specimens from Seears et al. [28] were identified based on the alignment of the partial SSU rDNA sequences with all existing data of Ujiié et al. [5]. Grey areas on the map show the ranges of water temperatures between 20-25, 25-28, and >28°C. The temperature data were obtained from the World Ocean Atlas 2005 [29]. The map is drawn by using the Ocean Data View [30]. events in the Plio-Pleistocene and obscure the relationships between the divergence times and specific environmental changes. All the previous molecular clock analyses used the partial SSU rDNA sequences. Although this~1000 base pair (bp) fragment (one-third of the entire length) contains useful variable regions (helices between 34 and 44) to identify the foraminiferal species and genetic types [21,36], it is difficult to achieve a precise estimation of the divergence time with such a short fragment. While the evolutionary rates vary among genes [37], the use of multiple genes can compensate this issue. Multiple gene markers like par actin-2, β-tubulin, and RPB1 are suitable markers for estimating the phylogeny at the family or higher levels of foraminifers [38,39]. Groussin et al. [40] used these genes in combination with the SSU rDNA to improve and succeed in estimating the divergence time with better accuracy for the early evolution of foraminiferal orders. Increasing the number of nucleotide bases and genes, which have enough resolution to show the phylogenetic relationship of the genetic types, facilitates more precise estimation of their divergence times in the past few million years. Thus, as the SSU rDNA and LSU rDNA have been used for phylogenetic analyses at the species level of foraminifers [41], we employed these genes to improve the precision of the divergence time estimation of the genetic types of P. obliquiloculata.
In the present study, we sequenced the entire coding regions of SSU rDNA and LSU rDNA from three morphospecies of planktonic foraminifers for the first time and investigated the phylogenetic relationships among the three genetic types of P. obliquiloculata. The divergence times estimated by using the Bayesian relaxed clock method based on different gene-datasets and fossil calibration sets were compared to examine the effects of increasing the number of nucleotide bases and different constraints on estimation. Then, we inferred a possible evolutionary scenario for P. obliquiloculata with respect to the environmental changes in the tropical ocean during the Plio-Pleistocene.

Materials and Methods
Nine P. obliquiloculata specimens derived from the samples that were previously collected at the three stations (St. C, D, and E) in the central Pacific equatorial area [5] were used in the present study (Fig 1). Four additional P. obliquiloculata specimens and two specimens each of Neogloboquadrina dutertrei and Globorotalia inflata, as outgroups, were newly collected from off Japan in the northwest Pacific. No specific permissions were required for the sampling location because it is the open sea. Our field studies did not involve endangered or protected species. All the examined specimens were individually cleaned in filtered seawater under a microscope.

DNA Extraction, Amplification, Cloning, and Sequencing
Genomic DNA was extracted from each specimen in accordance with the guanidinium thiocyanate protocol [5]. Nearly the entire coding region (~3500 bp) of SSU rDNA was amplified by polymerase chain reaction (PCR) analysis using the primers sA10 and sBf (Table 1). Two overlapping PCR-amplified fragments for LSU rDNA were assembled: the 5´region (~1250 bp) was amplified using primers L2f and L25fr, and the 3´region (~3500 bp) was amplified using primers L44 and LB ( Table 1) (Table 1, Fig 2). All the examined clones from a single specimen contained the identical rDNA sequence without single nucleotide polymorphisms. In total, 15 SSU rDNA sequences and 17 LSU rDNA sequences were deposited in the GenBank (accession numbers LC049304-LC049335).

Phylogenetic Analysis
Each nucleotide sequence of the SSU rDNA and LSU rDNA was manually aligned with the SeaView v4.3.4 program [42]. Published sequences of N. dutertrei and G. inflata (GenBank: EU199449 and EU199447) [43], which were amplified by the same primers sA10 and sBf, were also included in the SSU rDNA sequence dataset. After the ambiguously aligned sites had been excluded, 2685 sites of SSU rDNA sequences and 3591 sites of LSU rDNA were available for the phylogenetic analyses. Three analyses were performed: two with a single-gene (SSU or LSU) dataset and one with a two-gene (SSU+LSU) dataset.
For each single-gene analysis, the best-fit nucleotide substitution model was selected using MrModelTest 2.3 [44] and Treefinder [45] programs. For the SSU dataset, the general time reversible (GTR) [46] model was used with a gamma (Γ) [47] distribution for variable rates and a proportion of invariant sites (I) [48]. For the LSU dataset, a GTR model + Γ distribution were used. For the SSU+LSU dataset, we used separate model conditions that involved the individual optimization of all parameters for each gene.
Bayesian analyses were conducted on the three datasets with each of the optimal models, using MrBayes v 3.1.2 program [49]. The Markov Chain Monte Carlo (MCMC) process was set to enable the simultaneous functioning of four chains (three heated and one cold). Two independent runs were conducted for 1.6 × 10 6 (SSU), 1.1 × 10 6 (LSU) and 1.1 × 10 6 (SSU +LSU) generations. The trees and log-likelihood values were sampled at 100-generation intervals. The first 6 × 10 5 (SSU) and 1 × 10 5 (LSU and SSU+LSU) generations were excluded as burn-in. Pooled trees (1.0 × 10 6 generations) were used to obtain the Bayesian posterior probabilities for each dataset. The ML analyses for each of the three datasets were performed using Treefinder [45], and the bootstrap support was based on 1000 replicates in each dataset.

Bayesian Divergence Time Estimation
The divergence times of five lineages (three genetic types of P. obliquiloculata, N. dutertrei, and G. inflata) were estimated by the Bayesian relaxed clock method, as implemented in the MCMCTREE program (PAML 4.4b) [50]. A representative sequence was selected from each of three P. obliquiloculata clades (T94, KH279, and KH281). In the further analyses, we compared the CIs among the partial SSU rDNA between the helices 34-44 (the region between the primers s14F1 and sBf in Fig 2), two single-genes (SSU rDNA and LSU rDNA), and the two-gene datasets. The differences in the divergence estimates were also examined using four fossil calibration sets ( Table 2). The phylogenetic relationship of the fossils is inferred based on the morphological similarities of the shells (Fig 3A). Yet it is uncertain whether those morphological characters are inherited or not. Therefore, we considered four different combinations of calibration dates. The calibration set-1 used only the FAD of the extant species P. obliquiloculata for the minimum constraint of the node between P. obliquiloculata and N. dutertrei (Node 1) as the most relaxed constraint. The calibration set-2 employed the FAD of Pulleniatina praecursor, as the first divergence of the genus Pulleniatina, for the minimum constraint of Node 1. In accordance with the monophyletic origin of two lineages P. obliquiloculata and N. dutertrei as shown in the paleontological interpretation [23], the calibration set-3 added the FAD of common ancestor Neogloboquadrina acostaensis as the minimum constraint of the node between P. obliquiloculata + N. dutertrei and G. inflata (Root Node). Among the three extant species, P. obliquiloculata, N. dutertrei, and G. inflata, the lineage of G. inflata firstly diverged from the most common ancestor Neogloboquadrina continuosa (Fig 3A). Then, the tight constraint calibration set-4 used the FAD of N. continuosa as the minimum constraint of Root  Node and the FADs of N. acostaensis and P. praecursor as the maximum and minimum constraints of Node 1, respectively. In all datasets, the ML estimates of branch lengths were obtained under the GTR + Γ model with the gamma priors set at 0.5 using the BASEML program in the Phylogenetic Analysis by Maximum Likelihood (PAML) package [50]. The priors for the overall substitution rate were set at (a) G (1, 40.9) for the partial SSU dataset, (b) G (1, 30.9) for the SSU dataset, and (c) G (1, 34.3) for both of the LSU and SSU+LSU datasets. The prior for the rate-drift parameter was set at G (1, 2.32) for all the datasets. The independent-rates model [51] specified the prior rates among internal nodes. The parameters of the birth-death process for tree generation with species sampling [52] were fixed at λ = μ = 1 and ρ = 0. The root age and loose maximum bound for the root were set at 2.32 (23.2 Ma) and <3.72 (37.2 Ma) based on the FAD of Neogloboquadrina continuosa and Paragloborotalia pseudokugleri, which are the ancestral lineages of P. obliquiloculata + N. dutertrei and all three morphospecies studied including G. inflata [23] (Fig  3A). A hard and softbound version was used so that the probability of the true divergence time falls between the minimum and maximum bounds. The calibration nodes with minimum and maximum bounds assumed a heavy-tailed density based on a truncated Cauchy distribution of p = 0.1 and c = 1 as the default [50].
The MCMC approximation was obtained from a total of 1 × 10 4 samples, which were taken every 20 cycles after a burn-in period of 5 × 10 4 cycles. Two replicate MCMC runs with two different random seeds for each analysis were performed to confirm the convergence of the Markov chains to the stationary distribution. The MCMC samples from the two runs were combined after checking the distributions of the parameter values using the Tracer v 1.5 program [53]. The number of samples was large enough to reach the effective sample sizes (ESS > 200) for all parameters estimated in the present study.

Single-gene and two-gene phylogenies
Nearly complete SSU and LSU rDNA sequences were obtained from three morphospecies, P. obliquiloculata, N. dutertrei, and G. inflata. We confirmed that the LSU rDNA sequences in the present study branched from the published sequences of benthic foraminifers among the eukaryotes (data not shown). The phylogenies obtained from the two single-gene (SSU rDNA and LSU rDNA) datasets (S1 Fig) and the two-gene (SSU + LSU) dataset (Fig 4) showed the same topology, with three robust monophyletic clades associated with the three genetic types (Types I, IIa, and IIb) of P. obliquiloculata. The monophyletic clades based on the two-gene phylogeny showed the highest support levels (posterior probability of 1.00, bootstrap values of 99-100%).

Estimation of divergence times
The estimated ages and CIs were configured for assessing the prior probability at each node. The estimated ages were 27.0-24.7, 11.0-8.9, 4.2-3.1, and 2.1-1.2 Ma at Root Node, Node 1, Node 2, and Node 3, respectively (Table 3, Fig 3B). The differences between these ages at Root Node and Node 1 were relatively larger (~2.0 Ma) than those at Nodes 2 and 3 (~1.0 Ma). In the same fossil calibration set, the ages at each node ranged within 0.7 Ma among the four gene-datasets, except for the age at Root Node in the set-4 (Table 3).
In all estimations, the widths of CIs were found to be large in the partial SSU dataset, whereas they were found to be the smallest in the SSU+LSU dataset (Table 3, Fig 5). The CI widths tended to shorten gradually from the partial SSU, LSU, SSU to SSU+LSU datasets. Comparing the four fossil calibration sets, the largest differences (4.9 and 2.4 Ma) of the CI widths among genes at Nodes 2 and 3 were observed in set-2, whereas the smallest differences (3.2 and 2.0 Ma at Node 2 and 3) were found in set-4 ( Table 3). The shortest widths of CIs at Nodes 2 and 3 were 2.5 and 1.4 Ma in the SSU+LSU dataset in the calibration set-4.

Increasing phylogenetic resolution
All the phylogenies obtained in the present study clearly showed the presence of three genetically isolated species: the genetic types I, IIa, and IIb, in the morphospecies P. obliquiloculata (Fig 4, S1 Fig). In the previous study, the phylogeny based on the partial SSU rDNA sequences, corresponding with the helices 34-44 [36], showed that the Type IIb clade nested within a multidivergent clade of Type IIa [5]. In contrast, the phylogenies based on either the almost-complete SSU rDNA or LSU rDNA sequences showed the robust topology with distinct clades and appropriate phylogenetic relationships between the types IIb and IIa (S1 Fig). The combined analysis of the two genes increased the posterior probabilities and the bootstrap values (Fig 4). Therefore, we conclude that the increase in the number of nucleotide bases and genes is effective in elucidating the phylogenetic relationship between the species lineages.
The results of the divergence time estimations using four datasets showed that the width of CIs was shortened with the increase in nucleotide bases and genes (Fig 5, Table 3). In particular, the CI widths at Nodes 2 and 3 with the two-gene dataset were less than half of those with the partial SSU dataset, which has been used for the divergence time estimation in the previous studies [24,34,35]. The fossil calibrations also helped in achieving the precise estimation of the divergence times. This can be shown through the smallest width of CIs at Node 1 in the calibration set-4, in which both maximum and minimum time constraints restricted the prior probability distribution of Node 1 (Table 2, Fig 5). This strict calibration set affected the shortening of the CI widths at Nodes 2 and 3, as well (Fig 5). However, the CI widths have larger differences ranging from 3.2-4.9 and 2.0-2.4 Ma at Nodes 2 and 3 among the gene datasets than among the four calibration sets (1.5-3.0 and 0.6-1.3 Ma, respectively). This indicates that the precision of the divergence time estimation, in particular of the genetic types, is raised with increase in nucleotide bases and genes rather than by using the strict calibration set. The estimated divergence ages of Nodes 2 and 3 were found to be slightly different either among the four calibration sets (0.1-0.7 Ma differences) or the four gene datasets (0.6-0.7 Ma differences) ( Table 3). The estimated age varied slightly even if tight or relaxed calibrations were used. On the other hand, the ages estimated by the LSU rDNA dataset were relatively older than those estimated by the SSU rDNA dataset. These differences between the two single-gene datasets could be caused by the different evolutionary rates of genes [37]. However, all the estimated ages were within the CIs of SSU+LSU dataset in the most restricted calibration set-4 ( Table 3). The use of at least two different genes in addition to appropriate fossil calibration is crucial to obtain precise estimates of divergence times of planktonic foraminiferal genetic types. According to the SSU+LSU dataset with the calibration set-4, the lineage P. obliquiloculata diverged into Types I and II around 3.1 Ma during the Pliocene, and the lineage Type II split into Types IIa and IIb around 1.4 Ma in the Pleistocene.

Planktonic foraminiferal evolution and its relation to the ocean environmental changes
Recent phylogeographic studies have shown that the three genetic types (Types I, IIa, and IIb) of P. obliquiloculata are distributed from 24°13´S to 36°23´N of the world's oceans [5,28] (Fig  1). The genetic type I is widely found from the tropical to the temperate areas of the Atlantic, Indian, and Pacific Oceans, whereas Type IIa is found only in the Western Pacific Warm Pool (WPWP), which is the largest and warmest area in the world's oceans. The genetic type IIb is mainly found along the subtropical gyre margin of the Pacific Ocean, although one specimen was reported in the tropical Atlantic. Each of these genetic types is distributed in the mixed layer upper thermocline without any vertical separation [5]. All three types co-occur in the same water column in the WPWP. In the fossil record, the occurrence of the genus Pulleniatina has been found between 36.5°S and 41°N of the world's oceans since~7.0 Ma, as noted in the CHRONOS database (http://www.chronos.org). Since the appearance of P. obliquiloculata, its habitat has been restricted to the warm waters ranging from the tropical to temperate zones.
The estimation of divergence times from the molecular data suggest that the P. obliquiloculata lineage diverged twice, once at around 3.1 Ma during the late Pliocene and then at around 1.4 Ma during the middle Pleistocene (Fig 3B). The ocean environment was in the warm state from the Miocene to early Pliocene, in particular, the sea-surface temperature in the equatorial  areas was higher than the present state [12,13,18,[54][55][56]. In the Pacific, the WPWP was expanded in the latitudinal direction [54,55]. Moreover, high subsurface temperatures with a small east-to-west gradient occurred across the equator similar to the El-Niño event, although the depth of thermocline was entirely deeper than that of the modern El-Niño state [18]. The development of the warm ocean might have generated the wide and deep habitats for the warm marine species. Indeed, a high diversity of tropical coastal taxa (e.g. large benthic foraminifera, gastropods, fish, and corals) existed in the tropical oceans from the Miocene to Pliocene [10,57,58]. The warm state in the Pliocene possibly affected the expansion of habitat range for the pelagic species including P. obliquiloculata. However, the environmental conditions have shifted from the warm to cold states around~4 Ma, owing to the tectonic changes with the closure of the Panama Seaway [14] and the growth of the Indonesian Archipelago [15], and the atmospheric-oceanic changes [17,18]. This cooling involved a change in the water column structure by shoaling and cooling of the equatorial thermocline [18,56]. These late Pliocene changes steepened the environmental gradient towards the margins of the warm water area, which was the habitat of the warm water species. The tropical organisms such as P. obliquiloculata that once extended their habitats toward the higher latitude areas during the warm period, were confronted with different environmental conditions in their marginal habitat from the tropics in the late Pliocene.
After the Pliocene, the Walker circulation promoted the growth of Northern Hemisphere ice sheets [15]. The development of Northern Hemisphere glaciation was associated with an increase in the latitudinal temperature gradient, which could be one of the driving forces for strengthening the subtropical gyres [16]. Subsequently, large glacial/interglacial oscillations have occurred around 1.7 Ma [16,17]. Along with the subtropical gyre, the western boundary currents, such as the Kuroshio Current, transfer the ocean heat from the equator to higher latitudinal zones. However, the impact of the subtropical gyre and its surface current system weakened during the glacial period, and then the environment of the higher latitudinal zones (subtropical to temperate areas) tended to cool [59]. Following this latitudinal shift of the warm water, the tropical pelagic plankton P. obliquiloculata might have been able to migrate in the latitudinal direction during the interglacial period, whereas it was exposed to the steep environmental gradient during the glacial periods. Consequently, both divergences around 3.1 and 1.4 Ma probably coincided with steepening of the environmental gradient from the tropical to temperate oceans but not with any geographic separation between the ocean basins and/or areas, where P. obliquiloculata is distributed. Environmental diversification in habitat might cause ecological differences that could trigger divergent selection and adaptive differentiation in individuals leading to reproductive isolation [8,60]. The P. obliquiloculata population across the steep environmental gradient might have been ecologically differentiated between the core (tropical) and marginal (temperate) area. Such ecological changes may have forced the organisms to adapt to new environments if they were to survive. Although future studies are needed to investigate the ecological differences among the three genetic types of P. obliquiloculata, they have diverged when adaptive differentiation caused reproductive isolation between the core and marginal habitats.
The geographic distributions and divergence ages of P. obliquiloculata genetic types help in hypothesizing their evolutionary history, though it is impossible at the present moment to assign which oceans were their origins. The genetic type IIb is found in the tropical Atlantic and the Pacific, despite its divergence after the complete closure of the Panama Seaway between two oceans (Figs 1 and 3B). One possibility is that Type IIb could have migrated across the Pacific-Indian-Atlantic ocean basins through the global ocean circulation. However, no specimens of Type IIb were found in the Indian Ocean despite the thorough surveys [5,28] (Fig 1), which indicates a low possibility for inter-basin transfer through global ocean circulation.
Another plausible possibility could be that Type IIb is the descendant of the Type II lineage, which is globally distributed in the Pacific and Atlantic Oceans. The Panama Seaway was still open around 3.1 Ma, when the lineages of Types I and II diverged (Fig 3B). The eastward-flowing and westward-flowing currents facilitated the dispersal of marine taxa between the Atlantic and East Pacific Oceans until~2.6 Ma [57,61,62]. The populations of Type II might have remained in both oceans as the Type IIb lineage. Following this scenario, Type IIa diverged from the Type IIb lineage around 1.4 Ma. This young lineage have had narrow habitat range with temperature ranging between 25.4 and 29.2°C, compared with Type IIb having a wider temperature range between 13.9 and 29.2°C [5]. Taking into account the later divergence of Type IIa and its biogeography, Type IIa was considered as the residual population in the stable and warm areas of the tropical ocean. On the other hand, Type IIb could be the remaining survivor owing to the ecological differences generated among Type II populations along the sharp gradient in the warm waters.
The present study demonstrates that increasing the number of nucleotide bases of phylogenetically useful genes enables a precise estimation of the divergence times of planktonic foraminiferal genetic types. The genetic divergence of P. obliquiloculata has not caused by geographic (physical) separation but rather seems to be associated with ocean environmental changes in the warm water areas from the Pliocene to Pleistocene. An increase in precision of the divergence time estimation thus not only resolves the relationship between the evolutionary history and paleoenvironmental changes, enhancing our understanding of how the evolution of organisms is associated with environments.