Reconstructing the population history of the sandy beach amphipod Haustorioides japonicus using the calibration of demographic transition (CDT) approach.

Calibration of the molecular rate is one of the major challenges in marine population genetics. Although the use of an appropriate evolutionary rate is crucial in exploring population histories, calibration of the rate is always difficult because fossil records and geological events are rarely applicable for rate calibration. The acceleration of the evolutionary rate for recent coalescent events (or more simply, the time dependency of the molecular clock) is also a problem that can lead to overestimation of population parameters. Calibration of demographic transition (CDT) is a rate calibration technique that assumes a post-glacial demographic expansion, representing one of the most promising approaches for dealing with these potential problems in the rate calibration. Here, we demonstrate the importance of using an appropriate evolutionary rate, and the power of CDT, by using populations of the sandy beach amphipod Haustorioides japonicus along the Japanese coast of the northwestern Pacific Ocean. Analysis of mitochondrial sequences found that the most peripheral population in the Pacific coast of northeastern Honshu Island (Tohoku region) is genetically distinct from the other northwestern Pacific populations. By using the two-epoch demographic model and rate of temperature change, the evolutionary rate was modeled as a log-normal distribution with a median rate of 2.2%/My. The split-time of the Tohoku population was subsequently estimated to be during the previous interglacial period by using the rate distribution, which enables us to infer potential causes of the divergence between local populations along the continuous Pacific coast of Japan.

Introduction

Calibration of the molecular clock is one of the major challenges in marine population genetics. It is undoubtedly crucial to use an appropriate evolutionary rate to obtain demographic timelines as well as population parameters such as divergence time and migration rates [1, 2]. Although fossil records are most commonly used in phylogenetic rate calibration, they may not be applicable in population-level studies because little or no morphological difference among fossils is expected for younger episodes [3]. Geological dating is also unrealistic for some marine species because complete separation of populations is sometimes implausible because of gene flow among populations via passive transport (e.g., larval and egg dispersal) or active migration over long distances [4]. Phylogenetic evolutionary rates of related species are sometimes applied, but this is controversial because of the variation in molecular rates among organisms [5, 6].

Acceleration of the evolutionary rate is the other source of errors in molecular dating, and is particularly emphasized in population studies. Ho et al. [7] first summarized ideas about rate change by modeling the rate as an exponential function of time since calibration events and empirical studies subsequently verified this behavior in molecular clocks [8-10] (reviewed in [3]). This “time-dependency of the molecular clock” is partly explained by selection and genetic drift after deep phylogenetic events, which leads to reduction of the molecular rate [3, 7]. Grant [2] recommended caution in using a phylogenetic rate for population studies, especially for demographic reconstruction, because it can causes overestimation of the coalescence time, leading to misinterpretation of the results.

Calibration of demographic transition (CDT) [11] is one of the most promising techniques currently available for overcoming difficulties in rate calibration. CDT is a generalization of expansion dating [9] that assumes demographic expansion for some shallow-water invertebrate species in relation to an increase in habitat availability after the last glacial period. Both CDT and expansion dating use a two-epoch coalescent model (TEM), which assumes an ancient epoch of constant population size followed by a modern epoch of rapid population increase [12]. In expansion dating, the initial dates of the habitat increase (19.6 and 14.6 kya, based on two different assumptions) was used to obtain deterministic estimates of the evolutionary rate. CDT also assumes a post-glacial expansion but uses temperature as a demographic proxy in a calibration of species genealogy and obtains a stochastic function for the evolutionary rate. CDT is applicable for a wide range of species [11]. Regardless of its potential usability very few studies have used CDT, probably because its power in marine population genetics is currently not widely known.

We focused on the sandy beach amphipod Haustorioides japonicus Kamihira, 1977 (Amphipoda: Dogielinotidae) in the present study as a model for exploring population history because its biological traits potentially allow for reconstruction of the population history from molecular data. This species is characterized by an extremely low rate of effective migration and thus limited gene flow between local populations [13] because of the lack of a planktonic dispersal stage [14] and its specific habitat requirement for sandy beaches [15]. The species is distributed along the Sea of Japan and Pacific coasts of Hokkaido Island, but has also been recorded in the Tohoku region along the Pacific coast of northeastern Honshu Island, Japan (Fig 1) [16]. Considering the extremely low migration activity of the species, the populations in the Tohoku region are assumed to be genetically distinct from other known populations, with a local population history that is linked to the paleoceanography of the area.

Fig 1

Map showing sampling sites for the sandy beach amphipod Haustorioides japonicus.

Empty circles and numerals show locations and ID numbers for sampling sites in the present study. Filled and colored circles indicate the population assignment of the sites in a previous study (Takada et al. 2018 [13]; see also Fig 2). Shaded areas indicate the coastlines during the last glacial maximum (–120 m). Map insert shows previously reported sites with H. japonicus (stars, Kamihira 2000 [16]; circles, Takada et al. 2015 [15]). The map was created with QGIS v2.18.0 (http://www.qgis.org) using layers freely available at Natural Earth (https://www.naturalearthdata.com/downloads/10m-physical-vectors/).

Fig 2

Minimum spanning tree of mitochondrial COI haplotypes from the northwestern Pacific (NWP) population of the sandy beach amphipod Branch length and circle size are proportional to the T92 + Γ distance and haplotype frequency, respectively. The haplotype ID is beside each node (see also S1 Table). ECS, East China Sea; SJS, southern Sea of Japan; SJC, central Sea of Japan; NWP, northwestern Pacific; SJN, Northern Sea of Japan.

Here we incorporate CDT with a Bayesian Skyline Plot (BSP) [17] for the reconstruction of genealogy of the sandy beach amphipod Haustorioides japonicus, and demonstrate the importance of an appropriate molecular rate calibration in marine population studies. We exhaustively sampled the sandy beach amphipod in the Tohoku region to find the most peripheral populations of the species and assigned them to currently recognized sandy beach amphipod populations (Fig 1) [13]. We then inferred population histories by reconstructing demographic timelines, estimating the split-time of the local populations, and considering the biological traits of the species and the paleoceanography of the area.

Materials and methods

Field surveys and sample collection

A total of 32 specimens were collected from three sandy beaches in the Tohoku region and one site along the coast of Hokkaido (Fig 1, Table 1). Specimens were collected by using a 1-mm-mesh sieve, and parts or whole bodies of individuals were preserved in 6 M TNES (Tris HCl, EDTA, NaCl, SDS) urea buffer [18] (hereafter, “urea buffer”). The urea buffer contains a high concentration of urea, which allows cell lysis and DNA extraction and preservation at ambient temperature. Each specimen was identified to species in the laboratory according to Kamihira [14] using a stereomicroscopy.

Table 1

Sampling location, sample size (n), coordinates, collection date, and reference for sandy beach amphipod Haustorioides japonicus specimens used in this study.

ID	Location	n	Latitude	Longitude	Date	Reference
1	Oshamanbe	24	42.526	140.393	May 2015	Takada et al. [13]
2	Hakodate	20	41.770	140.743	May 2015	Takada et al. [13]
3	Tomari	3	43.031	140.524	May 2015	Takada et al. [13]
4	Usu	11	42.512	140.780	May 2015	This study
5	Ogawara	3	40.923	141.393	July 2018	This study
6	Misawa	3	40.674	141.439	July 2018	This study
7	Shichigahama	15	38.289	141.066	July 2018	This study

Laboratory procedures and data processing

Total DNA was extracted from the urea buffer-tissue lysate using a DNeasy Blood & Tissue Kit (Qiagen Inc., Valencia, CA, USA) according to the manufacturer’s instructions. Partial sequences of the mitochondrial cytochrome c oxidase subunit I (COI) gene were amplified by PCR using the primers LCO1490 and HCO2198 [19] (for details, see Takada et al. [13]). The products were directly purified by ExoSAP-IT (Affymetrix, Santa Clara, CA, USA) and sequenced in both directions at Food Assessment & Management Center (FASMAC) Inc. (Kanagawa, Japan).

The forward and reverse sequences were assembled by using Mesquite v3.01 [20] and visually inspected (accession numbers: LC474498–LC474506, also see S1 Table). The 105 sequences used by Takada et al. [13] (accession numbers LC224174–LC224278) were added to the dataset and multiple sequence alignment was subsequently performed using MAFFT v7.294 [21]. The best nucleotide substitution model for the dataset was determined using MEGA7 [22]. The sequence dataset was then collapsed into haplotype using custom Perl script.

To assign sequences from the Tohoku region to currently known sandy beach amphipod populations, we constructed a minimum spanning tree using the software packages Arlequin v3.5 [23] and SplitTree v4 [24]. Genetic differentiation between the populations was confirmed by the pairwise Φst test. We also obtained standard genetic diversity indices (H, number of haplotypes; h, haplotype diversity; π, nucleotide diversity) and neutrality indices (Tajima’s D [25]; Fu’s Fs [26]; R2 [27]). These statistical tests were performed with Arlequin except for R2 test done on DnaSP v6 [28]; significance of the pairwise Φst and neutrality indices was assessed with 1000 permutations.

Rate calibration using the CDT approach

Rate calibration was performed following standard CDT procedures [11], but with some minor modifications (hereafter ‘rate’ means substitution rate, not ‘divergence’ rate). To minimize coalescent errors and improve convergence, we built TEMs using all individuals from Tohoku as well as its parental population among currently known sandy beach amphipod populations. We created custom .xml files for BEAST v1.7 [29] for the combination of two clock models (strict clock, log-normal and exponential relaxed) and two TEMs (exponential and logistic; also see S1 Fig), referring to the example file of Hoareau [11] (Dryad Digital Repository; http://dx.doi.org/10.5061/dryad.3q24t). A total of 10,000 out of 5.0 × 108 Markov chain Monte Carlo (MCMC) steps were recorded by using BEAST. The results from the six models were then compared by using the Bayes factor test implemented in Tracer v1.6 [30]. After discarding the first 1000 records (i.e., 10%) as a burn-in, 9000 values of the mutation-scaled transition time from the best model were obtained. We subsequently downloaded 10,000 calendar years derived from the discretized change rate of temperature over the last 25 ky (Dryad Digital Repository; see Jouzel et al. [31] for original data) to obtain 9 × 107 mutation-rate values. The density distribution of the mutation rate was then modeled using log-normal and Gaussian functions, and the best model was selected based on the Akaike information criterion (AIC). These statistical calculations were performed using MASS v7.3 package [32] in R v3.4.3 [33].

Reconstruction of population history

We inferred the demographic history of the Tohoku and its parental population using a BSP [17], and the genealogy derived from the Bayesian skyline model. We first built a BSP with a total of 10,000 records out of 2.0 × 109 MCMC steps and discarded the first 1000 records (i.e., 10%) as a burn-in. The number of groups was set to 10 and a piecewise-linear BSP model was adopted. The molecular rate was modeled following the best clock model and statistical distribution from CDT. We also built a BSP under a strict clock model with a conventional evolutionary rate for crustacean mitochondrial genes of 0.7%/My [34] as a proxy of ‘general’ rate to compare results from different molecular clock assumptions. The site model was set to the best model inferred by MEGA. A generation time of one year [14] was assumed to convert generation time-scaled female effective population size into effective population size (Ne). The split-time of the Tohoku population from the parental population was subsequently estimated based on the genealogy of the mitochondrial COI and reconstructed along with the demographic changes in the BSP. TreeAnnotator v1.7 [28] was used for time-tree reconstruction. We used 9000 trees after discarding the first 1000 trees out of 10,000 and obtained a maximum clade credibility tree with a median node height. The posterior distribution of the split-time between the Tohoku and northwestern Pacific (NWP) populations was then summarized using TreeStat v1.7 [28].

Results

All of the 32 newly-obtained sequences were assigned to a group that is common along the coast of the northwestern Pacific Ocean (Fig 2), and the following analyses were therefore performed on 79 sequences for this genetic group, including the 47 previously reported (Table 1). The final dataset comprises 608 bp of COI sequences from 79 individuals, which were collapsed into 15 haplotypes defined by 17 polymorphic sites. The best-fit model was estimated to be HKY + Γ [35, 36]. Frequencies and accession numbers for the 15 haplotypes are given in S1 Table.

Among 79 individuals of the NWP group, 15 from the southern Tohoku region (site 7, Fig 1) are genetically distinct from the other members of the NWP group (Φst = 0.75, P < 0.05; Fig 2). We therefore treat the group of 15 individuals from site 7 as a distinct population (hereafter “Tohoku population”) and we explore its history in the following section. Table 2 includes summary statistics for the two group of populations (NWP and Tohoku). Neutrality indices for these populations were negative and significant except for the R2 value of Tohoku population.

Table 2

Genetic diversity and neutrality indices for two group of populations of the sandy beach amphipod Haustorioides japonicus.

Population	Standard diversity indices				Neutrality indices
	N	H	h	π	D	Fs	R2
NWP	64	11	0.62 ± 0.05	0.0014 ± 0.0011	-1.92*	-7.48*	0.042*
Tohoku	15	4	0.37 ± 0.15	0.0009 ± 0.0009	-1.82*	-1.72*	0.143

n, number of individuals; H, number of haplotypes; h, haplotype diversity; π, nucleotide diversity; D, Tajima’s D; Fs, Fu’s Fs. Asterisks (*) indicate significance of tests at P < 0.05.

A logistic TEM with exponential relaxed clock (model 3, Table 3) was chosen as the best model in the Bayes factor test, although there was no substantial difference between median estimates of transition time in the four models. For the exponential TEM with strict clock, we could not obtain the result because the runs did not converge. Hereafter we used the best TEM only. The density distribution of transition time was unimodal and the median estimate was 2.58 × 10−4 mutations/site in the TEM (Fig 3). We selected the log-normal model, which had a lower AIC than the Gaussian model for the density distribution of evolutionary rate (S1 Fig); LogMean and LogSD for the log-normal approximation were –3.89 and 0.53, respectively. The median CDT rate was 2.2%/My.

Table 3

Model comparison results.

Listed are the model ID, type of two-epoch model (TEM) and clock model, transition time (TT, 10−3 mutations/site), median marginal likelihood (MML), and 2×log Bayes factors (2×lnBF) for the models. Parameters for the best-fit model are shown in bold type. The result for the exponential TEM with strict clock is not shown because the run did not converged.

Model	TEM	Clock model	TT	MML	2×lnBF
					#1	#2	#4	#5	#6
1	Exponential	Exponential	2.60	-1012.3		7.2	-6.8	-4.8	8.0
2	Exponential	Lognormal	2.51	-1015.9	-7.2		-12.0	-10.0	-8.0
3	Logistic	Exponential	2.56	-1008.9	6.8	14		2.0	14.8
4	Logistic	Lognormal	2.53	-1009.9	4.8	12.0	-2.0		12.8
5	Logistic	Strict	2.59	-1016.3	-8.0	-0.8	-14.8	-12.8

Fig 3

Demographic timeline derived from the best two-epoch demographic model (model 3, Table 3).

An ancient epoch with a constant population size (θ1) is followed by the modern epoch of logistic population growth (θ2). The median estimate of the transition time is shown as a vertical red line. Horizontal and vertical blue lines and the blue shaded area show the median estimate for mutation-scaled female effective population size (θ), the time of the most recent common ancestor, and the 95% credible interval, respectively. The histogram shows the posterior density distribution for the two-epoch transition time.

The BSP was built with exponential relaxed clock using all sequences from NWP and Tohoku populations and the CDT rate, revealed rapid growth after the last glacial period, coinciding with warming temperatures, whereas the time of demographic expansion was estimated to be during the last glacial maximum (20–40 kya) with a phylogenetic rate of 0.7% (Fig 4). The divergence between NWP and Tohoku populations falls back to the previous interglacial period (median = 97 kya, but with a 95% credibility interval of 11–360 kya; Fig 5) using the CDT rate, whereas the estimate is 340 kya (95% credibility interval: 61–450 kya) based on the phylogenetic rate.

Fig 4

Bayesian skyline plots based on the calibration of demographic transition (CDT) rate and a conventional mitochondrial evolutionary rate of 0.7%.

Median estimates of female effective population size are shown as lines and vertical lines show median estimates of the time of the most recent common ancestor. Global trends of historical temperature are superimposed on the plots (Jouzel et al. [31]).

Fig 5

Posterior density distribution histograms of the split-time between Tohoku and other NWP individuals estimated by using CDT rate (A) and a conventional mitochondrial evolutionary rate of 0.7% (B), and global trends of historical temperature from Jouzel et al. [ Maximum clade credibility trees are also shown for each rate.

Discussion

It is generally difficult to explore the histories of marine populations because vicariant events and fossil records are rarely applicable to the rate calibration. Even if these calibration techniques happen to be available for a target species, they can cause misinterpretation of the results because of the time dependency of the molecular clocks [2]. The use of demographic transition as a global rate calibration point enables us to overcome these difficulties in marine population genetics and attribute recent coalescent events such as population splits to certain geological episodes. In the following section, we elucidate the history of the sandy beach amphipod population and compare our results from the CDT approach with those based on a phylogenetic rate. We then discuss the limitations of the current approach.

Fluctuations of temperature as a proxy for demography

We observed a strong concordance between fluctuations in temperature and historical demography of the sandy beach amphipod population (Fig 4). Post-glacial expansion in the NWP is much more plausible than population growth throughout the last glacial maximum, considering the life-history traits of the species and the paleoceanography of the northwestern Pacific Ocean. The sandy beach amphipod population grows rapidly in late spring to summer from a small number of overwintering individuals [37]. The extended winter season during the last glacial period probably shortened the reproductive season for the species, which could have caused a serious decrease in effective population size. Sea-ice cover in the area expanded both spatially and temporally during the last glacial period in this region [38], which could have also altered the environmental conditions on sandy beaches. After the last glacial period, the population seems to have rapidly recovered from the severe reduction in size in parallel with warming temperatures (Fig 4). Our results from the CDT rate seem quite reasonable considering that the North Atlantic intertidal gastropod Littorina saxatilis Olivi, 1792 whose habitat conditions are close to those of the sandy beach amphipod, also showed post-glacial expansion [11].

Divergence between the local populations

In the NWP, sea-level changes throughout glacial cycles have greatly influenced marine populations via changes in both inter-population connectivity and oceanographic conditions [39-43]. Because the Tsugaru Strait connecting the NWP and the Sea of Japan was narrow and shallow during the last glacial period, the Tsushima Current from the Sea of Japan (Fig 1) ceased [44], possibly causing lineage splits for marine species. Considering that the posterior distribution of the split-time shows a peak during the last glacial period (Fig 5), it can be assumed that gene-flow among the sandy beach amphipod populations was lost during the last glacial period as a result of the change in oceanic currents, resulting in genetic differentiation of the Tohoku population (Fig 1).

Alternatively, we also hypothesized that occasional migration events contributed to the southward range expansion of the species. We note that tsunamis are one of the possible transport mechanisms in the Tohoku region because they have repeatedly swept sandy beaches in this area [45]. Most recently, a megathrust M9.0 earthquake generated a huge tsunami on 11 March 2011, which caused serious damage in shallow-water environments and greatly altered benthic community structure in this region [46-49]. We discovered, however, that the local sandy beach amphipod population in the Tohoku region persisted after the tsunami, which strongly suggests the survival or resettlement of the local population. It is thus plausible that there were occasional migrations between local sandy beaches driven by episodes such the Great Tsunami. Unfortunately, it is difficult to attribute a population split to a specific geological event, owing to the large credible interval (Fig 5). Increasing the number of individuals, loci, and sites sampled would improve the precision of the split-time estimates and is thus a promising avenue for further exploration of the population history of this species.

Limitations of the demographic rate calibration for the sandy beach amphipod

Rapid growth of a population after the last glacial period is a fundamental assumption of CDT [11]. Although we have discussed the plausibility of the post-glacial expansion hypothesis for this species in the previous section, there still remains some uncertainty. The CDT rate from the present study exceeds the phylogenetic rate of 1%, but is much slower than that from Hoareau [11], possibly signaling an underestimation of the rate. We also note that there is a range among observed rate for the crustacean species [9, 50, 51]. The difficulties are mainly caused by the fact that ultimately it is impossible to know the real trigger for the population expansion. A false assumption of post-glacial expansion is a potential limitation of the demography-based rate calibration technique [9, 11], and this possibility must always be considered.

Additional sources of uncertainty in BSP and estimation of the divergence time are large coalescent errors caused by a small number of sequences, loci, and polymorphic sites [2]. We used only a single locus of COI with 608 bp from 79 individuals, which could lead to misinterpretation of the results. Also, using the probabilistic CDT approach increases the range of the credible interval for the BSP and split-time, which makes it difficult to interpret the results, as discussed in the previous section. The two-epoch demographic model is currently applicable only for a single gene in BEAST, and this becomes a more practical problem when one wishes to improve the credibility of the analysis by using a larger dataset. The implementation of multi-locus data in the program will help to deal with this issue.

Conclusion

We have demonstrated the applicability of rate calibration with global demographic transition by using the CDT approach. The median estimate of the CDT rate was approximately 2%/My for COI sequences from the sandy beach amphipod, which exceeds the conventional phylogenetic rate of 0.7%/My. The BSP for the entire population was more realistic using the CDT rate than the phylogenetic rate, which detected population growth throughout the last glacial maximum, in terms of the post-glacial expansion (Fig 4). We further inferred the drivers of the population split by estimating a divergence time from the CDT rate. Expanding the dataset will enable further exploration of the population history of marine organisms such as the sandy beach amphipod.

Histogram showing the frequency distribution of evolutionary rate derived from CDT.

Log-normal and Gaussian approximations of the rate are also shown.

(PDF)

Click here for additional data file.

List of haplotypes.

Haplotype ID, total number of individuals (n), number of individuals per site (Per site), and accession number are shown.

(PDF)

Click here for additional data file.

31 Jul 2019

PONE-D-19-16483

Reconstructing the population history of the sandy beach amphipod Haustorioides japonicus using the calibration of demographic transition (CDT) approach

PLOS ONE

Dear Dr. Sakuma,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

We would appreciate receiving your revised manuscript by Sep 14 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Tzen-Yuh Chiang

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating in your Funding Statement:

[This study was partly supported by “A project for development of assessment methods for the coastal environment” funded by the Ministry of Agriculture, Forestry and Fisheries, and Grants-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan (25340114 and 15H02265).].

* Please provide an amended statement that declares *all* the funding or sources of support (whether external or internal to your organization) received during this study, as detailed online in our guide for authors at http://journals.plos.org/plosone/s/submit-now.  Please also include the statement “There was no additional external funding received for this study.” in your updated Funding Statement.

* Please include your amended Funding Statement within your cover letter. We will change the online submission form on your behalf.

3. We note that [Figure 1] in your submission contains a map image which may be copyrighted. All PLOS content is published under the Creative Commons Attribution License (CC BY 4.0), which means that the manuscript, images, and Supporting Information files will be freely available online, and any third party is permitted to access, download, copy, distribute, and use these materials in any way, even commercially, with proper attribution. For these reasons, we cannot publish previously copyrighted maps or satellite images created using proprietary data, such as Google software (Google Maps, Street View, and Earth). For more information, see our copyright guidelines: http://journals.plos.org/plosone/s/licenses-and-copyright.

We require you to either (1) present written permission from the copyright holder to publish these figures specifically under the CC BY 4.0 license, or (2) remove the figures from your submission:

You may seek permission from the original copyright holder of Figure(s) [1] to publish the content specifically under the CC BY 4.0 license.

We recommend that you contact the original copyright holder with the Content Permission Form (http://journals.plos.org/plosone/s/file?id=7c09/content-permission-form.pdf) and the following text:

“I request permission for the open-access journal PLOS ONE to publish XXX under the Creative Commons Attribution License (CCAL) CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). Please be aware that this license allows unrestricted use and distribution, even commercially, by third parties. Please reply and provide explicit written permission to publish XXX under a CC BY license and complete the attached form.”

Please upload the completed Content Permission Form or other proof of granted permissions as an "Other" file with your submission.

In the figure caption of the copyrighted figure, please include the following text: “Reprinted from [ref] under a CC BY license, with permission from [name of publisher], original copyright [original copyright year].”

If you are unable to obtain permission from the original copyright holder to publish these figures under the CC BY 4.0 license or if the copyright holder’s requirements are incompatible with the CC BY 4.0 license, please either i) remove the figure or ii) supply a replacement figure that complies with the CC BY 4.0 license. Please check copyright information on all replacement figures and update the figure caption with source information. If applicable, please specify in the figure caption text when a figure is similar but not identical to the original image and is therefore for illustrative purposes only.

The following resources for replacing copyrighted map figures may be helpful:

USGS National Map Viewer (public domain): http://viewer.nationalmap.gov/viewer/

The Gateway to Astronaut Photography of Earth (public domain): http://eol.jsc.nasa.gov/sseop/clickmap/

Maps at the CIA (public domain): https://www.cia.gov/library/publications/the-world-factbook/index.html and https://www.cia.gov/library/publications/cia-maps-publications/index.html

NASA Earth Observatory (public domain): http://earthobservatory.nasa.gov/

Landsat: http://landsat.visibleearth.nasa.gov/

USGS EROS (Earth Resources Observatory and Science (EROS) Center) (public domain): http://eros.usgs.gov/#

Natural Earth (public domain): http://www.naturalearthdata.com/

4. We noted in your submission details that a portion of your manuscript may have been presented or published elsewhere. [Some nucleotide sequences (see Table S1) used in the present study were published previously and available online. The aim of the present study and thus analytical procedures are different from that of the previous one and also, we used additional sequences. This submission therefore does not constitute dual publication.]

* Please clarify whether this  publication was peer-reviewed and formally published. If this work was previously peer-reviewed and published, in the cover letter please provide the reason that this work does not constitute dual publication and should be included in the current manuscript.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: Sakuma and colleagues have reconstructed the demographic history of several populations of the amphipod Haustorioides japonicus using an approach based on demographic transition from an ancient epoch of constant population size to a modern epoch of rapid increase. The rapid population increase is assumed to have been triggered by warming temperatures following the Last Glacial Maximum. This approach is very interesting and has rarely been explored. As such, the present study is important as it explores the strengths and weaknesses of the method. The manuscript is well written, concise, the language is very good and the figures are of good quality. Overall, the analyses are done properly, although some clarifications are needed here and there. I have some remarks that should improve the presentation of the paper.

Line 41 “change of rate temperature” should be “rate of temperature change”.

Line 147. I recommend the authors to also use the R2 test (implemented in DnaSP) since it is more sensitive, especially regarding small sample sizes (Ramos-Onsins, S.E., Rozas, J., 2002. Statistical properties of new neutrality tests against population growth. Mol. Biol. Evol. 19, 2092–2100).

Line 154. “we built TEMs based on the whole sequence of the population to which the Tohoku population was assigned”. This needs to be rephrased as it is not clear. Does it mean that TEMs were build using all individuals from the NWP populations?

Line 156. Why did the authors not consider a strict clock as well for the BF comparisons?

Line 163. “9000 mutation values from the best model were obtained”. Please specify exactly to what parameter this refers to from the BEAST output file.

Lines 165-169. Please give more details here. What packages and their versions were used for these analyses?

Line 171. “inferred the history” should be “inferred the demographic history”. Write what BSP means when mentioned the first time. It is a bit confusing because the authors mention that they inferred the history of the Tohoku population, when in fact they analyzed this population together with the other NWP populations. This is made clear only in the results, but it should be clear here as well. Also in the BSP analyses it is not clear what kind of clock was used (strict vs relaxed). The authors mention the strict clock only in the second model based on the 1% evolutionary rate.

The authors compare their obtained rate with a phylogenetic (interspecific) rate of 1%. Although I see the use of such comparison, I do not see the point of the arbitrarily chosen 1% rate. I would recommend the authors to compare their intraspecific rate with interspecific (phylogenetic) rates that have been inferred for other crustaceans (e.g. Knowlton & Weigt 1998 https://doi.org/10.1098/rspb.1998.0568) and recently even amphipods (Copilas-Ciocianu et al. 2019 https://doi.org/10.1007/s13127-019-00401-7).

The authors should also clarify what they mean by “rate”. There is substitution rate and divergence rate (which is 2x substitution rate). In their MS, I suppose that the authors refer to the substitution rate. There is confusion in the literature about this, and that is why there should be more clarity (Schenekar & Weiss 2011 doi:10.1038/hdy.2011.48).

Line 325. Here the authors should also mention that several studies have shown fast intraspecific rates in crustaceans as well (Audzijonyte&Vainola 2006 http://dx.doi.org/10.1111/j.1365-294X.2006; Crandall et al. 2012 http://dx.doi.org/10.1093/molbev/msr227).

Reviewer #2: I found your work on the sandy beach amphipod Haustorioides japonicus really interesting. It represents an important and interesting contribution as it shows, for the first time, the power of the calibration of demographic transition approach and the importance of using an adequate evolutionary rate in an amphipod species.

In general, the organization of the manuscript is satisfactory and its easy reading. The TITLE clearly reflects the contents. The ABSTRACT is sufficiently informative and has a correct length. The INTRODUCTION is very interesting, clear and concise, and contains enough background to put the reader in context. The statement of the objective is adequate and appropriate. MATERIAL AND METHODS are clearly explained, being sufficiently informative to allow replication. The analyses are properly detailed, and they are solid. The RESULTS are clearly presented and described in a logical order, and their interpretation and posterior discussion are justified by the data and consistent with the objectives. The DISCUSSION and CONCLUSIONS are very clear and well organized. Authors explain properly their results obtained with many adequate references. Finally, the English is adequate.

Overall, I think the paper is interesting and worthy of publication. I think that this paper could be accepted in PLOS ONE as the subject of the manuscript falls within the scope of the journal. Some suggestions and corrections were made to improve it (see below) that should be relatively easy for the authors to fix. So, I consider that this manuscript is acceptable for publication after minor revisions.

IMPORTANT:

According to the “World Register of Marine Species” (WoRMS) database, Haustorioides japonicus is an unaccepted species. The accepted species name corresponds to “Eohaustorioides japonicus (Kamihira, 1977)”. Therefore, the authors must correct this along the whole manuscript including the corresponding title and the abstract.

INTRODUCTION

- Line 86: “(Dogielinotidae: Amphioda)”. The Order name must be before the Family name. Therefore, the authors should place “Amphipoda” before “Dogielinotidae”.

- Line 95: “genetically distinct”. Genetically distinct from what? The authors should clarify this to avoid confusion.

- Line 97: “BSP”. What does BSP mean? Bayesian Skyline Plot? When you use an abbreviation, its meaning should be specified the first time it appears in the manuscript. Therefore, the authors should include this information.

- Line 97: “genealogy”. Please, include “of the sandy amphipod Eohaustorioides japonicus”, just to make it clear.

MATERIALS AND METHODS

- Line 109 - 110: “three sandy beaches in the Tohoku region and two sites along the coast of Hokkaido”. In the present study, according to Figure 1 legend (empty circles) and locations in Table 1, only one location was sampled along the coast of Hokkaido (site 4, Usu). Therefore, the authors should replace “two sites along…” with “one site along…”

- Line 111: “parts”. Which parts were preserved? Did the authors preserve some specific parts for some reason? Please, specify this.

- Line 115: Please, include a “a” after “using”.

- Line 135: I will include a sentence saying that sequences were deposited in GenBank (accession numbers: LC474498–LC474506) (S1 Table).

- Line 137: “The 300 sequences”. Sequences used by Takada et al. [13] were 105 and not 300, as the authors can see in the Table S1 of Takada et al.’s paper. The authors should correct this.

RESULTS

- Line 207: “(site 7)”. The authors should include “Fig 1” after “(site 7)” to ensure that readers understand exactly which location is talking about.

- Line 211: “two populations”. Really, NWP is not only one population but a group of them. Therefore, I think, the authors should replace “two populations” with “two group of populations”. The same should be modified in Table 2 captions.

DISCUSSION

- Line 289: “Littorina saxatilis”. The authors should include the authority of this taxon.

REFERENCES:

- Line 412: Please, remove the comma after “Kamihira”

- Line 418: Please, remove the comma after “Kamihira”

- Line 465: “Haustorioides japonicus” must be in italics. Please, change it.

FIGURES AND TABLES:

- Figure 1: It is clear and appropriate for the data being presented. However, I think, it is necessary to include/specify two things:

o According to the manuscript, Tohoku region is located along the Pacific coast of northeastern Honshu Island, Japan. However, this region is not clearly indicated in the Figure. As this region is very important in the present study, it should be clearly indicated to put the readers in context and to ensure that they understand exactly which region corresponds to Tohoku region.

o This Figure is extremely similar to Figure 1 of Takamada et al. [13]. In fact, the maps used are the same. I think, the authors should include in the legend something like: “Modified from Takada et al. [13]”.

- Figure 2: In the Figure 2 (A) appears a total of 16 haplotypes. However, in Table S1 and in the manuscript (for example, line 194) is stated the existence of 15 haplotypes, not 16. Specifying, the haplotype Hja_033 appears in Figure 2 but is missing in Table S1. Could the authors explain this fact? Please, correct this properly in the paper.

- Table 1: Please, include “(n)” after “sample size”

- Table 2: Please, replace “two populations” with “two group of populations” in Table legend. See comments for Line 211.

- S1 Table: “total number (n)”. Total number of individuals? Please, indicate this in the Table caption.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Denis Copilas-Ciocianu

Reviewer #2: Yes: M. Pilar Cabezas

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

11 Sep 2019

Response to Reviewer #1

Phylogenetic substitution rate used in the comparison. The reviewer suggested to use inter- or intraspecific rate for other crustaceans such as 0.7% (Knowlton &Weight 1998) and 1.77% (Copilas-Ciocianu et al. 2019), instead of “general” evolutionary rate of 1% we used in the comparison. It is true that the rate of 1% is not estimated for crustacean species and we have reanalyzed the data with 0.7%/my by Knowlton &Weight (1998). The point in the figure 4 (and 5) is the disparity between general (or easy) rate and CDT rate, and also an uncertainty which is properly expressed by using stochastic clock model. Because showing the results with three (or more) different rates would increase the complexity of the figure (please check attached file Fig3-2.pdf), we would like to show BSPs with two rates. Ultimately a choice of the clock is arbitrary process and it is out of the focus.

Line 41: The statement has changed following the reviewer’s comment.

Line 153: We performed R2 test according to the reviewer’s suggestion and have added the statement regarding the analysis in Materials and methods (Line 153) and Results (Line 220 and Table 2) sections.

Line 159: We have added the statement that “(Hereafter ‘rate’ means substitution rate, not ‘divergence’ rate)” following the reviewer’s suggestion that it should be clarified what the rate means (substitution rate or divergence rate).

Line 160: As the reviewer suggested, “TEMs were built using all individuals from the NWP populations” is just what we mean. We however do not know to which of the five populations the Tohoku population assigned at this point. The statement has therefore been rephrased as follows: “we built TEMs using all individuals from Tohoku as well as its parental population among currently known sandy beach amphipod populations.”.

Lines 163, 231, 243: We have included the result for the strict clocks following the reviewer’s comment, while the result for the exponential TEM with strict clock model is not shown because the run did not converge.

Line 169: “9000 mutation values” has been replaced with “9000 values of the mutation-scaled transition time” for clarity.

Lines 175: We used ‘fitdistr’ function in the MASS v7.3 package, and we have included the name of the package in Line 175.

Line 99, 179, 252: Following the reviewer’s comment, ”BSP” in Line 99 has been spelled out as “Bayesian Skyline Plot” and reference added. Then the statement in Line 179 has been rephrased as follows: “We inferred the demographic history of the Tohoku and its parental population using a BSP.” For the clock model used in the BSP analysis, we stated that “The molecular rate was modeled based on the statistical distribution from CDT.” in Line 183, because the model is not determined at this point. For better clarity, we have replaced the phrase with “The molecular rate was modeled following the best clock model and statistical distribution from CDT”, also added the statement in Result section (Line 255) to show that exponential relaxed clock was applied to BSP.

Line 340: We have added the references for the faster evolutionary rate of crustacean species according to the reviewer’s suggestion. The point is the rate variation among crustacean taxa and we do not discuss each rate here and just introduced them as examples.

Response to Reviewer #2

Scientific name of the sandy beach amphipod. As suggested by the reviewer, generic name of ‘Eohaustorioides’ is also known for a single species of Haustorioides japonicus Kamihira 1977, although it seems not to be widely accepted. The genus Euhaustrioides was originally introduced in a book chapter without detailed description (Bousfield and Tzvetkova 1982), and the record of WoRMS is solely based on this information. Later Jo (1988) and Kamihira (1999), the authority of this amphipod group, refuted the establishment of Eohaustorioides, and the status of the genus seems still controversial. Our research team is now preparing a taxonomic paper with additional refutation on using this name. Furthermore, WoRMS is a web database and has nothing to do with scientific justification, while it argues that the database is controlled by “taxonomic experts”. We therefore avoided using Eohaustorioides in this manuscript in the first round. If it is further required to use Eohaustorioides in this manuscript, we would like to follow the editor/reviewer’s suggestion.

Line 86: We have corrected the statement following the reviewer’s suggestion.

Line 95: We have rephrased the statement as “genetically distinct from other known populations” following the reviewer’s comment.

Line 97, 98: “BSP” has been spelled out as “Bayesian Skyline Plot”, and the phrase “of the sandy amphipod Haustorioides japonicus” has been added after “genealogy”, according to the reviewer’s suggestion.

Line 111: “two sites” has been replaced with “one site”.

Line 116: “a” is inserted after “using” as suggested.

Line 142: GenBank accession numbers are added following the reviewer’s recommendation.

Line 143: Number of sequences are corrected.

Line 216: We have added “Fig 1” after “site 7” following the suggestion.

Lines 219, 223: The phrase “two populations” was replaced with “two group of populations”.

Line 303: The authority is added after the name of the species.

Lines 428, 434, 488: Formatting mistakes in the Reference section has been corrected according to the reviewer’s comments.

Figure 1 has been revised following the reviewer’s comment.

S1 Table and Figure 2: Haplotype Hja_033 has been removed from the figure because it is haplotype from the different genetic population (SJN, northern Sea of Japan) and not included in the analysis. The caption in the Table has been revised to show that “(n)” means total number of individuals.

Table 1 (Line 129): “(n)” is included after “sample size”.

Submitted filename: Response to reviewers.docx

Click here for additional data file.

25 Sep 2019

Reconstructing the population history of the sandy beach amphipod Haustorioides japonicus using the calibration of demographic transition (CDT) approach

PONE-D-19-16483R1

Dear Dr. Sakuma,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Tzen-Yuh Chiang

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: Authors perfectly addressed the comments made by the two reviewers. The manuscript is now acceptable for publication.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Denis Copilas-Ciocianu

Reviewer #2: Yes: M. Pilar Cabezas

30 Sep 2019

PONE-D-19-16483R1

Reconstructing the population history of the sandy beach amphipod Haustorioides japonicus using the calibration of demographic transition (CDT) approach

Dear Dr. Sakuma:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Tzen-Yuh Chiang

Academic Editor

PLOS ONE