Concordant geographic and genetic structure revealed by genotyping‐by‐sequencing in a New Zealand marine isopod

Abstract Population genetic structure in the marine environment can be influenced by life‐history traits such as developmental mode (biphasic, with distinct adult and larval morphology, and direct development, in which larvae resemble adults) or habitat specificity, as well as geography and selection. Developmental mode is thought to significantly influence dispersal, with direct developers expected to have much lower dispersal potential. However, this prediction can be complicated by the presence of geophysical barriers to dispersal. In this study, we use a panel of 8,020 SNPs to investigate population structure and biogeography over multiple spatial scales for a direct‐developing species, the New Zealand endemic marine isopod Isocladus armatus. Because our sampling range is intersected by two well‐known biogeographic barriers (the East Cape and the Cook Strait), our study provides an opportunity to understand how such barriers influence dispersal in direct developers. On a small spatial scale (20 km), gene flow between locations is extremely high, suggestive of an island model of migration. However, over larger spatial scales (600 km), populations exhibit a clear pattern of isolation‐by‐distance. Our results indicate that I. armatus exhibits significant migration across the hypothesized barriers and suggest that large‐scale ocean currents associated with these locations do not present a barrier to dispersal. Interestingly, we find evidence of a north‐south population genetic break occurring between Māhia and Wellington. While no known geophysical barrier is apparent in this area, it coincides with the location of a proposed border between bioregions. Analysis of loci under selection revealed that both isolation‐by‐distance and adaption may be contributing to the degree of population structure we have observed here. We conclude that developmental life history largely predicts dispersal in the intertidal isopod I. armatus. However, localized biogeographic processes can disrupt this expectation, and this may explain the potential meta‐population detected in the Auckland region.


| INTRODUC TI ON
A wide variety of factors act to determine genetic structure within populations of marine organisms. For example, variation in life-history traits, such as habitat specificity or dispersal ability can result in different opportunities for gene flow. Population genetic structure can also arise as a consequence of the presence of biogeographic barriers such as ocean currents, land masses, or continental shelves limiting gene flow across the barrier. Alternatively, genetic structure can result from selection acting in different populations, for example, differential adaptation in heterogeneous environments (divergent selection).
In marine organisms, dispersal ability is highly influenced by developmental mode, which can be direct or biphasic. Biphasic species generally exhibit a pelagic larval stage, during which dispersal over large distances can occur via ocean currents. In contrast, direct developers have juveniles that resemble adults. Dispersal in most direct developers is achieved via relatively small-scale mechanisms such as floating, rafting, creeping, or hopping (Winston, 2012).In these species, the stepping stone model of dispersal (where populations exchange migrants most frequently with sites that are close proximity) is the most common pattern of dispersal (Palumbi, 1994).
This leads to a frequent isolation-by-distance pattern of population genetic variation in which genetic differentiation increases with geographic distance. This geographically limited dispersal of direct-developing species predicts that they should exhibit greater population genetic structure relative to their biphasic counterparts (Ayre et al., 2009;McMillan et al., 1992;Pelc et al., 2009;Puritz et al., 2017;Waples, 1987). However, this is not true in all cases (Palumbi, 1994). For example, a comparative study on the phylogeography of Australian marine invertebrates showed no effect of a biogeographic barrier on genetic structure for two direct-developing species, the banded periwinkle, Austrolittorina unifasciata, and the carpet sea star, Meridiastra calcar (Ayre et al., 2009), while the same barrier had a strong effect in six biphasic species (Ayre et al., 2009).
Gene flow across this biogeographic barrier in these two direct developers was enabled by small patches of habitat across the barrier functioning as stepping stones for dispersal. This barrier consisted of a 300 km stretch of highly variable environmental conditions, alongside habitat and oceanographic discontinuity (Ayre et al., 2009).
Thus, it has been proposed that habitat availability may be a better predictor of genetic structure than life history (Ayre et al., 2009).
Even within restricted taxonomic divisions, population genetic structure is not easily predicted from developmental life history.
Marine isopods, like all peracarid crustaceans, are direct developers. Previous work has established that marine isopods often exhibit strong genetic structure over small spatial scales, on the order of tens of kilometers or less (e.g., Idotea chelipes (Jolly et al., 2003), Austridotea lacustris (McGaughran et al., 2006), and Jaera albifrons (Piertney & Carvalho, 1994)). This small-scale structure is congruent with the hypothesis of reduced dispersal in direct developers and may be responsible for the widespread occurrence of multiple cryptic species of isopods (Ligia and Tylos spp.) on the Southern Californian coastline (Hurtado et al., 2010(Hurtado et al., , 2013Markow & Pfeiler, 2010).
However, the mangrove boring isopod, Sphaeroma terebrans, is also a direct developer, but is widely distributed across both the Atlantic and Indian Oceans. The relationship between life history and population structure is complicated by the fact that multiple factors can act to influence population structuring. For example, Riginos et al. (2011) observed that, in fishes, both life-history traits (i.e., egg type) and biogeography (i.e., biogeographic regions delineated on a range of contemporary and historical factors) were significant predictors of population structure. Furthermore, the New Zealand intertidal isopod Limnoria, exhibits population genetic structuring that is associated with both kelp genetic structuring, and large-scale ocean currents in the sub-Antarctic-suggestive of a strong influence of rafting on genetic structuring (Nikula et al., 2010). Thus, the specific mechanism of dispersal (e.g., rafting) can shape population structure, alongside dispersal mode (i.e., larval dispersal vs. direct development). Thus, it is clear that developmental life history alone is not a reliable predictor of population genetic structure.
Population genetic structure can also be influenced by geography (such as isolation-by-distance) or selection (such as isolation-by-adaptation). In particular, geographic distance has the potential to influence population structure, especially for coastally restricted direct-developing species with purported limited dispersal (Hellberg, 1994;Palumbi, 2003). This can produce correlations between population genetic differentiation and geographic distance between populations (i.e., isolation-by-distance). However, selection often cannot be ruled out as a contributor to this correlation because geographic distance can also create environmental gradients over which selection can occur. For example, geographic distance over a latitudinal range can be associated with a temperature gradient, which may promote selection on certain loci. Disentangling the selective and nonselective causes of population structure can be achieved by examining neutral versus non-neutral genetic variation (Kirk & Freeland, 2011). This is because greater genetic divergence at non-neutral loci compared to neutral loci is predicted between populations experiencing divergent selection.
Additional work is required to elucidate the complex interplay between life history (dispersal potential and development mode), biogeography, and selection in determining population genetic structure. To shed light on the relative importance of these factors on population genetic structuring, we examined population structure across various spatial scales in a direct-developing marine isopod endemic to New Zealand, Isocladus armatus. Isopods can be K E Y W O R D S evolution, genetics, genomics, genotyping-by-sequencing, isolation-by-adaptation, isolationby-distance, isopod, population, radseq particularly informative for studying these questions in direct-developing organisms, as their limited dispersal potential is expected to drive divergence between populations.
New Zealand presents an excellent opportunity to study marine isopods because it is a hotspot of marine isopod biodiversity (Bruce, 2009;Hurley & Jansen, 1977). Additionally, many New Zealand isopod species are abundant and easily sampled in intertidal zones across extensive geographic ranges (Bruce, 2009;Hurley & Jansen, 1977;Wells & Dale, 2018). Here we explore possible mechanisms responsible for maintaining this diversity using the intertidal isopod, I. armatus. This species exhibits several characteristics that make it an ideal candidate for study. First, populations are found in abundance on easily sampled semisheltered rocky shorelines.
Second, they are found across a wide geographic range throughout the country (Jansen, 1971). Third, it is a highly color polymorphic species. This is interesting both from a natural history point of view, and because it is an easily identifiable trait that may be under strong selection, possibly affecting genetic diversity and population structure. Fourth, I. armatus is highly mobile, and has a strong swimming ability, often found swimming within the incoming tide (Jansen, 1971;Morton & Miller, 1973). This has the potential to either increase or decrease gene flow, as strong swimming ability may encourage dispersal, but may also help limit the tidal displacement of individuals. The mechanism of dispersal in I. armatus is not known.
Although I. armatus is known to swim more in disturbed water, they will rapidly exit the water column and settle onto seaweed when available (Jansen, 1971). This behavior could facilitate their dispersal via rafting on dislodged seaweeds in the dynamic littoral zone, especially as many Sphaeromatid isopods, including those within Isocladus, are found in macroalgal holdfasts (Hurley & Jansen, 1977).
A previous population genetic study of I. armatus found no significant population genetic structure between two sites separated by 11 km of coastline (Wells & Dale, 2018), but strong population genetic structure over a larger geographic scale (1,000 km). However, because no sampling was conducted at intermediate spatial scales, it is unclear at what scale this population genetic structure begins to break down. Furthermore, the stretch of coastline sampled in Wells and Dale (2018) is intersected by multiple biogeographic breaks such as the East Cape and the Cook Strait (Figure 1), and crosses multiple biogeographic regions. These biogeographic regions, defined by Shears et al. (2008) describe marine regions across New Zealand that exhibit discrete and identifiable differences in community assemblages and structure to each other. Recent evidence from a multispecies analysis of population genetic variation in New Zealand points to correspondence of the borders of these biogeographic regions with genetic structuring (Arranz Martinez, 2017). Thus, it is not well-established how population genetic structure in this isopod species is affected by biogeographic barriers.
In order to shed light on these factors, we used genotyping-by-sequencing (GBS) to resolve population genetic structure in I. armatus for both neutral and non-neutral loci. We sampled populations over a range of spatial scales and at locations intermediately located between previously sampled sites. These data allowed us to test (a) how population structure and dispersal is influenced by geographic distance, (b) whether known biogeographic barriers such as the East Cape and the Cook Strait influence gene flow and dispersal, and (c) whether there are divergent selective pressures among populations that are responsible for population genetic divergence.

| Sample collection
We collected specimens of I. armatus between May and July 2018, from around the North Island, New Zealand, from locations where F I G U R E 1 Map of sampling localities within New Zealand (colored dots). The prevailing ocean currents and biogeographic breaks proposed by Shears et al. (2008)  I. armatus had previously been recorded (Hurley & Jansen, 1977;iNaturalist, n.d.). These sites were Stanmore Bay, Browns Bay, Opito Bay (Coromandel), Mt Maunganui, Māhia Peninsula, and Wellington ( Figure 1). At each site we collected a minimum of 32 individuals.
Because I. armatus is extremely color polymorphic, we opted to sample evenly across morph type to allow us to test for any influence of morph type on population structure. When possible, we collected specimens larger than 5 mm in order to provide enough tissue for DNA extraction (see Table A1 for details). We ensured that the maximum distance between individuals collected at any site did not exceed 30 m. We stored samples at −80°C in 100% ethanol until DNA extraction. For all analyses, here, we included previously collected samples from 2015 from Hatfield's Beach, Stanmore Bay (sampled again in 2018 for this study), and Kaikoura (see Wells and Dale (2018) for sampling methods) to increase the number of sampled individuals and the spatial sampling range. In addition, the two sampling events of Stanmore Bay were also used to explore any short-term changes in allele frequencies within the population.

| DNA extraction
We extracted DNA following a modified Qiagen DNEasy Blood and Tissue protocol from Wells and Dale (2018). Briefly, we used 178 µl of 0.5 M EDTA (pH 8) and 22 µl of 20% SDS in each extraction (rather than varying volumes by weight of tissue). Additionally, we eluted DNA from the spin column three times, using 50 µl of nuclease-free water for each elution. We let the eluent sit on the column for 15 min before centrifugation for one minute at 7,000 rcf.

| Data collection and processing
DNA samples were processed by Diversity Arrays Technology (DArT) Ltd using DArTseq, a genotyping-by-sequencing (GBS) approach.
The restriction enzymes PstI and SphI were chosen for complexity reduction, the complete methodology for this approach is outlined in Wells and Dale (2018), but see (Kilian et al., 2012). DArT performs SNP calling using a proprietary pipeline. SNPs are only called if both homozygous and heterozygous genotypes can be identified (Wells & Dale, 2018).

We analyzed the dataset provided by Diversity Arrays
Technology together with the data from Wells and Dale (2018). To ensure the datasets were compatible, we filtered each dataset separately based on the conditions described below using the R packages dartR (Gruber et al., 2018) and radiator (Gosselin, Lamothe, Devloo-Delva, Grewe, 2020). We then used only the SNPs shared across both datasets for the remainder of the analyses.
We required SNPs to have a call rate ≥ 0.9 (i.e., a genotype was identified for at least 90% of individuals), a minor allele count of at  (Luikart et al., 2003). We identified SNPs under selection using BayeScan (Foll & Gaggiotti, 2008), with population as the grouping factor. SNPs exhibiting a q-value of ≤ 0.05 were excluded from any further analyses.
While it is common practice to filter SNPs based on being For some analyses, we separately analyzed the SNPs that were deemed under selection based on the BayeScan results. To do this, we performed the same filtering steps above, except that we did not filter on minor allele count or heterozygosity, and we removed SNPs that were considered neutral (q-value ≥ 0.05) -this was because we wanted to retain any SNPs under putative selection, and both heterozygosity and minor allele count can signal non-neutrality (Hernandez et al., 2019;Oleksyk et al., 2008).

| Data analysis
We calculated population F statistics using StAMPP (Pembleton et al., 2013). We used pairwise F st as the primary measure of population genetic differentiation (Whitlock, 2011). p-values for the pairwise comparisons were adjusted for multiple comparisons using the Benjamini and Yekutieli correction in R (Benjamini & Yekutieli, 2001).
We conducted principal component analyses (PCA) using the R package adegenet (Jombart & Ahmed, 2011). This approach apportions genetic variation between individuals without any a priori assumptions of sampling location or reliance on a specific population model. Therefore, PCA is a useful first step in describing population structure. In order to understand the correspondence between the principal components and geography, we performed a Procrustes transformation of the first two principal components using MCMCpack in R (Martin et al., 2011). Procrustes transformations scale, stretch, and rotate the PCA in order to minimize the differences between two matrices (in this case, the difference between principal components and geographic coordinates).
We also examined population genetic structure using the Bayesian clustering approach implemented in the software STRUCTURE v.2.3.4 (Falush et al., 2003). We performed this analysis with all populations using the full set of neutral loci (including the repeated samples of Stanmore Bay). This method implements a model-based clustering approach which probabilistically assigns individuals to one or more populations under an admixture model.
These clustering analyses are particularly useful as they are naïve to sampling locations (thus avoiding the conflation of sampling locations and genetic populations) but incorporate population genetic models that can help describe population structure. For these analyses, we assumed an admixture model with correlated allele frequencies. This model assumes that individuals can have shared ancestry from one or more of K genetic clusters, rather than the no admixture model which assumes no shared ancestry between any of K genetic clusters. The admixture model is thus appropriate as we see clear instances of admixture in the PCA. We ran the Markov Chain Monte Carlo simulations with 100,000 iterations and a burn-in of 50,000.
We conducted ten replicates of each run and varied K from 2 to 9.
We performed the final population inference by consolidating the results for each level of K in CLUMPP (Jakobsson & Rosenberg, 2007).
Additionally, we performed a separate STRUCTURE analysis on the Auckland populations with the implementation of the locprior model at a K of 3, in order to test for fine-grain population structure within Auckland. Due to concerns regarding the inferences made when defining K, we chose to present a range of realistic values for K (Lawson et al., 2018;Pritchard et al., 2010;Verity & Nichols, 2016).
In the case of hierarchically arranged populations, different values of K can provide different, but relevant, biological meaning (Hahn, 2019;Lawson et al., 2018;Pritchard et al., 2010). The assumption that there is a true single value of K is rarely correct and may lead to misinterpretations-especially in the absence of reliable demographic information, or where true discrete populations do not occur (Lawson et al., 2018;Meirmans, 2015). As a result of these situations, some authors have advocated for the presentation of multiple values of K and the use of penalized log likelihood for selection of K (Hubisz et al., 2009;Meirmans, 2015;Rosenberg et al., 2002).
In order to test whether the observed population genetic structure could be explained by an isolation-by-distance model, we tested for a correlation between genetic and geographic distances using standard and partial Mantel tests in the R package vegan (Oksanen et al., 2019). In a standard Mantel test we used a matrix of Slatkin's linearized F st (transformed using 1/1 − F st (Rousset, 1997)) as a measure of genetic differentiation, and an overwater distance matrix was used as an indicator of geographic distance. We calculated overwater distance using the marmap (Pante & Simon-Bouhet, 2013) and fossil (Vavrek, 2011) R packages, finding the minimum distance between populations around the coast within a depth range of 150 m.
For the partial Mantel test, we repeated the analysis above, but included a matrix of linearized pairwise non-neutral F st as a covariate.
This analysis allows us to examine whether any signature of isolation-by-distance could be due to potential selection across environmental gradients preferentially affecting loci under selection, or the result of nonadaptive differentiation generated by geographic distance itself.
We conducted a hierarchical Analysis of Molecular Variance (AMOVA) using the R packages ade4 (Dray & Dufour, 2007) and poppr (Kamvar et al., 2014). This approach partitions the variance in genetic data among and within each level of a predefined hierarchy (e.g., individuals, populations, or large geographically defined groups) (Excoffier et al., 1992). We tested the significance of the AMOVA by performing 1,000 random permutations using the R package pegas (Paradis, 2010

| RE SULTS
To examine population structure in I. armatus populations, we isolated 261 individuals distributed across eight populations across New Zealand (Figure 1). We obtained DArTseq SNP data for these 261 individuals, which identified 78,927 SNPs as being polymorphic.
All analyses showed clear evidence of population structure between most populations. We first conducted F statistics on this filtered SNP data to test for population differentiation. We found ev- Population structure was examined using principal component analysis, which identifies the combinations of SNPs that vary the most between individuals. We found that PC1 accounted for 19.5%, PC2 for 3.65%, and PC3 for 2.11% of the variance in SNP allele frequency (Figure 2a which is also seen for greater values of K ( Figure A3) (Evanno et al., 2005). Therefore, we rely on K = 4 as the most appropriate value of K in our discussion, based on penalized log likelihood ( Figure A4) (Hubisz et al., 2009).
We next tested for a correlation between geographic and genetic distance. The Procrustes transformation of PC1 and PC2 revealed high correspondence between genetic distance and geographic distance, suggestive of an isolation-by-distance (IBD) effect ( Figure 4). This was confirmed by the use of a Mantel test, which supported the hypothesis of IBD, revealing a significant positive correlation between Slatkin's linearized F st and geographic overwater

| D ISCUSS I ON
In this study, we quantified population genetic structure in a directdeveloping marine invertebrate across a wide range of spatial scales.
By employing a GBS approach with thousands of genomic markers, we aimed to increase the resolution of the spatial scale at which population structure can be detected. TA B L E 2 AMOVA analysis of both neutral (n = 8,020) and non-neutral (n = 392) SNPs shows population structure between populations and regions based on genetic variation, with regions accounting for the largest source of variation in non-neutral SNPs how this could have occurred is from rafting events, where individuals can disperse between populations on floating debris (Baratti et al., 2011;Nikula et al., 2010). While rafting against prevailing currents may seem counter-intuitive, Fraser et al. (2018) showed that wave-driven surface currents such as Stokes drift explain the counter-current dispersal of kelp rafts. I. armatus can be found on the same rocky coastlines as New Zealand kelps such as Ecklonia radiata and Durvillaea antarctica, and are known to cling to seaweeds in disturbed water (Jansen, 1971). Thus, rafting on kelp using Stokes drift could provide a mechanism of counter-current dispersal in I. armatus.
While these results suggest that isolation-by-distance is responsible for small-to medium-scale population structure in I. armatus, similar patterns of population structuring can also arise from a lack of migration-genetic drift equilibrium within populations. For example, if a population expansion or migration event occurred recently in the past, then the observed genetic structure should be interpreted as a historical estimate of population genetic processes (Bohonak, 1999).
Thus, the lack of population structure seen among the Auckland populations could be due to a migration-drift disequilibrium that has arisen from a recent range expansion, rather than extensive gene flow between populations. This possibility could be further investigated by conducting population expansion tests on each population.
The hypothesized biogeographic break produced by the East Cape Eddy (Figure 1) has been shown to affect population structure in a range of species. This includes direct developers such as the anemone, Actinia tenebrosa, and two species of amphipods (Stevens & Hogg, 2004;Veale & Lavery, 2012), as well as biphasic species with larval stages, such as the pāua, Haliotis iris  the marine gastropod, Buccinulum vittatum (Gemmell et al., 2018).  (Stevens & Hogg, 2004), which may experience reduced habitat availability, or a sessile species whose limited motility may inhibit successful rafting (Veale & Lavery, 2012).
Our study does not support the existence of genetic barriers across either the Cook Strait (Goldstien et al., 2006) or the East Cape (Stevens & Hogg, 2004) because we observe unremarkable genetic differentiation between populations either side of these proposed barriers. Instead, we found a strong north-south genetic disjunction between Māhia and Wellington in all analyses, with the largest F st differences being between these two locations, and genetic clustering of the Wellington and Kaikōura populations.
This north-south break is congruent with the placement of a proposed border between bioregions in this area (Shears et al., 2008) ( Figure 1). Shears et al. (2008) observed clear north-south differences across both macroalgal and invertebrate community assemblages, but no clear differences between bioregions either side of East Cape. Indeed, this north-south break appears to be broadly significant across marine invertebrate taxa, and for species with low dispersal potential is found directly south of Māhia (Arranz Martinez, 2017). Alternatively, this genetic barrier could be the result of environmental gradients imposing selective pressures on I. armatus, that is, isolation-by-adaptation through divergent selection (Nosil et al., 2009;Van Wyngaarden et al., 2017). This is supported by the observation that there was greater genetic variation between these regions in an AMOVA of non-neutral SNPs compared to an AMOVA of neutral SNPs only. While the removal of non-neutral SNPs should go some way to alleviating the effect of selective pressures in our analyses, genomic hitchhiking of linked SNPs may explain the large amount of variance between regions even in neutral loci Nosil et al., 2009).
If divergent selection is occurring, it could result in I. armatus forming either a species complex, or divergent lineages undergoing a speciation-with-gene-flow process . Within New Zealand, divergent lineages of the brooding brittle star, Amphipolis squamata, have been associated with strong north-south divergence, similar to our observations of I. armatus (Sponer & Roy, 2002).
Cryptic species have also been frequently observed in isopods (Hurtado et al., 2016;Leese et al., 2008;Markow & Pfeiler, 2010), and the degree of genetic divergence between the northern and southern group in I. armatus is similar to that found between other divergent lineages of isopods based on mitochondrial DNA (Leese et al., 2008). The potential for a species complex, or lineages undergoing speciation, is further supported by the observation of an individual from Browns Bay (which was excluded from all analyses) that, despite appearing morphologically similar to I. armatus, lacked 93% of SNPs that were present in other samples. While missing data is a common feature of reduced representation datasets, excessively high missingness in GBS data has been associated with divisions between species rather than populations (Tripp et al., 2017).

| CON CLUS ION
Isocladus armatus exhibits high levels of gene flow across small spatial scales. However, at distances greater than 20 km the level of population structure is consistent with the expectation of reduced dispersal in direct-developing species, and the presence of IBD.
Interestingly, the strongest genetic break we observed was between the Māhia Peninsula and Wellington, with populations forming a clear northern and southern grouping on either side of this break.
This was unexpected, as other well-known biogeographic barriersthe East Cape and the Cook Strait-appeared to have little effect on population genetic structure. Our results suggest either a strong geophysical barrier to gene flow occurs between these regions, or that I. armatus represents divergent lineages undergoing speciation.
Additional phylogeographic analysis and fine-scale sampling across this genetic break would help determine whether the genetic divergence we observe is the result of genetic barriers to gene flow (such as selection), or the effects of a geophysical barrier that prevent dispersal in this region.

ACK N OWLED G M ENTS
We would like to thank Stephen Shuster, Libby Liggins, and Vanessa Arranz for feedback on early drafts of the manuscript and providing valuable insight into previous research. We would also like to thank four anonymous reviewers for the comments and insights which significantly improved the manuscript.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.

E TH I C A L A PPROVA L
No ethics approval or permits were required to conduct this research.

DATA AVA I L A B I L I T Y S TAT E M E N T
Sequence data are available in the SRA under accession number: PRJNA643849, and associated metadata are deposited in the Ira Moana Project on GeOMe, under the expedition title "Isopod  Observed Heterozygosity < 0.5 8,020 Loci out of Hardy Weinberg Equilibrium within populations 4,552 Note: Datasets were first filtered separately and shared SNPs between the datasets were retained for the combined dataset. This was performed to minimize batch effects as a result of samples being sequenced at two different times. Additionally, the increased data output in the 2018 dataset (as a result of the inclusion of more individuals and populations) meant that almost twice as many SNPs were identified in this dataset, than in the 2015 dataset. By retaining shared SNPs the effect of this is minimized. Call rate refers to the proportion of individuals for which a genotype can be identified for the loci.

F I G U R E A 1 Principal Component
Analysis of Auckland populations, colored by morphotype. No substructuring or stratifying effect of coloration is observed F I G U R E A 2 Admixture barplots generated through STRUCTURE and CLUMPP for a value of 3 for K. This analysis was done with the locprior model on just Auckland populations and was used to identify fine grain structure. Mean log likelihood is −4.6 × 10 5 and variance is 8,324.2 F I G U R E A 3 Admixture plots based on Bayesian clustering analyses generated in STRUCTURE for K 7-9. Relatively little structure is observed, and a much higher amount of noise is observed as K increases. In all instances, Kaikoura and Wellington retain a high degree of separation from other peoples, and Māhia also retains some distinctiveness

F I G U R E A 4
Mean log likelihood for STRUCTURE analyses penalized by variance. As K increased, the penalized log likelihood increases, but stops increasing at a K of 4-thus, we have primarily relied on a K of 4 for most results