Combining genotypic and phenotypic variation in a geospatial framework to identify sources of mussels in northern New Zealand

The New Zealand green-lipped mussel aquaculture industry is largely dependent on the supply of young mussels that wash up on Ninety Mile Beach (so-called Kaitaia spat), which are collected and trucked to aquaculture farms. The locations of source populations of Kaitaia spat are unknown and this lack of knowledge represents a major problem because spat supply may be irregular. We combined genotypic (microsatellite) and phenotypic (shell geochemistry) data in a geospatial framework to determine if this new approach can help identify source populations of mussels collected from two spat-collecting and four non-spat-collecting sites further south. Genetic analyses resolved differentiated clusters (mostly three clusters), but no obvious source populations. Shell geochemistry analyses resolved six differentiated clusters, as did the combined genotypic and phenotypic data. Analyses revealed high levels of spatial and temporal variability in the geochemistry signal. Whilst we have not been able to identify the source site(s) of Kaitaia spat our analyses indicate that geospatial testing using combined genotypic and phenotypic data is a powerful approach. Next steps should employ analyses of single nucleotide polymorphism markers with shell geochemistry and in conjunction with high resolution physical oceanographic modelling to resolve the longstanding question of the origin of Kaitaia spat.

www.nature.com/scientificreports/ than half of the samples. No evidence of significant linkage disequilibrium was detected between locus pairs. Site-specific values of mean number of alleles, observed and expected heterozygosities, F IS and also A P were all very similar across sample locations (Table 1; Table S2). Allelic (genic) differentiation was observed in 9 of 15 pairwise sample comparisons after correction for multiple testing (Table 2). Overall, across all 10 loci, a significant difference in allelic frequencies was observed amongst mussels from the six sites (χ 2 = 57.5, DF = 20, P < 0.001). Genotypic differentiation was observed for only 1 of 15 pairwise comparisons, after correction for multiple testing (Table 3). Overall, across all 10 loci, a non-significant difference in genotypic frequencies was observed (χ 2 = 29.3, DF = 20, P = 0.082).
Pairwise estimates of F ST were low (0.000-0.011, mean ± SD = 0.003 ± 0.004) ( Table 3). Eight of 15 F ST values were statistically significant after correction for multiple testing.
Population genetic structure. The PCoA plot of multilocus microsatellite variation (site-specific means ± SE) revealed clear separation amongst most samples (Fig. 1). There was no pronounced separation of samples based on geography, although the two northern-most sites (Scott Point, Ahipara) were both located to the left of the plot (negative PCoA 1 values). Three sites (Tanutanu, Mitimiti, Oakura-listed north to south) were located close to one another with some degree of overlap on both the x and y axes. The Whatipu sample (mid-point of our sampling range) was differentiated from all other samples. Table 1. Descriptive measures of genetic variation in Perna canaliculus for 10 microsatellite loci and six site samples. Number of mussels (N), mean number of alleles per site (Na), observed heterozygosity (H O ), expected heterozygosity (H E ), the fixation index (F IS ) and the number of private alleles across all loci (A P ). Sites are arranged north to south. N = number of mussels analysed; Na = number of alleles per locus (mean ± SD); H O = observed heterozygosity; H E = expected heterozygosity; F IS = fixation index; A P = total number of private alleles. www.nature.com/scientificreports/ Assignment testing. Mean ± SD assignment probabilities for mussels from the six sites were 49.8% ± 18.5 ( Table S3a). The two most northerly sites of Scott Point and Ahipara (expected to be recipient sites of mussel spat washed up at Ninety Mile Beach) had two of the three lowest assignment values. Mussels of the four northernmost sites all had total assignment values of ~ 55% (mean ± SD = 56.65 ± 7.01), whereas mussels from the two most southerly sites, which are 300 km and 600 km south of the main spat collection region, had overall assignment values of 79% (Whatipu) and 53% (Oakura). Only 12 first generation (F 0 ) migrants were identified. Scott Point contributed and received the greatest number of F 0 migrants (4 and 3, respectively), whereas Tanutanu contributed no F 0 migrants and, along with Ahipara and Oakura, received only one migrant each (Table S3b).
Geochemical marker variation. Fifty shells from each site were processed and analysed. However, due to suspected metal contamination only 49 shells from Scott Point and Ahipara (Jan) and 47 for Ahipara (Feb) were used in subsequent statistical analyses.
Running the Q-DFA with a step-wise procedure revealed that six elemental ratios (B:Ca, Co:Ca, Li:Ca, Mg:Ca, Mn:Ca, Ni:Ca) explained the majority of the classification success (96.7%) of the Q-DFA (Table S4). These results were supported by the overall per cent classified figures of the Q-DFA model whereby 92-100% of shells were classified to their site of collection when using the above six ratios. The addition of further elemental ratios i.e., Ba:Ca, Cu:Ca, Sr:Ca, Ti:Ca and Zn:Ca, provided marginal to no improvement in classification success (Table S5).
Whilst the ability of the Q-DFA to correctly assign mussels to their site of collection varied from sample to sample, classification success was high (Table 4). Mussels collected from Whatipu and from Ahipara in February and March were correctly assigned to their site of collection with 100% classification success, whereas mussels from the remaining sites (Ahipara in January, Scott Point, Tanutanu, Mitimiti and Oakura) were correctly classified to their site of collection with between 92 and 98% success. Typically, misclassified shells were assigned to sites immediately geographically adjacent to the actual collection site (Table 4), the exception being Oakura with shells assigned to Scott Point and Ahipara. The Receiver Operating Characteristic (ROC) curves showed a high degree of sensitivity for this model, with values in the range 0.995-1.000 (Table 4).
The plot of the Q-DFA canonical scores revealed no overlap of the 95% confidence interval ellipses for the ratios of B, Co, Li, Mg, Mn and Ni to Ca in shells collected in January 2015 from the six sites, as well as the additional temporal samples collected from Ahipara in February and March, 2015 (Fig. 2). These results highlight how pronounced elemental variation contributes to differentiation of the spatial and temporal samples.
Ratios of elements to Ca within the shells varied significantly amongst the sites ( Fig. 3; Table S6) for multiple comparisons following Kruskal-Wallis tests. The ratio for Boron was significantly higher in shells from Scott Point (H (1,7) = 253.9, P < 0.001) than at other sites. Lithium (H (1,7) = 190.1, P < 0.001) and manganese (H (1,7) = 230.2, P < 0.001) ratios were significantly higher in shells from Whatipu than all other sites. Cobalt ratios were significantly higher in Ahipara mussel shells in January 2015 (H (1,7) = 248.4, P < 0.001) but not in February or March. Finally, the ratio for nickel was significantly higher in shells from Tanutanu (H (1,7) = 162.7, P < 0.001) than in shells from all other sites, and although magnesium was highly variable in shells collected from each site (particularly Ahipara in January), it was greatest in Mitimiti shells (H (1,7) = 259.2, P < 0.001). These results highlight www.nature.com/scientificreports/ the spatial (all six sites) and temporal (three samples from Ahipara) variation that exists in shell geochemistry of green-lipped mussels.
Geneland results. Genotypic variation in a geospatial context. In 20 independent runs Geneland resolved K = 2 (4 times), K = 3 (12 times) and K = 4 (4 times), where K is the number of clusters, based on genotypic variation of the mussels from the six sites. The top four models were K = 3 (twice) and K = 4 (twice). These results indicate that Geneland cannot resolve the K = 3 or K = 4 question, that the most likely interpretation is K = 3, but also that the mussels from the six sites do not constitute one homogenous group. Whatipu was placed into its own cluster in 16 runs, Scott Point in 14 runs, and Mitimiti in six runs. Ahipara, Oakura and Tanutanu clustered together in 18 of 20 runs. Pairwise F ST values were low and ranged from < 0.001 to 0.036 (mean ± SD = 0.018 ± 0.013, n = 64) (Table S6).
Spatial phenotypic variation in a geospatial context. In 10 independent runs Geneland resolved K = 6 in every instance. Based on shell geochemical variation the mussels from all six sites were very different and could easily be differentiated. Table 4. Classification success of a quadratic discriminant function analysis (Q-DFA) model run on elemental ratios in the shells of juvenile green-lipped mussels (P. canaliculus) collected in January 2015 from six sites along the west coast of the North Island of New Zealand. Additional mussels were collected from Ahipara in February and March 2015. Sites listed from north to south. Numbers in bold, on the diagonal, indicate the number of mussels from a collection site that were correctly classified as being from that site (the predicted collection site).  www.nature.com/scientificreports/ Temporal phenotypic variation in a geospatial context. In 10 independent runs Geneland resolved K = 3 in every instance for the monthly samples (January, February, March 2015) of mussels from Ahipara. Likewise, in 10 independent runs Geneland resolved K = 8 in every instance when all eight samples (i.e., six spatial samples collected in January 2015 plus the February and March 2015 samples from Ahipara) were analysed together. These results indicate that geochemical differences amongst shells collected from Ahipara at monthly intervals are sufficient to permit identification of month-specific samples: there is pronounced temporal variability in shell geochemical signatures. www.nature.com/scientificreports/ Combined genotypic and phenotypic variation in a geospatial context. In 10 independent runs Geneland resolved K = 6 in every instance. These results indicate that the mussels from all six sites are very different.

Discussion
Analysis of genotypic and phenotypic variation in a geospatial framework 34 may help with the identification of previously undetected clusters or the novel identification of source populations. The ability to utilise both phenotypic and genotypic information in assignment testing is important because phenotypic variation may be high even when genotypic variation is low within the same individuals. This phenomenon is often reported 35,36 and may occur because the genetic markers being used do not adequately reflect phenotypic variation or because genotypic variation and phenotypic variation are influenced by different processes 37 . Thus, given the differing levels of information they can provide, a combined genetic and geochemical approach may be particularly informative in identifying green-lipped mussel source populations. This approach has been used to delineate stocks of Mediterranean hake (Merluccius merluccius) and it was reported that microsatellite methods provided broad scale geospatial information of gene flow on evolutionary timescales, whereas geochemical methods could provide fine scale geospatial information on ecological relevant timescales 38 .
In terms of larval connectivity, new mussel recruits to any site are expected to originate from source sites to the south: given the weak mean northwestward flow of water along the coast 31 and the 28-35 day PLD 32 the natal sites could be 10s to 100s of km away. Balanced against this are the high levels of self-recruitment (in the range 40-70%) reported for green-lipped mussels 39 , and in addition factors such as offshore winds, storminess and low swell height in the onshore direction will all contribute to variation in natural spat deposition on Ninety Mile Beach 32 . Previous work in the North Island of New Zealand (including in the far north region) has been unable to resolve any population genetic structure (spatial differentiation) despite using a range of different genetic marker types 33,[39][40][41] . These studies did not employ geospatial analyses, and the apparent absence of spatial genetic differentiation suggested that it would be very hard, if not impossible, to identify source populations for the Kaitaia spat, given that all mussels in the region belonged to one homogeneous 'northern' group. The new analyses reported here resolve genetic differentiation between/amongst mussels from six sites in this northern region. Whilst different results were resolved by Geneland for the genotypic data alone, the most frequently reported was K = 3 clusters, with the sites of Ahipara, Mitimiti, Oakura and Tanutanu forming one group, and the sites of Scott Point and Whatipu each forming separate groups. There is no obvious association between a site's geographic location and its group membership. In terms of determining the origin of spat collected from the Ninety Mile Beach sites (Scott Point in the north, Ahipara in the south) these results indicate that we have not collected source mussels for the spat arriving at Scott Point, but that spat arriving at Ahipara may be derived from two neighbouring sites of Tanutanu and Mitimiti, which are 15 km and 48 km south, respectively, and also from the most distant site of Oakura, which is ~ 640 km south. Mussels from the northernmost Kaitaia spat site of Scott Point were consistently identified as being genotypically distinct, as were mussels from the geographically intermediate site of Whatipu, which suggests that both receive limited larval supply (gene flow) from the other four sites. Pairwise F ST estimates amongst the K = 3 clusters confirm this interpretation. Based on the microsatellite markers alone, these results indicate that Kaitaia spat are most probably not derived from a single source population/site, but are a mix of mussels derived from multiple different sites. The estimate of the total number of first generation migrants amongst all sites was 12 of 288 individuals analysed. Estimates of the number of first generation migrants between pairs of sites were all in the range 1 to 3. The two Kaitaia spat sites of Scott Point and Ahipara had the highest (n = 3) and the lowest (n = 1) number of first generation migrants, respectively. When grouped according to the Geneland results, the number of first generation migrants was 4 for Ahipara, Tanutanu, Mitimiti and Oakura, 3 for Scott Point, and 3 for Whatipu. There is no obvious geographic explanation for the patterns of gene flow (source populations versus destination sites) revealed by these analyses, but for the first time we have evidence that microsatellite variation can be used to accurately differentiate amongst samples collected 10s to 100s km apart within what was previously recognised as one large genetically homogenous 'northern' group 33,[39][40][41] .
Analysis of the shell geochemistry of mussels collected from the six sites revealed very high levels of elemental variation. In total, 17% (739 of 4396) of analysed geochemical values were imputed which could potentially introduce bias in the analysis. However, there was no correlation between importance of an element within the Q-DFA and the level of imputation (unpub. data). Assignment test success results were ~ 95% (range 90-100%), and are much higher than those reported for the same mussels when based on genotypic variation (mean ~ 50%, range 30-79%). High levels of assignment success based on shell geochemistry suggest that mussels from each site experience environmental conditions that result in different shell chemical fingerprints. Geneland recognised K = 6 clusters (i.e., mussels from each site were different from all other sites) in every run for the shell geochemistry data, confirming the distinct nature of the six samples. However, the PCoA plot revealed that shells from Scott Point and Ahipara (the sites at Ninety Mile Beach) exhibited the greatest variation in shell geochemistry. Thus, mussels from sites that are only a few kilometres apart can be differentiated with a very high degree of confidence. Clearly, this finding needs to be tested against samples from other sites to determine the extent of site specificity versus some sort of regional specificity. In a management sense, when trying to trace the source of Kaitaia spat, such high levels of shell geochemical variability are really only informative when mussels collected from Ninety Mile Beach cannot be differentiated from mussels collected from another site (the putative source) because absence of differentiation implies similarity of mussels from source and recipient sites. To resolve mussel sources would require high resolution subsampling of adult shell at the umbo using laser ablation ICP-MS. The composition of the umbo represents larval shell geochemistry which can be compared to the spatial variation in overall shell geochemistry shown in this study. Recently, this approach was used to determine level of larval www.nature.com/scientificreports/ spill over from mussel aquaculture to adjacent mussel beds 42 . Future studies using LA-ICP-MS analyses of shells along the northwest coast of New Zealand are planned to address this. Geneland analysis of combined genotypic and phenotypic variation revealed K = 6 clusters, with each sitespecific sample being different from all other samples. That is, the addition of the genotypic data set to the phenotypic data set did not increase the power of the analysis to detect differences amongst the clusters. At the moment, our analyses suggest that shell geochemical signal is more informative than the genetic signal in assigning mussels to site. Nonetheless, we view the combined application of genotypic and phenotypic data sets to be an important tool and a major step forward in terms of information (e.g., if phenotypic and genotypic traits are under different selection regimes) and management options 43 . This will be particularly true when the next generation of molecular markers is developed and employed (see subsequent section).Whilst we may not have been able to identify the source population(s) of Kaitaia spat our results indicate that combined genotypic and phenotypic data sets, when analysed in a geospatial framework, probably have the necessary discriminatory power to identify source populations. However, it is equally apparent that these analyses are, to a large extent, only as good as the source sites that are sampled. In the present study we focussed on six sites, two of which represent the Kaitaia spat collecting region and four of which are further south and could reasonably have been sources. Subsequent work will require more focussed spatial coverage of sample sites, and ideally this work will be a targeted approach to site selection based on, for example, physical oceanographic models 8 that can be used to help identify source populations, rather than the haphazard site selection process involved here. Subsequent work may also be informed by our results. For example, mussels from Whatipu were identified consistently as being different from all other mussels. This may be because the site's location on Manukau Harbour, a large body of seawater contained in a coastal embayment to the west of the city of Auckland, is expected to be characterised by a chemical composition that is very different from that of all other sites sampled. The results from Whatipu are helpful in that they suggest that further detailed examination of the geochemistry signals in the shells of mussels from Scott Point and Ahipara may permit the identification of environmental signatures that can be used to help identify source sites.
Six elemental ratios (Mg, B, Mn, Li, Ni, Co:Ca) provided the majority of classification success in discriminating amongst the six west coast mussel samples along the northern North Island of New Zealand. These six elemental ratios assigned mussels to their site of origin with 92-100% success: this information is valuable because it limits the number of elements needing to be assayed in future geochemical tag studies on this species along this coast. Reliable discrimination of mussels from sites encompassing similar distances to those in our study (10s-100s of km) have been recorded using geochemical tags for a number of marine taxa, including blue mussel (M. edulis) populations at ~ 50 km distance 44 , populations of M. californianus and M. galloprovincialis 20 km apart 20 , and New Zealand green-lipped mussel populations at a distance down to 12 km 21 . These studies highlight the power and utility of the shell geochemistry approach in identifying source populations, and also in addressing questions of spatial and temporal settlement variability at ecological relevant scales.
Populations occupying estuaries or those near streams or rivers are often easily discriminated using geochemical tags 20,45 because localised terrigenous input of elements within freshwater is important in controlling the availability of at least some elements 46 . Two of our sites, Oakura and Whatipu, are located near sources of freshwater and were readily distinguished from other sites. Although there was a tendency for a small number of mussels to be assigned to Oakura despite being 640 km to the south, mussels from the remaining sites were also able to be discriminated from one another with equal if not greater success. One explanation for this could be that the west coast of the North Island of New Zealand is an area of active terrigenous sedimentation due to the active plate boundary running through New Zealand 47 . This results in considerable uplift with erosion through a varied mix of volcanic and sedimentary rock catchments 48,49 . In the far north, the geology comprises a complex mix of allochthonous igneous rock, Miocene sedimentary rock and mixed accumulated sediments from the Pleistocene and Holocene periods 50 . Accordingly, it may well be that the elemental composition of surrounding coastal waters is complex due to being fed by riparian inputs draining a complex geological landscape.
Understanding the spatio-temporal stability of geochemical tags is vital if such tags are to be informative. Some species, such as the cockle, Cerastoderma edule, demonstrate significant temporal variability in geochemical tags among years 51 , whereas other species show seasonal stability in geochemical tags e.g., the hardshell clam, Mercenaria mercenaria 52 . Analysing geochemical tags in mussels collected from Ahipara on a monthly timeframe (January, February, March, 2015) revealed that mussels could be assigned back to their site of origin (Ahipara) with a high degree of confidence (0.995-1.000), both in the context of the other five sites and also in terms of the month of collection (i.e., spatially and temporally). There was, however, a comparatively high degree of variability in elemental ratios (e.g. Mg:Ca) amongst the mussels of the three-monthly samples. Comparing shell lengths of mussels sampled during these three months reveals a significant increase in shell length between January to February (one-way ANOVA F (2,149) = 5.931, P < 0.003), which presumably contributed to differential incorporation of elements into growing shells. The three samples from Ahipara illustrate the extent of shortterm temporal variability in elemental composition (e.g., cobalt ratios were significantly higher in Ahipara shells collected in January than in February and March) at any one site. Such variation may be highly informative in identifying spat source populations, but only as long as this shell geochemistry signal can be linked back to the source population(s).
The shell geochemistry data in particular, and the population genetics data to a lesser extent, highlight spatial variability (and temporal variability for the geochemistry) that may, in the future, be informative in identifying Kaitaia spat source populations for the New Zealand green-lipped mussel industry. Temporal variation in geochemistry likely reflects changing environmental conditions on short temporal timeframes at the sampling sites. These results present a challenge to the future identification of the source populations, given that the geochemical and genetic signals of the spat at Ninety Mile Beach will be mixed, reflecting the differing spatial and temporal contributions from different source populations. Interpreting this mixed signal will be difficult, www.nature.com/scientificreports/ but not impossible, and is akin to the use of mixed source models employed to test the origin of stable isotopes in diet analysis 25,26 . The ICP-MS method is an accurate approach to be able to discriminate populations separated by as little as 5 km along an open coast. In this study, the elemental ratios of B, Co, Li, Mg, Mn, and Ni:Ca appear to consistently discriminate mussel populations along the west coast of northern New Zealand. It is now vital that the spatio-temporal stability of these geochemical tags along this coastline be confirmed. Furthermore, if determining provenance of harvested mussels is of interest then shell geochemistry has a demonstrated role to play here as has been shown for other bivalve species e.g. Mytilus edulis 53 . The use of LA-ICP-MS methods combined with a thorough sampling regime (both temporally and spatially) will be necessary to trace green-lipped mussel spatfalls back to their site(s) of origin.
From a management perspective the identification of the Kaitaia spat source populations remains important, even as the industry moves towards commercial-scale hatchery production of spat. Hatchery production will help to secure year-round availability of spat and will permit selective breeding of spat with fitness characteristics of interest to the industry 54,55 . However, in the meantime, the supply of spat to Kaitaia that underpins the New Zealand green-lipped mussel aquaculture industry remains unpredictable, even if it is reasonably cheap to collect and distribute around the country. It seems likely that over the next several decades many mussel growers will continue with Kaitaia spat as their preferred source of juvenile mussels for stocking their farms. In addition, Aquaculture NZ has plans to expand coastal and open water aquaculture in New Zealand, and a large part of this strategy is focussed on growth of the green-lipped mussel sector 56 . Where the spat supply to meet this demand will come from is unclear, again suggesting that Kaitaia spat will remain an important resource for the greenlipped mussel aquaculture in New Zealand in the foreseeable future.
Future research to identify source populations of Kaitaia spat will include the use of more informative genetic markers (single nucleotide polymorphisms, SNPs) than the microsatellite markers employed here. SNPs are highly polymorphic co-dominant markers that can now be developed reasonably cheaply and in very large numbers (thousands to hundreds of thousands). They have proven to be useful in a range of different approaches [57][58][59] and are expected to outperform microsatellite markers for information content. This is expected to translate into far greater power to discriminate amongst samples in a spatial and temporal context. In addition, the results of the present study highlight the need for a more focussed approach to site selection. Future work will therefore be guided by newly developed high definition physical oceanographic models that will provide much more guidance for site selection than was available to us at the time of this study. Finally, Norrie et al. (2019) 22 analysed shell geochemistry in selectively bred green-lipped mussel and demonstrated a genetic component in the incorporation of elements within the shell of this species. If such population differences could be uncovered in wild mussels this could further strengthen future analyses seeking to determine site of larval origin.

Materials and methods
Collection and processing of juvenile mussels. Juvenile (sexually immature) green-lipped mussels, Perna canaliculus, 14-20 mm shell length, were collected in January 2015 from six intertidal reef sites along the west coast of the North Island of New Zealand (Fig. 4, Table S1). At these sites mussels of this size are expected to have settled within the last twelve months 60 . Mussels were cleaned of epibiota and had a periostracum from umbo to shell edge. Mussels were stored in plastic zip lock bags and frozen at − 20 °C within 1 h of collection. In the laboratory, mussels were thawed and photographed to record shell length (Table S1) via ImageJ software. Using plastic tweezers, mussels were opened and the flesh removed and placed in tubes containing 4 mL of 95% EtOH for genetic analysis. Shells were stored in plastic vials and refrozen for ICP-MS elemental analysis. For these mussel samples we consider that the genetic data (see below) relate to the site from which the larvae originated (i.e., the source population) whereas the shell microchemistry data (see below) relate to the site at which the individual mussel has resided the longest (i.e., not necessarily the source population, but if the mussel in question is a self-recruit then the site at which the mussel resides is actually the source population).
Microsatellite marker methods. DNA extraction and genotyping. Total DNA was extracted from mantle or gill tissue using Geneaid genomic DNA kits following manufacturer´s instructions. DNA concentrations and A260/A280 ratios were quantified using a NanoDrop ND-1000 (Thermo Scientific). Specimens were genotyped at 10 polymorphic microsatellite loci 61 64 . False discovery rate (FDR) control 65 was applied to P-values in all statistical analyses that included multiple comparisons. A principal coordinates analysis (PCoA) of allelic data with site as the grouping variable was carried in PRIMER using the MDS function 66 to examine variation in the multilocus data set. www.nature.com/scientificreports/ Assignment testing. To determine the most likely source of individuals at each site, self-assignment tests were performed in the software GeneClass2 67 using a Monte-Carlo resampling approach and the leave-one-out procedure, in which the individual under consideration is removed during computations for their site of collection 68 . The probability of assignment was based on 10,000 simulated individuals and an exclusion threshold of P < 0.05. Individuals that were excluded from their site of collection were assigned to another sampled site when P > 0.1. For detection of first-generation migrant individuals, we used the 'L_home' likelihood computation method to allow for the possibility that not all source sites were sampled, which is highly likely in this situation. Individuals identified as first-generation migrants were removed from the data set, and the remaining sampled individuals were used as a reference group for further assignment tests to identify the most likely source site for each first generation migrant. Where the results indicated more than one possible sampled source site (P > 0.1), the individual was assigned to the site with the highest probability. Thus, self-assignments are 'home' individuals (self recruits), whereas the first-generation migrants are 'away' individuals (a migrant from one known site to another).
Shell geochemistry. Preparation of shells. Individual shells were thawed and placed into 5 ml plastic vials and sonicated for 2 min before the water was drained. All plastic ware used to handle and store mussels after the sonication step was acid washed for 12 h in 10% HNO 3 (Sigma-Aldrich trace element grade) and rinsed three times in de-ionised water. To standardise analyses amongst sites, only the left shell of each individual was analysed. Shells were cleaned 69 21 the elemental suite analysed consisted of: 43 Ca, 66 Zn, 55 Mn, 11 B, 88 Sr, 25 Mg, 138 Ba, 63 Cu, 7 Li, 60 Ni, 47 Ti, and 59 Co. Raw elemental concentrations were converted to molar concentrations i.e. mmol mol −1 , and trace element:calcium (TE:Ca) ratios were calculated to standardise elemental concentrations.

Statistical analyses.
To assess the efficacy of geochemical markers in discriminating amongst shells from the six sites, a quadratic discriminant function analysis (Q-DFA) was performed. A Q-DFA was selected due to differing www.nature.com/scientificreports/ covariance estimates amongst samples. Due to non-normality, data were rank transformed and the Q-DFA was performed on the rank values 70 . An imputation procedure was applied when the concentration of an element was below detection limits of the instrument. This consisted of substituting the recorded 0 value with half of the minimum value of that element recorded in all shells, from all sites combined (K. Ruggiero, Dept. Statistics, University of Auckland, pers. comm. 2015). The Q-DFA was performed stepwise by firstly running the analysis on all elemental ratios and recording the overall classification success. Then the elemental ratio with the least discriminatory influence was removed and the Q-DFA run again with the difference in classification success noted. This was repeated for all elemental ratios to identify the point at which no increase in the classification success was observed and thus which combination of elements would prove to be successful in discriminating amongst shells from the different sites. Mean canonical scores were calculated for shells from each site and a classification success table generated in order to gauge the success of the Q-DFA in assigning shells back to their site of collection based on elemental ratios. Areas of Receiver Operating Characteristic (ROC) curves were calculated to provide a measure of the sensitivity of the Q-DFA by comparing true positive classifications against the false positive classifications. Kruskall-Wallis ANOVAs were used to test elemental ratios of shells amongst the sites and the temporal samples. All analyses were carried out in JMP 13.0 software 71 .
Site-specific clustering-genotypic and phenotypic variation in a spatial context. The software package Geneland v4.0.8 was employed to identify the number of clusters (K) within the six-site data set 34 .
Geneland is a spatial Bayesian analysis method that can include not only genotypic data (co-dominant microsatellite variation) but also phenotypic data (shell geochemical variation) and geographic data (latitude and longitude). It is implemented in the R language code 72 . We conducted four separate sets of analyses in Geneland: (1) the genetic data only were analysed in a geospatial setting, (2) the shell geochemistry data only were analysed in a geospatial framework, (3) the three samples from Ahipara (January, February, March 2015) were analysed in a geospatial context as a separate test of temporal variation in shell geochemistry i.e., we tested the hypothesis that shell geochemistry varied enough amongst the three samples to permit their identification as K = 3 clusters in Geneland, and we also analysed all eight samples together in a geospatial context (the six spatial samples collected in January 2015 plus the February and March 2015 samples collected from Ahipara) to test the extent of temporal variability in the range of spatial variability, and (4) the genetic and the shell geochemistry data were analysed together in a geospatial framework. This stepwise approach permits identification of clustering within the genetics dataset, within the shell geochemistry dataset, and then within the combined genotypic and phenotypic dataset. Geneland assumes approximate Hardy-Weinberg equilibrium (HWE) and linkage equilibrium between loci for the genotypic data. Because molluscan populations are rarely in HWE and are often characterised by null alleles 39,73 we accounted for this in the analysis using the 'null alleles' module. We employed the 'spatial' and 'correlated allele frequencies' modules because previous work has shown that allele frequencies in one population are similar to those in other neighbouring populations 39 . Finally, spatial uncertainty was set to zero because mussels are sessile and we know their exact site of collection. We ran multiple independent runs (100,000 iterations, thinning set at 100, minimum number of populations set to 1, maximum number of populations set to the maximum number of samples included in the analysis) to test for the number of clusters (K) and to ensure that the results converged. In all three cases (genotypic data, phenotypic data, genotypic plus phenotypic data) the results of these runs were assessed based on the frequency of the different K values and also on a ranking of the run-specific probability density values. For the genotypic data analysis we calculated the pairwise F ST values between the different clusters to help quantify the extent of differentiation amongst clusters.