Episodic marine incursion has been a major driving force in the formation of present-day diversity. Marine incursion is considered to be one of the most productive ‘species pumps’ particularly because of its division and coalescence effects. Marine incursion events and their impacts on diversity are well documented from South America, North America and Africa; however, their history and impacts in continental East Asia largely remain unknown. Here, we propose a marine incursion scenario occurring in East Asia during the Miocene epoch, 10–17 Ma. Our molecular phylogenetic analysis of Platorchestia talitrids revealed that continental terrestrial populations (Platorchestia japonica) form a monophyletic group that is the sister group to the Northwest Pacific coastal species Platorchestia pacifica. The divergence time between the two species coincides with Middle Miocene high global sea levels. We suggest that the inland form arose as a consequence of a marine incursion event. This is the first solid case documenting the impact of marine incursion on extant biodiversity in continental East Asia. We believe that such incursion event has had major impacts on other organisms and has played an important role in the formation of biodiversity patterns in the region.
Marine incursions are the inundation of continental land by ocean waters, generally resulting from rises in sea levels and regional tectonic subsidence . The extensive lowland flooding produced by marine incursion has had a dramatic influence on landscape and ecosystem evolution . Historical marine incursion events have exerted profound effects on the dynamics of biodiversity, as best demonstrated in Neotropical communities [2–4]. Marine incursions produce progressively desalinized habitats and episodically allow aquatic connections between freshwater systems and the marine realm . These estuarine environments of varying salinity provide pathways for evolutionary habitat transition of many marine-derived organisms [6–8]. Marine incursions may have also functioned as impenetrable dispersal barriers for terrestrial species, fragmenting their previously continuous distribution ranges and leading to vicariant speciation . Because its division and coalescence effects often produce bursts of diversification, incursion is considered to be one of the most productive ‘species pumps’ . In addition, ecological impacts do not cease after termination of marine incursion, because soil properties and river water characteristics are permanently altered. This enhanced edaphic and ecosystem heterogeneity has sustained present-day high terrestrial diversity . There is little doubt that the combined effect of Andean uplift and multiple historical marine incursions have contributed significantly to the extraordinarily high biodiversity of the Neotropics [3,4,10]. Interestingly, despite increasingly available evidence of marine incursion for other continents, such as North America and Africa [11–13], few attempts have been made to document the impact of marine incursion on the East Asian continent .
Compared with the Neotropics, continental East Asia experienced similar, if not even more dramatic, topographic changes in the Neogene. Prior to the end of the Palaeocene, the eastern Asian foreland was higher than central and western Asia . The India–Asia collision occurring approximately 50 Ma led to the uplift of the Himalayan–Tibetan orogen . During the same time period, eastern Asian sedimentary basins experienced continuous subsidence . Approximately 23 Ma during the Neogene, uplift and subsidence collectively changed the landscape from an eastwest tilt to a westeast tilt . This reconfiguration dramatically altered drainage systems ; for example, the course of the Yangtze River, which initially flowed southwest to the Paratethys Sea, switched eastwards to the Pacific Ocean . These topographic changes also set the stage for marine incursions during the Early-to-Middle Miocene (23–11 Ma), a time of high global sea levels . Indeed, marine and brackish foraminifer, ostracod and gastropod fossils have been reported in eastern China [19,20], suggesting the occurrence of brief large-scale episodic marine incursions during the Miocene [7,21]. Surprisingly, very few investigations have been carried out to address the potential impact of such incursions on extant biodiversity; an effect on river dolphins is currently the sole documented case .
Testing historical marine incursion events and their biogeographic impact on extant biodiversity ideally requires a group of organisms that currently span both coastal and inland ranges. Amphipods of the genus Platorchestia (Amphipoda: Talitridae) are predominantly intertidal forms, with many intertidal species being common and widespread along the Northwest Pacific coast. Several terrestrial species, however, occur further inland on the East Asian continent. With distributions spanning both the coast and inland areas of East Asian continent, Platorchestia talitrids can thus serve as an excellent model system to test the marine incursion hypothesis.
We reconstructed a phylogeny of all terrestrial Platorchestia talitrids from continental East Asia and their potential coastal relatives using nuclear and mitochondrial DNA sequence data. We hypothesized that historical marine incursions are probably responsible for the origin of the inland terrestrial species. During an incursion episode, ancestral marine intertidal organisms invaded the inland; upon subsequent regression, they became isolated from their marine relatives and evolved a terrestrial lifestyle. We used a molecular phylogenetic approach to test this marine incursion hypothesis and we predicted that (i) terrestrial Platorchestia talitrids are nested within coastal species, with their most closely related coastal species occurring along the Northwest Pacific coast, the source of marine incursion and (ii) temporally consistent with incursion-related scenarios, divergence between terrestrial talitrids and their closest coastal relatives occurred during the Miocene, an epoch of high global sea levels and subsidence of the eastern Asian continental margin.
2. Material and methods
(a) Sample collection
Extensive coverage, particularly of coastal species, is essential for determining the closest coastal relative of terrestrial Platorchestia species. Samples of Platorchestia were collected (see the electronic supplementary material, table S1) from 49 sites covering most known distribution ranges along the Northwest Pacific coastline and continental East Asia. A total of 64 specimens of five species were examined, including one terrestrial species (P. japonica) and four coastal species (Platorchestia pacifica, Platorchestia pachypus, Platorchestia munmui and Platorchestia joi). Two species from genera Talorchestia and Sinorchestia were collected as outgroups. Detailed sample information is presented in the electronic supplementary material, table S1. All collected samples were immediately preserved in 95 per cent ethanol and stored at −20°C.
(b) Laboratory protocols
Total genomic DNA extractions and amplifications were performed following published protocols , with annealing temperatures optimized at 50°C–53°C. An 864-bp mitochondrial cytochrome oxidase I (COI) fragment was amplified using species-specific primers (forward primer 5′-TCTACTAATCACAAAGATATTGG-3′ and reverse primer 5′-CCAWCTAAAAACTTTRATACCKGT-3′). An approximately 410-bp portion of mitochondrial 16S was amplified with primers Y16F (5′-GGTAATTTGACCGTGCTAAG-3′) and Y16R (5′-CCGRTTTGAACTCARATCATGT-3′). A 327-bp histone H3 region was amplified with primers H3F and H3R , and an approximately 1340-bp nuclear 28S fragment was amplified with primer pair 28S-3311F and 28S-4434R . DNA sequencing was carried out using BigDye technology on an ABI 3730 automated sequencer (Applied Biosystems, Foster City, CA, USA).
(c) Construction of data matrices
DNA sequences were checked and edited with sequencher v. 4.1 (Gene Codes Corporation). Sequences of coding fragments COI and H3 were aligned using ClustalW  with default multiple alignment settings. Upon subsequent visual inspection, no insertions or deletions were found in either of these two regions. For sequences of length-variable 16S and 28S rDNA, we used two different alignment programs, ClustalW and mafft v. 5 , to assess sensitivity of our data to different alignment parameters . The preferred alignment for each locus gave the lowest score of incongruence with the protein-coding sequences. The number of variable sites and parsimony informative sites was calculated using MEGA5 .
(d) Phylogenetic analyses
Both maximum-likelihood and Bayesian inference were used for phylogenetic tree reconstruction, with analyses conducted on both individual and combined datasets. COI and 16S mitochondrial DNA (mtDNA) gene regions were concatenated as a single locus but partitioned in the phylogenetic analyses. Maximum-likelihood analyses were conducted using RAxML v. 7.0.3  starting with 1000 rapid bootstrap replicates followed by a thorough maximum-likelihood tree search. Each individual locus was analysed separately under the GTRCAT model of evolution; the concatenated dataset was analysed under the GTRMIX model, with each gene partition assigned a distinct model but with a joined branch length optimization. Bayesian phylogenetic analyses were conducted using MrBayes v. 3.1.2 . The best-fit models for each gene partition were selected based on the Akaike information criterion using MrModeltest v. 2.3 . For single locus, Bayesian analyses consisted of two runs of five million generations, each with four simultaneous Markov chains. Default priors were used and all states were considered equally probable. Tracer v. 1.4  was used to plot the resulting likelihood values and to determine when the Markov chains reached convergence. Trees were sampled every 500 generations; after discarding the first 3000 sampled trees as burn-in, we used the remaining trees to estimate a consensus tree and Bayesian posterior probabilities. For the combined data analysis, the data were separated into gene partitions with each partition following its best-fit evolutionary model. MrBayes was set to estimate parameters independently and simultaneously for each partition, and default priors were used. Each analysis consisted of two independent runs of 10 million generations, with sampling carried out every 500 generations. The first 5000 sampled trees were discarded as a conservative burn-in.
(e) Divergence time estimation
To specify interspecific levels of variation when estimating ages of divergence events, single- and multiple-threshold general mixed Yule coalescent (GMYC) analyses were optimized using the SPLITS package  in R. An ultrametric tree was constructed from all COI haplotypes of Platorchestia samples and outgroups under an uncorrelated lognormal relaxed molecular clock model using BEAST v. 1.6.2 , with the substitution model selected by MrModeltest and the substitution rate under the clock parameter setting set to one. A coalescent prior with constant population size was used to invoke the GMYC null model, with default options employed for all other parameters. We ran three independent analyses of 20 million generations each, sampling every 1000 generations; results were combined using LogCombiner v. 1.6.2 , with the first five million generations of each run discarded as burn-in. Effective sample size (ESS) values were used to determine the Bayesian statistical significance of each parameter (ESS > 200). Tracer was used to determine run convergence and stability.
We used two different strategies to estimate divergence times. First, divergence dates were estimated under a gene tree framework from the mtDNA data. A single exemplar representing each independent entity was selected for divergence dating analysis at interspecific levels. The sparse fossil record complicated the assessment of ages of certain lineages; the arthropod mitochondrial clock of 2.3 per cent divergence per million years  was therefore applied as calibration. The analysis implemented with BEAST was run using separate substitution models for COI and 16S partitions under an uncorrelated lognormal relaxed molecular clock model, and the Yule process was used to describe the speciation. Given the variation in topologies derived from the nuclear DNA data, we constrained the topology to that generated from maximum-likelihood analysis of the combined mitochondrial and nuclear DNA. Four independent simultaneous runs of 30 million generations were conducted with sampling carried out every 1000 generations. The resulting files from the independent runs, with the first 10 million generations excluded as burn-in, were combined using LogCombiner.
A more accurate estimate of divergence times is achievable via a fully Bayesian implementation of the multispecies coalescent carried out using a species tree framework [35,36]. This multispecies coalescent model co-estimates multiple gene trees embedded in a shared species tree, and incorporates uncertainty of the coalescent process. It has been suggested that this method has much better estimation accuracy for species tree topology than simple gene concatenation . Consequently, we also used a Bayesian Markov chain Monte Carlo method implemented in *BEAST v. 1.7.3 [36,37], to simultaneously estimate both a species tree and divergence times from our multilocus data. Because analysis in *BEAST requires a priori designation of taxonomic units, we defined the independent evolutionary entities identified by the GMYC model as distinct taxonomic units. We selected the best-fit substitution models calculated by MrModeltest, and unlinked substitution models, clock models and trees among loci. We applied Yule process for the species tree prior and specified a random starting tree for each locus. As described earlier, we used an uncorrelated relaxed clock model, with a mitochondrial rate of 2.3 per cent divergence per million years as the informative prior. We conducted four independent runs of 200 million generations each, sampling every 1000 generations and discarding the first 25 per cent as burn-in. Acceptable convergence to the stationary distribution was assessed by inspecting posterior samples using Tracer. Because all runs produced the same topology with very similar posterior probabilities, we combined the independent runs using LogCombiner and subsampled every third tree.
A total of 262 sequences of COI, 16S, 28S and H3 regions were successfully generated from 66 specimens in this study. Because 28S sequences obtained from two samples could not be interpreted unambiguously, only 64 sequences of 28S were used in subsequent analyses. Length of aligned protein-coding sequence datasets were 864 bp for COI and 327 bp for H3. Alignment parameters for length-variable fragments were selected to reduce incongruence with protein-coding sequences (see the electronic supplementary material, table S2). The final 16S and 28S alignments contained 416 and 1346 nucleotide sites, respectively. Detailed information on sequence variation is presented in the electronic supplementary material, table S3. All sequences were deposited in GenBank (for accession numbers, see the electronic supplementary material, table S1).
(b) Phylogenetic analyses
For the mtDNA data, best-fit models of evolution were determined to be GTR + I + G for COI and HKY + G for 16S. In the topologies recovered from both maximum-likelihood and Bayesian analyses (figure 1a), terrestrial Platorchestia talitrids appear as a well-supported monophyletic group (bootstrap value = 97 and Bayesian posterior probability = 100). This P. japonica clade is the sister species to coastal P. pacifica, although with low support. Both analyses strongly support each of the remaining species as monophyletic taxa, but the relationships among them are not strongly supported (figure 1a).
For the nuclear genes, best-fit models were GTR+I for 28S and HKY+I for H3. Overall, the two nuclear loci exhibited few informative characters, leading to unresolved gene trees in many cases, particularly with regard to interspecies relationships (see figure 1b,c and electronic supplementary material, table S3). There was no well-supported conflict between mitochondrial and nuclear gene trees. With respect to observed monophyletic taxa, both 28S and H3 trees are largely congruent with the mitochondrial phylogeny, except that P. japonica and P. pacifica are weakly supported as paraphyletic in the 28S tree (figure 1c).
The combined nuclear and mitochondrial gene dataset comprised 2953 bp. Maximum-likelihood analysis and Bayesian inference both support the same phylogenetic relationships between terrestrial and coastal species. All terrestrial lineages cluster as a well-supported (bootstrap value = 98 and Bayesian posterior probability = 100) monophyletic group, the P. japonica clade, which is nested within the coastal species (figure 2a).
The species tree estimated using *BEAST is identical to that generated from concatenated mitochondrial and nuclear data, again showing support for the monophyly of terrestrial lineages (Bayesian posterior probability = 100) and a sister relationship to coastal P. pacifica (Bayesian posterior probability = 96; figure 2b). For all nodes, support values obtained from species tree inference analysis are considerably higher than those derived from analysis of concatenated data (figure 2).
(c) Divergence times
To minimize the potential risks of mixing inter- and intraspecific levels of variation when estimating divergence times, we applied a GMYC model to specify independent entities at the interspecific level. Compared with the null model that assumes uniform branching rate for all lineages, both single- and multiple-threshold GMYC analyses provided significantly better fits to the ultrametric tree (likelihood ratio test, p < 0.05). Single- and multiple-threshold models delimited 22 (confidence interval = 18–24) and 24 (16–24) entities, respectively. Because the multiple-threshold model did not provide a significant improvement for the data (likelihood ratio test, p = 0.88), a single exemplar representing each independent entity delimited from the single-threshold model was selected for estimation of divergence times (figure 3).
Overall, our two dating approaches produced similar divergence time estimates. For most nodes, the 95% highest posterior density (HPD) of divergence events from the two analyses overlapped, but the means of estimates derived from the species tree analysis were slightly closer to the present (see the electronic supplementary material, figure S1). This is expected because divergence timing based on gene trees might overestimate than those derived using multilocus coalescent species tree-based approaches . Divergence time estimates derived from the species tree suggest that the terrestrial P. japonica diverged from its closest coastal relative at approximately 16.9 Ma (95% HPD: 13.6–20.3 Ma; node A in figure 4b and electronic supplementary material, figure S1), with initial divergence among the terrestrial forms occurring approximately 9.7 Ma (95% HPD: 7.5–12.1 Ma; node B in figure 4b and electronic supplementary material, figure S1).
The recovered phylogenies, molecular dating, and the biogeography of endemic Platorchestia from East Asia are consistent with our predictions of the marine incursion hypothesis. The generated phylogenies strongly support the monophyly of all terrestrial Platorchestia lineages and their derived relationship to coastal species (figure 2). This phylogenetic pattern suggests that a single land colonization event by marine intertidal ancestors gave rise to all current endemic terrestrial populations (P. japonica). Platorchestia pacifica, the sister species to P. japonica, has the largest geographical distribution of all coastal relatives from the Northwest Pacific. It is therefore reasonable to suggest that the most recent common ancestor of P. japonica and P. pacifica was a widespread species in the Northwest Pacific, which would have had the greatest opportunity to participate in a marine incursion event. According to molecular clock calibrations, terrestrial P. japonica diverged from its closest coastal species, P. pacifica, at approximately 16.9 Ma. This time coincides with a period of high global sea levels during the Early-to-Middle Miocene (figure 4b). The time of initial divergence among the terrestrial forms coincides with low global sea levels during the Late Miocene (figure 4b). In summary, these results are consistent with our temporal and geographical predictions under a marine incursion hypothesis.
Considering all the evidence, we reconstructed the following marine incursion scenario. During the Early-to-Middle Miocene, a time of high global sea levels, populations of the P. pacifica-like ancestor invaded inland areas of the East Asian continent, and became isolated from their coastal cousins. These populations evolved into the current P. japonica, while the coastal populations evolved into the current P. pacifica. Subsequent low sea levels in the Late Miocene led to fragmentation of the inland populations and vicariant divergence within P. japonica.
Several potential alternative explanations exist for the origin of terrestrial Platorchestia talitrids in East Asia. The occurrence of inland species could be the result of active and direct land penetration by ‘palustral’ ancestors . For instance, the land invasion may have been facilitated along coastlines by high rainfall blowing in from the sea. In such cases, invasions would be essentially opportunistic . Given the wide distribution of inland talitrids, however, we would expect them to have arisen on multiple independent occasions, which is not consistent with our results. Alternatively, surface streams might have indirectly facilitated passive dispersal of estuarial forms that then gave rise to the inland species . This explanation seems unlikely because talitrid amphipods exhibit direct development without a planktonic larval stage , and thus have limited long-distance-dispersal ability.
It has been recognized that marine incursion could play an important biogeographic role in organizing biodiversity . Incursion conditions produce progressively desalinized habitats that have been proposed as possible ‘competition troughs’ that can be exploited by newly arriving marine lineages . Consequently, this unique set of conditions may have facilitated ecological shifts of marine-derived taxa, such as molluscs, ostracods, fishes and dolphins [4,7,12,19]. Because marine incursion events also function as obvious dispersal and ecological barriers, particularly for resident terrestrial and freshwater organisms [9,40], congruent biogeographic patterns would reflect shared responses by different species to the same major vicariant event. Previous studies suggested that some resident organisms, including freshwater fishes, birds and frogs, may have been pushed into upland areas by elevated salinity in lowland regions [11,41,42]. Such situations would have isolated populations and over time led to allopatric speciation . Sediment composition, geochemistry and water characteristics were largely altered following the termination of marine incursions. This enhanced ecological heterogeneity provided an ideal setting for divergent natural selection in terrestrial and aquatic taxa, which in turn would generate and sustain high diversity in the region , such as in a recent study of croaker . Over the long term, the complex ecological and evolutionary trends initiated by Neogene tectonic events (e.g. topographic deformation and marine incursion) and continued by the Pleistocene climatic changes would be important for the generation and maintenance of biodiversity on the East Asian continent.
Species diversity in the terrestrial group may be much higher than currently recognized. Because of their morphological conservativity, all populations are currently circumscribed under the same name (P. japonica). Results from a recent study, however, suggest that terrestrial forms in Japan and Taiwan probably represent different species . The genetic diversity of the continental forms is high; for example, COI per cent divergence ranges from 2.0 to 14.1 per cent. The recent development of new species delimitation methods using multilocus data under a coalescent framework  enables a more efficient and less subjective approach to species delimitation. Application of these new methods and a more thorough examination will probably reveal multiple species within the continental terrestrial forms.
In conclusion, the phylogeographic patterns exhibited by Platorchestia talitrids provide the first solid evidence of a marine incursion event onto the East Asian continent leading to diversification of extant taxa. We expect that the shared signature of this incursion event is preserved in a diverse cross-section of East Asian continental organisms. Preliminary findings that point to an important role for marine incursions suggest that further biogeographic investigations of this palaeogeographic event will likely prove fruitful. Our understanding of the impact of historical events such as glaciation and marine incursion on this region of rich biodiversity is just beginning. Our study should open a new window to the evolution of ecosystems, as well as the origin and maintenance of endemic biodiversity patterns in East Asia.
We thank Jinzhong Fu, Associate Editor and two anonymous reviewers for constructive comments which greatly improved this manuscript. We thank Bryan C. Carstens, Anna Papadopoulou and Jean-François Flot for assistance with analyses. Ruiying Zhang helped with artwork. This study was supported by the National Natural Science Foundation of China (China National Funds for Distinguished Young Sciencists-31025023 to S.L. and NSFC-31172049 to Z.H.).
- Received December 4, 2012.
- Accepted February 7, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.