Genetic structure and shell shape variation within a rocky shore whelk suggest both diverging and constraining selection with gene flow

Variation in snail shell shape has provided evolutionary biologists with excellent material for the study of local adaptation to local environments. However, the assumption that shell shape variation is evidence of distinct lineages (species) might have led to taxonomic inflation within some gastropod lineages. Here, we compare shell shape variation and genetic structure of two independent lineages of New Zealand rocky shore whelks in order to understand the process that leads to an unusual disjunct distribution. We examined the Buccinulum vittatum complex (three subspecies plus Buccinulum colensoi) using mitochondrial DNA sequences, 849 single nucleotide polymorphisms and geometric morphometric data on shell shape. Specimens within the range of B. colensoi shared nuclear markers and mitochondrial DNA haplotypes with the southern subspecies Buccinulum vittatum littorinoides, whereas the northern samples (Buccinulum vittatum vittatum) had distinct genotypes. We infer that gene flow between B. colensoi and B. v. littorinoides is greater than between either of these taxa and B. v. vittatum. The distinctive shell sculpturing of B. colensoi is maintained despite gene flow, suggesting a role for selection. Although the shell form of B. colensoi appears to be an adaptation to local conditions, we did not observe convergent evolution in the sympatric whelk species Cominella maculosa.


INTRODUCTION
The geographical distribution of related species is used by biologists to infer evolutionary processes (Warren et al., 2014) and is fundamental to the study of biogeography and speciation.The geographical distributions of species have been used to make inferences about speciation mechanisms, competitive exclusion and niche evolution (Webb et al., 2002;Anacker & Strauss, 2014).For example, a widespread species with its range interrupted by the occurrence of a congeneric species might be interpreted as evidence of competitive exclusion (Emerson & Gillespie, 2008).However, when species classification is based on morphological traits under selection we can be misled about the evolutionary relationships of populations, and thus the distribution of these putative taxa could lead to erroneous conclusions.Species distributions should be based on taxonomy that accurately reflects evolutionary relationships, essential for inferences about the processes of biogeography.
Many snail lineages exhibit adaptation to local environments, resulting in the formation of ecotypes (e.g.Davison, 2002;Teshima et al., 2003;Quesada et al., 2007;Stankowski, 2013;Wada et al., 2013;Butlin et al., 2014).The phenotype of a snail shell responds to This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons. org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 a wide variety of genetic and environmental influences.In marine environments, factors such as water depth, temperature, currents, wave energy and predation have been shown to affect plastic shell traits (Palmer, 1990;Baker et al., 2004;Olabarria & Thurston, 2004;Shelmerdine et al., 2007;Hollander & Butlin, 2010;Hills et al., 2012;Avaca et al., 2013;Butlin et al., 2014;Haig et al., 2015;Irie & Morimoto, 2016; reviewed by Bourdeau et al., 2015).Different ecological conditions across the geographical range of a species can generate polymorphisms that can lead to the erroneous taxonomic over-splitting of a lineage (Hills et al., 2012;Donald & Spencer, 2015).Admixture after range expansion of populations previously isolated can also create spatial patterns of morphological variation (Sivyer et al., 2018).Conversely, constraining selection that maintains traits associated with certain conditions can lead to similarities among taxa from separate lineages.Where species are closely related, this can result in cryptic species, another cause of problems in taxonomy (Allmon & Smith, 2011).Landmark-based geometric morphometric techniques allow for the empirical comparison of complex shell shapes and can remove the confounding effects of size.By decomposing shape variation within a dataset into linearly uncorrelated variables [principal components (PCs)], model-based cluster analysis provides an unbiased and objective tool to cluster and assign specimens to groups naively, based on shape.This combination of methods to describe shape variation and to determine the most efficient model for clustering variation among specimens adds empiricism and statistical independence to studies of natural phenotypic variation.
In New Zealand, the marine gastropod genus Buccinulum Deshayes, 1830 comprises species that are morphologically similar and, in some cases, difficult to distinguish (Ponder, 1971;Marshall, 1998).Buccinulum vittatum (Quoy & Gaimard, 1833) is a common rocky shore whelk with an unusual distribution.The species occurs the full length of New Zealand but is absent from a ~400 km stretch of eastern coast, where it is interposed with the congeneric Buccinulum colensoi (Suter, 1908) (Fig. 1).The exclusion of one species of rocky shore whelk by another is possibly attributable to spatially abrupt differences in environmental conditions that facilitate interspecific competitive exclusion.Alternatively, the local environment might have selected for an ecotype of B. vittatum that has mistakenly been assigned a separate taxonomic status.Here, we examine morphological and genetic variation to investigate these competing explanations.An ecotype would not be expected to be differentiated strongly at neutral genetic loci from other populations, whereas a separate species would show concordance between shell shape and genetic markers.Habitatspecific convergent evolution provides a framework for testing adaptive hypotheses (Minards et al., 2014), because similar responses in separate lineages indicate that traits are adaptations resulting from local selection (Butlin et al., 2014).To determine whether the local environment can explain the particular shell phenotype recognized as B. colensoi, we also examined an independent lineage from the same environment.The rocky shore snail Cominella maculosa (Martyn, 1784), although superficially similar in appearance, is not closely related to Buccinulum (Donald et al., 2015;Vaux et al., 2017).Both species lay eggs that hatch as small snails (Ponder, 1971;Morley, 2013;Donald et al., 2015), with no pelagic larval stage, and so are expected to have similar dispersal rates.The direct larval development of a marine species can facilitate the maintenance of biogeographical patterns (Lee & Boulding, 2009) when water currents are considered (Schiel, 2004;Sponaugle et al., 2005).On the eastern coast of North Island New Zealand north of East Cape, the coastal ocean current flows south until reaching East Cape (where the distributions of B. colensoi and B. vittatum meet), where the current moves offshore.In contrast, south of East Cape the coastal current flows north (Heath, 1985).Both lineages studied here are intertidal predatory whelks and co-occur in the same rock pools along the same coastline (Ponder, 1971).Matching spatial distributions of shell shape variation in C. maculosa and the B. vittatum complex would indicate a convergent adaptive response to local environmental conditions and help us to understand the distinctive shell form of B. colensoi.

Taxonomic background
Buccinulum and Cominella species are small (10-40 mm shell length) members of the globally distributed gastropod family Buccinidae (Bouchet et al., 2005(Bouchet et al., , 2017)), which in New Zealand also includes Aeneator Finlay, 1926, Antarctodomus A. Adams, 1863, Austrofusus Kobelt, 1879, Euthrenopsis Powell, 1929, Pareuthria Stebel, 1905and Penion Fischer, 1884(Powell, 1979;Spencer et al., 2016).Most of the 14 currently recognized species of Buccinulum are readily distinguished from each other, but a group of taxa endemic to New Zealand, referred to here as the B. vittatum complex, are not so easily distinguished.The key traits used to classify these taxa are shell shape, spiral sculpturing and location (Ponder, 1971).The B. vittatum complex comprises taxa that have all, at some time, been considered subspecies of B. vittatum (Ponder, 1971).Recognized in this study are : B. v. vittatum, B. v. littorinoides (Reeve, 1846), B. v. bicinctum (Hutton, 1873) andB. colensoi. Recently, B. v. littorinoides has been viewed as a synonym of B. vittatum (Spencer et al., 2016;following Marshall, 1998), and B. colensoi as a discrete species (Spencer et al., 2016;following Powell, 1979).These taxa are separated primarily on the basis that they form distinct geographical clusters (Fig. 1), but the B. colensoi shell morphotype is generally broader, more robust and with stronger spiral sculpturing than the other members of this complex and has rugose external texturing (Ponder, 1971; Fig. 1).Buccinulum v. vittatum is found north of B. colensoi, whereas B. v. littorinoides is found south of B. colensoi, on the east coast of North and South Island New Zealand, and on the coast of more southern islands within the New Zealand region (Stewart, Snares and subantarctic islands; Fig. 1).Buccinulum v. bicinctum is restricted to the Chatham Islands (Fig. 1).
Nine Cominella species are found in New Zealand waters (Powell, 1979;Donald et al., 2015;Spencer et al., 2016).The species C. maculosa is sympatric with the B. vittatum complex (Fig. 1), and recent mitochondrial DNA (mtDNA) sequencing has revealed fine-scale changes in population sample haplotype frequencies along the east coast of North Island, suggestive of the mixing of northern and southern haplotypes of this species (Dohner et al., 2018;Fleming et al., 2018).

Sample collecTion
Specimens of Buccinulum and C. maculosa were collected by hand from rocky shore locations around the coast of New Zealand.Freshly caught live snails were frozen and transferred to the laboratory, where they were thawed, extracted from shells and preserved in ample 99% ethanol.Shells were retained for morphometric analysis (Table 1).Specimens were identified using gross shell morphology and collection location (Ponder, 1971;Powell, 1979).

geomeTric morphomeTric daTa
To produce a geometric morphometric dataset, digital images of shells were obtained following the approach of Dowle et al, (2015) and the recommendations of Collins & Gazley (2017).Shells were placed in a bed of contrasting coloured sand with ventral surface upwards and positioned so that the aperture was horizontal.Photographs were taken using a Canon EOS 600D with EF100 mm f2.8 USM macro lens mounted on a high-precision Kaiser stand.Two digital 'combs' were positioned over images of each shell using Adobe Photoshop CS6 (see Supporting Information, Fig. S1).Shell morphology was analysed using a twodimensional, landmark-based geometric morphometric approach with a combination of fixed landmarks and semi-landmarks.Biologically homologous points at the tip of the spire and base of the columnar and at the edges of whorls provided positions that are interpreted Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 as type 1 landmarks (sensu (Bookstein, 1991;Gunz et al., 2005; see Supporting Information, Fig. S1)).Sliding semi-landmarks were positioned along aperture and outline curves using the digital combs and slid computationally to minimize bending energy, effectively shifting them along the curve in order to minimize effects of arbitrary placement on the curve (Zelditch et al., 2004).Shell shape analysis using this method is effective in discovering variation between ecotypes, species and genera (Hills et al., 2012;Collins et al., 2013;Dowle et al., 2015;Vaux et al., 2018).Digitization was undertaken in tpsdig2 v.2.17 (Rohlf, 2013) on a Wacom Cintiq 22HD digitizing tablet.Digitized semi-landmarks were slid using SEMILAND7, part of the imp714 package, implementing the Procrustes distance method.Landmark x-y coordinates were then imported into morphoJ 1.05f (Klingenberg, 2011) for nonparametric statistical analysis.
A set of five fixed and 21 semi-landmarks were digitized on the curves of a shell for morphometric analysis (see Supporting Information, Fig. S1).Higher numbers of landmarks were trialled but provided no noticeable improvements in resolution based on principal component loadings.The lower number of landmarks also helped to reduce degrees of freedom that result from having more landmarks than samples, thus reducing the likelihood of type I errors that are more likely with datasets with high numbers of landmarks.
To quantify the amount of methodological error introduced into analysis from the placement of the shell for photography and digitizing photographs, the disparity of a set of replicates was calculated and compared with the variation found in the dataset as a whole.Photographic variables, such as camera height and the position of the shell within frame, have been shown to have negligible effect (Collins & Gazley, 2017).Disparity tests were performed in the R package Geomorph (Adams & Otárola-Castillo, 2013) implemented in the R programming environment (R Core Team, 2017).
Shell shape variation was examined using the javabased program morphoJ 1.06f (Klingenberg, 2011) and assessed using principal components analysis (PCA) across all individuals and all landmarks.A brokenstick test using Eigenvalues produced in morphoJ was implemented using the vegan 2.4-1 package (Oksanen et al., 2015) in R, to determine how many statistically significant, informative principal components were present.Canonical variates analysis (CVA) was used to test the discrimination of groups based on the current taxonomic treatment.Discrimination success was estimated using cross-validation scores (number of individuals correctly assigned to each group with 1000 permutations).
For each genus, morphometric data were examined using a model-based clustering approach to explore the distribution of variation and to test for natural clusters without a priori classification.We used the MCLUST v.5.0.2 package (Fraley et al., 2012) in R that implements Gaussian modelling where the total dataset is considered as a mixture of multivariate normal datasets, with a selection of covariance structures and vectors of expectation (Nanova, 2014;Scrucca et al., 2017).Unlike discriminant analysis (CVA), MCLUST analysis does not require prior information about specimen identity to classify sample data (Fraley & Raftery, 1999, 2002, 2003).The optimal model and number of clusters in the data are selected based on Bayesian information criteria (BIC), using the value of the maximized log likelihood, with a penalty for the number of parameters in each model (Fraley & Raftery, 1999, 2002, 2003;Cordeiro-Estrela et al., 2008;Nanova, 2014).In MCLUST, BIC scores are multiplied by minus one, so that higher BIC scores indicate lower global average and median classification uncertainty, and better model fitting (Fraley & Raftery, 1999;Cordeiro-Estrela et al., 2008).

miTochondrial dna Sequence daTa
DNA was isolated from snail foot tissue using either the Geneaid column extraction kit or a modified DNA extraction method using cetyltrimethyl ammonium bromide (CTAB) protocol involving the removal of unwanted organic compounds with chloroform:isoamyl alcohol (Doyle & Doyle, 1990;Trewick et al., 2009).In most cases, the CTAB method produced higher yields of genomic DNA than kit extraction.DNA quantity and quality were assessed using agarose gel electrophoresis and a Qubit broad-range DNA assay.
Partial sequences of the mitochondrial cytochrome oxidase c subunit I (cox1) were amplified using the primers HCO2198 and LCO1490 (Folmer et al., 1994).Polymerase chain reactions (PCRs) were performed in 20 μL volumes containing: 200 μM dNTPs, 2.5 µL NEB thermopol 10× buffer with 1.5 mM MgCl, 1 μM primers, 0.20 U of NEB taq DNA polymerase and 1-10 ng of template DNA.Standard thermal cycling conditions were followed, with 50 °C annealing temperature and 35 cycles, in a Biometra T1 thermocycler.The PCR products were sequenced using primer HCO2198 with BigDye Terminator v.3.1 on an ABI 3730 genetic analyser.The DNA sequences were visualized, checked for errors and ambiguities, and aligned in Geneious v.7 (Kearse et al., 2012).
Mitochondrial cox1 fragments for the two whelk lineages were separately aligned using the Geneious v.7 alignment tool with default settings.Alignments and protein translations were checked by eye for anomalies.Minimum spanning networks (Bandelt et al., 1999) were inferred for mtDNA sequences from both the Buccinulum complex and C. maculosa datasets in PopART v.1.7 (Leigh & Bryant, 2015).Phylogenetic analysis for the Buccinulum complex used MrBayes plugin 2.2.2 for Geneious 9.1.3with 10 000 000 iterations and a burn in of 1 000 000 (10%) using a GTR inv-gamma model (Tavaré, 1986) as selected by Jmodeltest (Posada, 2008).The species Buccinulum pallidum powelli Ponder, 1971 was used as an outgroup (Vaux et al., 2017).
Evidence of isolation by distance was sought using a Mantel test of the correlation of pairwise logarithmically transformed geographical distances and pairwise Φ ST /(1 − Φ ST ), with 100 000 permutations in genepop v.4.2 (Raymond & Rousset, 1995;Rousset, 2008).Population genetic differentiation is estimated using the correlation of random mtDNA haplotypes within a population relative to that of random pairs of haplotypes drawn from the whole sample (Φ ST ).Geographical distances between coastal sites were estimated by following the coastline, taking the shortest distance between locations, using the measurement tool in Google Earth.Pairwise Φ ST was estimated from the cox1 sequences in Arlequin v.3.5 (Excoffier & Lischer, 2010).

nuclear geneTic markerS
To detect population genetic structure, we chose to maximize the number of nuclear markers studied rather than maximizing the number of specimens, because this has been shown to be an efficient strategy for detecting structure (Landguth et al., 2012;Prunier et al., 2013).Anonymous single nucleotide polymorphic nuclear markers (SNPs) were generated for members of the Buccinulum complex using a double digest restriction enzyme DNA sequencing (ddRAD-seq) protocol, following the methods outlined by Peterson et al. (2012) with minor modifications.The restriction enzymes used for digestion of genomic DNA were NsiI-HF and MboI (NEB), which were selected after a trial of potential enzymes to determine the best cutting efficiency.Although MboI can cut the Buccinulum mitochondrial genome many times (sequences available on GenBank: MH198158-MH198162), NsiI-HF has only one cutting site per mitochondrial genome.Thus, only a single mtDNA fragment could be included downstream, but size selection would remove it.Whole genomic DNA of 78 Buccinulum specimens was digested separately.Barcode sequences were ligated to fragments using Invitrogen T4 ligase to enable identification of individuals after pooling of samples.Fragments between 250 and 350 bp long were selected for sequencing.Single-end highthroughput sequencing was performed using an Illumina Hi-seq (NZGL), and resulting data were processed using the STACKS 1.0.1 pipeline (Catchen et al., 2013).
Selection of nuclear markers was undertaken so that analyses would be performed on loci likely to be single copy and for which the maximum number of individuals could be genotyped (Harvey et al., 2015).In STACKS, a range of parameter settings were implemented relating to read coverage, individual number and population coverage.Initial exploration of the data suggested that depth coverage varied for the 73 individuals.Recommended read depth coverage settings vary in the literature (Peterson et al., 2012;Buerkle & Gompert, 2013); therefore, we experimented with parameter optimization.After initial tests using the STACKS pipeline to assess information content, samples with file sizes of < 2 megabytes were excluded from further analysis.These samples had low DNA read numbers, and their inclusion reduced analytical power downstream.Low coverage combined with high error rate has the potential to reduce the number of true loci by a substantial amount, up to 51% (Catchen et al., 2011).
The number of mismatches allowed between alleles when processing a single individual (−M) was explored.If −M is too low, some real loci are not formed, and subsequent alleles will be treated as different loci (undermerging).If −M is too high, repetitive sequences and paralogues will form large nonsensical loci (overmerging).We tested a range of values (−M = 5-25) and found that catalogue construction failed at high values (e.g.−M = 25).When alleles were allowed to differ by a maximum of 15 substitutions (−M = 5-15), some variation in F ST (the amount of genetic variance that can be explained by population structure, used here as a proxy for identifying population differentiation) was recorded, but the relationships among population samples remained the same among trials.
The number of mismatches allowed between alleles within a locus (amongst individuals) when building the catalogue (−N) was explored with a range of values (−N = 2-25).If alleles amongst individuals within a locus were not allowed additional sequence variation (−N = 0), there would be loci represented independently across individuals that are, in fact, alleles of the same locus.If N > 0, the consensus sequence from each locus is used to attempt to merge loci.This parameter is important for population studies where monomorphic or fixed loci may exist in different population samples.Merging fixed alleles as a single locus can increase the probability of assembling real loci and therefore decrease the allele error rate.However, erroneous loci will be created if −N is too high.It is recommended that the setting of −N should consider the genetic similarity expected among samples within each study.When −N = 25 was used, catalogue construction failed.When values for −N were between five and 15, some variation in F ST was recorded, but the relationships among population samples remained the same.
The de novo map pipeline for STACKS was used to build a catalogue using 28 Buccinulum individuals from Oakura, Waiheke Island, East Cape, Mahia, Castlepoint and Kaikoura (Table 1).Here, we report analysis using a minimum of five reads per individuals (−M) as providing a reliable set of markers for downstream analysis and excluding all stacks with a lower coverage.Potentially spurious, highly repetitive stacks were removed.Within an individual, we allowed a maximum of three mismatches between alleles (−M) and five mismatches between primary and secondary reads (−N = M + 2).Ten mismatches were allowed between alleles in different individuals (−N) when generating the SNP set.The final parameter settings were as follows: USTACKS, −M 5, −N 5, −M 3; CSTACKS, −N 10.Analysis was restricted to a single SNP site per putative locus (always the first), avoiding potential problems of data non-independence.
A range of values for the number of populations a SNP marker was required to be present in before being recorded were tested in POPULATIONS, which is part of the STACKS pipeline allowing computation of population-level summary statistics, and the output of site-level SNP calls for subsequent analysis in STRUCTURE v.2.3.4 (Pritchard et al., 2000;Hubisz et al., 2009).The resulting number of loci returned varied considerably when tightening or relaxing the stringency at which the presence of an SNP was required in a number of putative populations or proportion of individuals in a population.The dataset required that each putative locus included was genotyped in individuals from one or more of the six population samples (-r 1) and genotyped in ≥ 40% of individuals within those population samples (−r 0.4).This resulted in allelic data for ≥ 12 of the 28 specimens for each locus.The minor allele frequency was set at 0.1 (-min_maf 0.1), removing rare alleles from the dataset.This resulted in the following final POPULATION settings: -p 1 -r 0.4 -f p_value --write_single_snp --fstats --min_maf 0.1-k.
The resulting SNP dataset was analysed using Bayesian clustering in STRUCTURE to identify population differentiation.Ten replications of the admixture model with independent allele frequencies, Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 with a burn-in of 200 000 and 1 000 000 generations, were used.Potential populations (K) was set from one to six (the number of population samples).To determine the optimal number of clusters, we examined estimates of the posterior probability of the data for a given K (Pritchard et al., 2000) and ΔK, the rate of change in logarithmic probability of the data (Evanno et al., 2005) implemented in STRUCTURE HARVESTER v.0.6.8 (Earl & vonHoldt, 2012).Pairwise F ST as an indicator of population differentiation was estimated from 849 nuclear SNPs using the genepop R package (Rousset, 2008).Evidence of isolation by distance was sought using a Mantel test of the correlation of logarithmically transformed pairwise geographical distances and pairwise F ST /(1 − F ST ), with 1000 permutations, in genepop v.4.2.Geographical distances between coastal sites were estimated by following the coastline, taking the shortest distance between locations.

RESULTS
Shell Shape variaTion uSing geomeTric morphomeTricS Shell shape variation was studied with fixed and sliding landmarks on 88 Buccinulum snail specimens from 11 locations.A disparity test indicated that errors introduced from the digitization process, comb and landmark placement (0.3% of total shape variation), and total variation and variation from shell placement and digitization (1.4%) were negligible, and therefore these sources of error were ignored in subsequent analyses.
The first five principal components obtained from a PCA across all Buccinulum specimens provided statistically significant information (broken-stick test) and explained 84.4% of shell shape variation (see Supporting Information, Table S1).When the first two principal components (63.8% variation) were used to Figure 2. Shell shape variation within two sympatric marine snail lineages: Buccinulum and Cominella.Principal components analysis of landmark-based geometric morphometrics of shell shape, with the colour of the specimen coded by collection region (coastal New Zealand).Ellipses are 95% confidence ellipses on group means and indicate some regional partitioning of morphological variation.The thin plate spline warp grid diagrams show the shape variation captured from shells by principal component (PC)1 and PC2.visualize shape variation, there was some clustering of specimens, but overlapping shell shape distributions (Fig. 2).Buccinulum specimens from the eastern region (the range of B. colensoi; Fig. 1) were the most distinct for PC1.Clusters of shell shape for other regions were also apparent, with non-overlapping confidence regions on group means, but without clear separation of the four a priori groups.Not all individuals were allocated to their original sample groupings in cross-validation tests (P-values for difference between means ranging between 0.37 and 0.97; Table 2).A single morphological cluster of specimens was inferred as optimal for the Buccinulum data using the model-based clustering approach with principal components 1-5 (diagonal multivariate normal model).
Given that the majority of shell shape variation was explained by PC1 (51.7%;Fig. 2), we explored the data further using a model-based clustering of specimens performed with only PC1.The Bayesian information criterion selected a model with two morphological clusters as being the optimal clustering strategy for PC1 (Fig. 3).One cluster contained 15 specimens collected from the eastern region (geographical range of B. colensoi) and two specimens from the northern region (Oakura).These shells were characterized by relatively short spires and wide apertures (Fig. 2).The other cluster comprised seven specimens from the eastern region (Mahia, Akitio and Castlepoint) with samples from the northern, southern and Chatham samples (assignment probabilities for each specimen are shown in Fig. 3).
Shell shape variation within the sympatric whelk C. maculosa was examined using 46 shells from six sampling locations.The first three principal components provided statistically significant information (broken-stick test) and explained 73.74% of shell shape variation (see Supporting Information, Table S2).When the first two principal components (65.26% of variation) were used to visualize shape variation, little clustering of regional samples was evident.There was overlap of regional samples in morphospace, although eastern region shells had a 90% mean confidence ellipse that did not overlap the ellipses for the other two regions (Fig. 2).Snails collected from the eastern region had shells that were, on average, tall and narrow (Fig. 2).Nevertheless, a single morphological cluster was identified by modelbased clustering of PC1 and from the first three principal components using MCLUST.

miTochondrial dna diverSiTy and STrucTure
A fragment of cox1 sequence (443 bp) was used for phylogenetic inference of B. v. vittatum, B. v. littorinoides, B. v. bicinctum and B. colensoi, with the outgroup B. powelli.Within the B. vittatum complex, three mtDNA clades were resolved (see Supporting information, Fig. S2), one representing the northern subspecies B. v. vittatum, which was sister to a clade of Chatham Island specimens (B.v. bicinctum), and a third clade with the southern subspecies B. v. littorinoides and B. colensoi (maximum distance between haplotypes ~10.6%; see Supporting information, Fig. S2).An alignment of partial cox1 DNA sequence (443 bp) was compiled for 92 individuals from the B. vittatum complex, with 361 identical sites (81.3%).We excluded the Chatham Island samples (B.v. bicinctum), given their geographical and genetic separation (Supporting information, Fig. S2).
Variation was expressed in a haplotype network (Fig. 4).Mitochondrial DNA haplotypes from B. v. vittatum formed a cluster, and haplotypes from B. v. littorinoides and B. colensoi grouped together (Fig. 4).Identical mtDNA haplotypes were observed in whelks with the shell shape and sculpturing of B. colensoi and B. v. littorinoides but collected from different locations.Correlation between geographical distance and mtDNA genetic distance (pairwise Φ ST ) was significant, providing evidence of isolation by distance [Pr(correlation < observed correlation) = 0.0191; one tailed P-value; see Supporting Information, Fig. S3].
Partial cox1 (516 bp) was sequenced from 22 C. maculosa individuals from the same geographical range as the B. vittatum complex (excluding Chatham Islands).These sequences were added to 40 haplotypes from GenBank (Fleming et al. 2018), to produce a dataset representing 555 C. maculosa individuals, each coded into one of three geographical collection regions of New Zealand (North, East, South).No clear geographical pattern was observed, although the common haplotypes in the East were rare elsewhere.The Northern sample contained greater haplotype diversity than East or South, but one mtDNA haplotype was found in all three regions (maximum distance between haplotypes ~2.1%;Fig. 4).A set of 849 anonymous nuclear loci (SNPs) were produced for 28 individuals from the B. vittatum complex, from six locations (Table 1).The genetic clustering model with the best fit to the data (optimal K value = 2; see Supporting Information, Table S3 and Fig. S5) divided the Buccinulum specimens into two clusters (STRUCTURE HARVESTER; Fig. 3).Ten specimens (from Oakura, Waiheke Island and Hicks Bay) grouped with assignment probabilities between one and 0.69, and 18 specimens (from eastern and southern areas) grouped together with assignment probabilities between one and 0.79.This clustering of specimens was consistent with the mtDNA results, which indicated that the specimens from the northern region (B.v. vittatum) form a distinct genetic clade from eastern and southern population samples (B.colensoi and B. v. littorinoides).However, significant allele sharing was apparent in the genotypes of whelks from Hicks Bay, revealing gene flow between the northern and the eastern/southern cluster.In models with more genetic clusters (higher K values), B. colensoi specimens become separated from B. v. littorinoides, but with evidence of gene flow in their low assignment probabilities (see Supporting information, Fig. S4).Allele sharing of whelks from Hicks Bay with both northern and southern population samples means that these individuals cluster most closely with different populations at different K values (e.g. at K = 2 they cluster with northern samples; at K = 3 they cluster with southern samples; at K = 4 they form their own group; Supporting Information, Figs S4, S5).All populations showed substantial genetic differentiation (estimated with pairwise F ST ; Tables 3 and 4), with the exception of Mahia and Castlepoint (F ST close to zero; 0.032).Geographically adjacent population samples were generally more similar to each other than more distant samples, and the only population sample from the west coast of New Zealand (Oakura) was genetically the most distinct from all others (maximum pairwise F ST = 0.62).The correlation between geographical distance and genetic distance (pairwise F ST ) was significant, providing evidence of isolation by distance [Pr(correlation < observed correlation] = 0.0260, onetailed P-value; see Supporting Information, Fig. S3).
Pairwise genetic differentiation suggests that gene flow between all populations sampled is low, but with slightly higher gene flow between B. colensoi and  4).

DISCUSSION
The distribution of currently recognized taxa within the B. vittatum complex suggests that shell phenotype might be misleading us about evolutionary relationships.Despite the current taxonomy recognizing two distinct species and three subspecies within B. vittatum, initial geometric morphometric data did not resolve distinct (non-overlapping) morphological clusters.Specimens from within the range of B. colensoi have the most distinct shell shape.However, the most identifiable characteristic of the shell of B. colensoi is the strong surface sculpturing (Ponder, 1971), which was not captured by the landmark-based morphometric analysis of shell shape.We noticed within our sample that even these shell surface traits formed a cline over the range of B. colensoi, but we have not quantified this observation.Although our geometric morphometric analyses do not take account of sculptural differences used in traditional taxonomic treatments, the subtle shape differences measured here are expected to reflect genetic differences and phylogenetic relationships.Using principal components generated from our two-dimensional, landmark-based morphometric data, a single morphological cluster fitted the model better than a set of discrete clusters.However, when only PC1 was analysed, a cluster predominantly from the eastern region (corresponding to B. colensoi, but including some west-coast specimens of B. v. vittatum), was distinguished.This discrepancy between multiple PC and single PC analysis indicates that although the shape variation within each PC is statistically informative, the variation captured by them is not necessarily taxonomically informative or geographically constrained.Members of the eastern region morphological cluster are not exclusively from this region; thus, our analyses of shell shape variation did not show taxonomically diagnostic characteristics of shell shape that would be consistent with recognizing a separate species from the eastern North Island (B.colensoi).However, we did not include thickness of shell or surface sculpturing in these analyses.One possible explanation for the recognition of the East Coast species, B. colensoi morphotype, is that the local environment is distinctive and selects for a particular shell shape.For example, geological characteristics of the region could potentially affect sediment loads, substrate types or algae species present, which might influence shell shape (Martín-Mora et al., 1995;Bourdeau et al., 2015).The southern east coast of North Island is primarily sedimentary rock, although the characteristic broad intertidal platforms are probably not fundamentally different from many other New Zealand coastal regions (Healy & Kirk, 1992).Nevertheless, we sought evidence of convergent evolution by examining the shell shape of a sympatric whelk collected from the same coastal habitat.These two whelk taxa might have different microhabitat preferences, but their widespread sympatry is expected to facilitate similar responses to the same environmental variation.Convergent regional phenotypic differentiation was not supported by study of the sympatric whelk C. maculosa.Although, on average, specimens from the eastern region did  There is a morphotype of Buccinulum currently identified as B. colensoi present on the east coast of the North Island, but the shell traits that are distinctive to this taxon were not found to be diagnostic of a separate lineage.Variation within population samples suggest that the B. colensoi form might be better considered an ecotype.There is no clear indication of environmental conditions that might be influencing phenotype.In contrast to the limited shell shape variation, we found strong genetic clustering in Buccinulum, revealing the western and northern samples of B. vittatum as genetically differentiated from the southern samples.However, almost all population samples were genetically differentiated from others within these two genetic clusters.Relative rates of gene flow as inferred from F ST estimates suggest that B. v. vittatum has experienced less exchange with the other two taxa than B. colensoi and B. v. littorinoides have with each other, but isolation by distance of both nuclear markers and mtDNA suggest divergence with gene flow.The two genetic clusters resolved by nuclear and mtDNA markers suggest that a taxonomic revision is necessary.Our data indicate that B. colensoi and B. v. littorinoides are likely to be a single polymorphic cluster and that there is sufficient genetic differentiation to consider uniting them as a separate genotypic cluster (species), distinct from B. v. vittatum.We thus recognize the species B. littorinodes (the name has precedence over colensoi) as comprising both the southern form and the distinctive east-coast form.At this stage, we cannot provide a diagnostic morphological trait for distinguishing southern B. littorinodes from B. vittatum.This example illustrates that shell shape and sculpture can be misleading when identifying evolutionary lineages, and thus inferences drawn from species distributions should consider whether taxonomy reflects evolutionary relationships (Warren et al., 2014).
Studies of marine species have previously recognized a biogeographical break on the North Island of New Zealand at East Cape (Gardner et al., 2010), with several species having distributions restricted to this area (Marshall, 2014).The pattern of genetic structuring found here in the B. vittatum complex is similar to that found in Haliotis iris Gmelin, 1791 (pāua) (Will et al., 2011(Will et al., , 2015)), where specimens from the north and west coast of North Island New Zealand form a cluster, with a genetic break near East Cape, and a separate Chatham Island cluster.The abrupt change in genetic markers at East Cape might be related to ocean currents in the area that reduce the opportunity for gene flow between adjacent populations.Here, the southerly moving east Auckland current becomes the east coast current and moves slightly offshore.South of East Cape, a branch of the southern current moves northwards close to the coast (Heath, 1985).A detailed sampling of C. maculosa populations found mtDNA haplotype variation over the same coastline and detected abrupt changes in haplotype frequencies at both the southern and northern limits of the range of B. colensoi (Walton et al., 2018;Dohner et al., 2018;Fleming et al., 2018).However, the mtDNA diversity within C. maculosa and H. iris over the same geographical range was an order of magnitude smaller than that seen within the B. vittatum complex, suggesting that divergence of populations in the latter might have begun much earlier.This contrasts with studies of sympatric marine snails in the Sea of Japan, where genetic divergence within unrelated lineages was almost identical (Iguchi et al., 2007).The two New Zealand whelk lineages are ecologically similar and sympatric over much of the same coastline; however, they do not share similar patterns of shell shape and sculpture variation.Thus, we have no evidence of Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 convergent evolution and cannot determine whether there is an environmental factor that selects for phenotypic differences between B. colensoi and B. vittatum.One might infer from the signature of population differentiation of neutral nuclear genetic markers that shell shape and texture variation could simply be the result of genetic drift.However, the ecotype B. colensoi shares mtDNA haplotypes with B. v. littorinoides, suggesting that selection is preventing the homogenizing effect of gene flow.
Local adaptation among animal populations has been inferred from phenotypic traits, but determining the role of natural selection in shaping geographically partitioned variation is not simple (Merilä & Hendry, 2014).Shell shape variation of snails can result from a combination of genetic drift in isolation or selection on functional difference (Conde-Padín et al., 2009) and/ or represent ecophenotypic plasticity (Iguchi et al., 2005) and historical phylogeographical effects.If variation results from divergence owing to drift, then phenotype and neutral genetic markers are expected to show similar patterns.This is not observed in Buccinulum.Phenotypic variation within a single species maintained by selection or ecophenotypic plasticity would explain the pattern we observed.The factors that distinguish the niche of these two similar whelks might be key to understanding why B. littorinoides has a distinct East Coast ecotype, but C. maculosa does not.

Figure 3 .
Figure 3. Assignment of whelks of the Buccinulum vittatum complex into morphological or genetic clusters, ordered by geographical collection location on New Zealand coastline.A, assignment probabilities of snails into two morphological clusters assigned by MCLUST based on the first principal component of shell shape variation (51.7% of variation; optimal model had two clusters).B, assignment of snails into three clusters based on the minimum spanning network inferred from mtDNA (cox1).C, assignment probabilities of snail genotypes (from 849 nuclear loci) into two clusters from STRUCTURE output.Optimal K = 2, as indicated by STRUCTURE HARVESTER.

Figure 4 .
Figure 4. Minimum spanning networks of mitochondiral DNA haplotypes from two sympatric whelks from the New Zealand coast.Colours indicate the same sampled regions (red, north; yellow, east; blue, south).Buccinulum vittatum complex (above; cox1; 443 bp; N = 92) have two distinct haplotype clusters associated with geographical regions.Cominella maculosa (below; cox1; 573 bp; N = 555) contains less sequence diversity over the same geographical range.Hash marks represent base substitutions between haplotypes, and circles are proportional to sample size.Shells photographed represent shape diversity within these samples.

Table 1 .
Sampling location, regional designation and sample size of New Zealand coastal whelks from the Buccinulum vittatum complex and Cominella maculosa used in morphometric and genetic analyses © The Biological Journal of the Linnean Society 2018, 2018, XX, 1-17

Table 2 .
Failure of shell shape to identify New Zealand whelks within the Buccinulum vittatum complex Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 Cross-validation scores from canonical discriminant analysis for shell shape variation of Buccinulum vittatum bicinctum (Chatham), B. v. vittatum (Northern), B. v. littorinoides (Southern) and B. colensoi (Eastern) specimens labelled by geography.The proportion of individuals assigned to incorrect populations is shown based on sampling location.© The Biological Journal of the Linnean Society 2018, 2018, XX, 1-17 geneTic STrucTure of nuclear markerS

Table 3 .
Genetic population differentiation is evident from pairwise F ST .values between population samples of New Zealand whelks within the Buccinulum vittatum complex based on 849 nuclear markers

Table 4 .
Ponder (1971)rentiation among New Zealand whelks within the Buccinulum vittatum complex, pairwise F ST values based on 849 nuclear markers, presented for pooled population samples for the three taxa Downloaded from https://academic.oup.com/biolinnean/advance-article-abstract/doi/10.1093/biolinnean/bly142/5139735 by Massey University user on 25 October 2018 have slightly different shell shapes compared with conspecifics, they did not show the relatively short spire and wide aperture that characterizes B. colensoi (Fig.2).In addition, C. maculosa from this region do not exhibit the differing shell surface sculpturing observed in Buccinulum.In contrast to the single morphometric cluster identified with five principal shape components, the mitochondrial DNA sequence (cox1 gene) resolved distinct partitioning within the B. vittatum complex.Individuals from the Chatham Islands (identified as B. v. bicinctum from their location) formed a distinct mitochondrial clade.Individuals from within the range of B. v. vittatum (north) formed a distinct genetic cluster, and individuals from within the range of B. v. littorinoides (south) and B. colensoi (east) together formed a cluster of haplotypes.Analysis of 849 anonymous nuclear loci (without the Chatham Island population sample) showed this same pattern, resolving two clusters[northern B. v. vittatumdistinct from samples from eastern and southern locations(B.colensoiandB.v. littorinoides)].Thus, the genetic partitioning of both mitochondrial and nuclear markers was in conflict with the recognized taxonomy.In contrast to the mtDNA data, the nuclear markers suggest that gene flow between populations from within the range of B. v. vittatum (north) and B. colensoi (east) is ongoing.This finding is in agreement with observations byPonder (1971)of rare 'intermediate specimens' on the north-east side of East Cape, where the two taxa meet.The three snails in our study from Hicks Bay, East Cape, which grouped with B. v. vittatum in mtDNA analysis, shared nuclear (SNP) alleles with both B. vittatum complex genetic clusters.