MATERIALS AND METHODS

2.1 Population sampling

We sampled populations of the mud-tidal snail Batillaria attramentaria (G. B. Sowerby I, 1855) from two sites in Japan and one site in the USA (details in Appendix 1: Table A1). We also include here published data which we previously collected from a South Korean population (Ho et al., 2019a). Population locations comprised Hacheon, Cheollabuk-do, South Korea (on June 2016 at 35°32’N, 126°33’E, Ho et al., 2019a); Nemuro city, Hokkaido Prefecture, Japan (June 2017 at 43°15’N, 145°28’E); Matsushima Bay, Miyagi Prefecture, Japan (May 2018 at 38°22’N, 14 1°4’E); and Monterey Bay, Elkhorn Slough, CA, USA (February 2017 at 36°49’N, 121°45’W) (Figure 1). For each of the three native populations (one Korean and two Japanese sites), 100 individuals were collected. These native sites all had high surface salinities of 29-33 PSU at the time of collection. For the introduced (USA) population, 50 individuals were collected. This site had a low surface salinity of 4 PSU at the time of collection.
Due to river discharge, topography, and tides (e.g., Yoon & Woo, 2013), different estuaries can experience wide fluctuation in daily salinity or very little fluctuation. We characterized salinity fluctuation profiles for each sampling site as low (Nemuro-Japan site; 27-34 PSU), moderate (Hacheon-Korea site; 16-30 PSU), or high (Elkhorn Slough-USA site; 0-30 PSU) based on publicly available data on salinity fluctuation collected over the past several years. These data were obtained from http://www.nemuro.pref.hokkaido.lg.jp for Nemuro (Japan), http://www.khoa.go.kr for Hacheon (Korea), and http://www.mbari.org for the Elkhorn Slough (USA). We could not find recent records of salinity fluctuation at Matsushima Bay (Japan); however, past data on average monthly salinity levels inside Matsushima Bay indicates that salinity at this site fluctuates from 27 to 34 PSU (Ventilla, 1984) and is similar to the Nemuro site.
In animal locomotion experiments, it can sometimes be advantageous to use individuals of similar size. However, we noticed in the field that the typical body size of the introduced population was larger than the native populations. Since our primary question centers on how different populations respond to salinity stress, and since body size might play a role in both locomotion and salt tolerance, we allowed body size to differ between the native and introduced populations. This decision somewhat constrains our ability to conclude whether body size is a driving factor behind salinity responses; however, it allows us to use a representative sample of each population rather than using specimens whose size is not typical of their population, and thus makes our study more ecologically relevant. Ultimately, the native and introduced samples we collected differed by about 1 cm (average native shell length = 2.1 cm; introduced shell length = 3.1 cm). All collected specimens were maintained in a plastic aquarium with a constant air temperature of 25°C, a water salinity that was the same as their collection site, and a 12h L:12h D photocycle for two days prior to the experiments, in order to reduce the effects of transportation stress and to allow for acclimation to the lab.

2.2 Salinity stress experiment

For each population consecutively, we conducted a 30-day experiment examining locomotor behavior under different salinity treatments. We randomly divided each population into five treatment groups, with 20 individuals per group for the native populations and 10 individuals per group for the USA population. These groups were maintained in separate plastic aquaria (40 × 23 × 21 cm3, with an inclined layer of sea sand set up on the bottom and one liter of aerated artificial saline water, Supplemental Figure S1, Ho et al., 2019b) at salinities of 43, 33, 23, 13, and 3 PSU for 30 days. Saline water was freshly prepared every two days from overnight-aerated distilled water and Instant Ocean Sea Salt (United Pet Group Inc., Cincinnati, OH, USA). Snails in each group were marked with nail polish (Eco Nail color, Innisfree, South Korea) to keep track of their identity. All animals were fed to satiation every two days with a commercial brand of fresh, chopped seaweed (Ottogi, South Korea) throughout the 30-day experiment.

2.3 Locomotor behavior tracking

We recorded the movement of each snail for one hour every two days throughout the 30-day experiment following the protocol in (Ho et al., 2019a). Briefly, we used a Sony NXCAM camera (AVCHD Progressive MPEG2 SD) to film individual snails in the center of a disposable Petri dish (diameter: 9 cm) filled with artificial seawater which had the same salinity as the snail’s assigned treatment group. We increased the video playback rate using AVS Video Editor v.7.1.2.262 and cropped videos using Avidemux v.2.6.12. We used the Spectral Time-Lapse (STL) toolbox (Madan & Spetch, 2014) implemented in Matlab release R2014a (MathWorks Inc., Natick, MA, USA) to quantify movement distance of the snails.

2.4 Shell measurements

We measured shell length of all individuals using images extracted from the recorded videos. We used K-Multimedia Player software (KMP Player, PandoraTV, Pankyo, Korea) to extract one video frame (resolution about 1200 × 1200 pixels) for each snail, and then used tpsDig2 to digitize the most anterior and posterior points of the shell (Rohlf, 2008; Appendix 1: Figure A1). The diameter of the petri dish (9 cm) was used as the conversion scale to calculate snail length.

2.5 Sequencing and phylogenetic analysis

After the salinity experiments, we extracted genomic DNA from fresh foot tissue of all snails using a Dneasy Blood & Tissue kit (Qiagen, Hilden, Germany). We PCR-amplified the mitochondrial CO1 gene using published CO1 primers (Ho et al., 2015) and a Fastmix/Frenchetm PCR kit (IntronBio, South Korea). PCR products were purified using a Dr. Prep kit (MGmed, Seoul, South Korea). Sequencing reactions were performed using a Bigdye Terminator V3.1 Cycle Sequencing kit (Bionics; Seoul, South Korea).
We performed a Bayesian phylogenetic analysis based on the CO1sequences of five representative specimens which were identified as Kuroshio and Tsushima haplotypes and from Korea, Japan, and the USA (GenBank accession no.: MG241503-06 and MT800763), and 53 previously sequenced specimens from the shorelines of Korea (GenBank Accession: No. HQ709362– 81, Ho et al., 2015) and Japan (GenBank Accession: No. AB164326-58, Kojima et al., 2004). We used a closely related species,Batillaria multiformis, as an out-group (GenBank Accession: No. AB054364, Kojima, Ota, Mori, Kurozumi, & Furota, 2001). We employed MegaX (Kumar, Stecher, Li, Knyaz, & Tamura, 2018) to predict the best substitution model for the CO1 data, and found that the Hasegawa-Kishino-Yano + Gamma + Invariable (HKY85+G+I) substitution model was the best fit based on its corrected Akaike information criterion (AICc) value of 1500.65. We then ran the analyses applying Maximum Likelihood statistical method (Nei & Kumar, 2000) using MegaX with 1000 replications of bootstrap and Neighbor-Joining statistical method (Saitou & Nei, 1987; Studier & Keppler, 1988) using Geneious tree builder incorporated in Geneious v.6.1.8 (Kearse et al., 2012) with 1,000,000 of random seed to construct a phylogenetic tree.

2.6 Model terms

Our main purpose was to assess the effect of differential salinity exposure (es) on the locomotion of B. attramentaria (see the description of the 30-day experiment above). In addition to es, we also included in our analyses four other model terms that might influence snail response to salinity: origin (o), location (lo), population (p), and CO1 lineage (li). Origin was defined as either native (pooled three populations from Korea and Japan) or introduced (one population from the USA). The location term indicated the countries where the snails were collected (Korea, Japan, or the USA). The population term described the sampling site (Hacheon in Korea, Nemuro city and Matsushima bay in Japan, and Elkhorn Slough in the USA). Lineage was defined as either Tsushima (comprising the Hacheon and Nemuro populations) or Kuroshio (comprising the Matsushima and Elkhorn Slough populations) based on the individual’s position in the CO1 phylogenetic tree.

2.7 Statistical analyses

For all four populations, snails exposed to 3 PSU did not move at all and died within the first 16 days of the experiment, apparently due to the extreme osmotic stress represented by such low salinity. Locomotor data for all 3-PSU groups were therefore excluded from all analyses. All other individuals survived the entire duration of the experiment and were included in the analyses.
We applied a linear mixed-effect model (LMM) to assess the impacts of multiple predictor variables (see “Model terms” above) on the locomotor response of B. attramentaria to different levels of salinity stress using the package nlme version 3.1-140 (Pinheiro, Bates, DebRoy, Sarkar, & R, 2018) implemented in R version 3.0.2 (R Development Core Team, 2011). Since we wished to investigate specifically the impacts of salinity stress, geographic distribution, and genetic composition, we fit a set of competing models with separate single predictor variables (es, o, t, lo, p, and li) and their interactions (es × o, es × o + li, es × lo, es × lo + li, es × p, es × p + li, es × li, es × li + o, es × li + lo, and es × li + p). We treated snail identity as a random variable. For the sake of simplicity, we limited the maximum number of predictor variables to three per model. Since body size closely corresponded to the origin factor, which was already in our LMM analysis (i.e., we knew native snails were smaller than introduced snails), we did not include body size as a factor in the LMM analysis. We centered all the predictor variables to mean = 0 and standardized to s.d. = 0.5 to remove multicollinearity and to directly interpret the results in terms of effect size (allowing us to compare predictors). The response variable was the movement distance of snails, which was measured every two days throughout the 30-day acclimation experiments. Prior to the analysis, we square-root transformed the movement data using the package rcompanion 2.2.1 (Mangiafico, 2017) implemented in R 3.0.2 (R Development Core Team, 2011) to meet the assumption of normal distribution. To determine the best covariance structure for the LMM tests, we tested our response variable against several covariance structures: First Order Autoregressive (AR1), Compound Symmetry (CS), and Unstructured (UN). We compared their corrected AICc values using the package MuMIn 1.9.13 (Barton, 2009) for R 3.0.2. The AR1 model was the best-supported covariance structure based on the AICc value (3407.24 versus 3710.22 (CS) and 3708.22 (UN)) and was therefore chosen for the LMM tests as the best available compromise between bias and lack of precision.
We next applied a multimodel inference procedure (Burnham & Anderson, 2003) to the set of competing linear mixed-effect models (LMMs) to select the most parsimonious model that best described snail locomotor response. The models were compared based on AICc values using the aforementioned MuMIn package. The model with the lowest AICc value and those satisfying a ΔAICc ≤ 6 cut-off rule (Richards, 2005, 2008) were considered the best-fit or most parsimonious models. We then conducted post-hoc multiple comparison tests of the best-fit model to examine the effects of each explanatory factor. Additionally, we used MuMIn to perform model averaging and estimate the importance of predictor variables by summing the weights of models where the variables appeared.
In parallel, we examined differences in shell length between native and introduced groups and among native populations using a two-way ANOVA. The first level involved comparisons of snails from four populations, three locations (i.e., countries), and two origins (introduced vs. native). The second level involved comparing individuals belonging to the two mitochondrial CO1 lineages of Kuroshio and Tsushima.
2.8 Data deposition
Data available from the Dryad Digital Repository:https://datadryad.org/stash/dataset/doi:10.5061/dryad.455mv2mand the Mendeley Data repository:http://dx.doi.org/10.17632/jjjmh26c2g.1.