Population structure and dispersal across small and large spatial scales in a direct developing marine isopod

Marine organisms generally exhibit one of two developmental modes: biphasic, with distinct adult and larval morphology, and direct development, in which larvae resemble adults. Developmental mode is thought to significantly influence dispersal, with direct developers expected to have much lower dispersal potential. However, in contrast to our relatively good understanding of dispersal and population connectivity for biphasic species, comparatively little is known about direct developers. In this study, we use a panel of 8,020 SNPs to investigate population structure and gene flow for a direct developing species, the New Zealand endemic marine isopod Isocladus armatus. On a small spatial scale (20 kms), gene flow between locations is extremely high and suggests an island model of migration. However, over larger spatial scales (600km), populations exhibit a clear pattern of isolation-by-distance. 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. Our results indicate that I. armatus exhibits significant migration across these barriers, and suggests that ocean currents associated with these locations do not present a barrier to dispersal. Interestingly, we do find evidence of a north-south population genetic break occurring between Māhia and Wellington, two locations where there are no obvious biogeographic barriers between them. We conclude that developmental life history largely predicts dispersal in intertidal marine isopods. However, localised biogeographic processes can disrupt this expectation.

is important to carefully quantify population structure and gene flow in direct developing 56 marine species across a range of spatial scales and biogeographic contexts. This will help 57 elucidate the complex relationships between life history, biogeography, and dispersal. 58

Population Genetics in Isopoda 59
Marine isopods offer a tractable model system for investigating the dispersal potential and 60 population structure of direct-developing species. Many isopod species are abundant and 61 easily sampled in intertidal zones across extensive geographic ranges. Nevertheless, the 62 forces acting to maintain population genetic structure in marine isopods are not well 63 understood. Some species exhibit genetic structure over small spatial scales, on the order of 64 tens of kilometres or less (e.g Idotea chelipes (Jolly & Rogers, 2003), Austridotea lacustris 65 (McGaughran et al., 2006), and Jaera albifrons (Piertney & Carvalho, 1994)). This is 66 congruent with the hypothesis of reduced dispersal in direct developers, and may be 67 In contrast, other species of isopods exhibit transoceanic distributions, indicative of 71 historically high dispersal rates. One example is Sphaeroma terebrans, which is distributed 72 across both the Atlantic and Indian Oceans. The apparently high dispersal rate and wide 73 distribution of this wood-boring isopod may be result of its reliance on "rafting" for dispersal 74 (a process in which individuals use large rafts of algae or other floating debris to disperse) 75

171
To examine population structure in I. armatus populations, we isolated 261 individuals 172 distributed across 8 populations across New Zealand (Fig. 1) We first used this filtered SNP data to test for population structure using F statistics (see 177 Methods). We found evidence of genetic structure (likely due to decreased rates of genetic 178 exchange) between populations that were more than 18 km apart (Browns Bay and 179 Stanmore Bay), with Fst values ranging from 0.06 to 0.25 (Table 1). 180 Kaikōura populations from the other populations. PC2 primarily differentiated between the 190 northern populations, and revealed three potential cases of migration between Mt 191 Maunganui and the Māhia Peninsula ( Fig. 2A). Finally, PC3 differentiated the southern 192 populations, Wellington and Kaikoura. All Auckland populations clustered together, 193 suggesting these individuals form one large population. We found no difference between the 194 temporal samples from Stanmore Bay, indicating that allele frequency variation over short 195 timescales is relatively stable. 196 always clustered together across all ranges of K (Fig. 4). Again we observed no fine-grain 224 structure within Auckland when sampling location was incorporated into the analysis 225 (through use of the locprior model) (Suppl. Fig. 3). 226 As we increased K, additional population genetic structure became apparent. When

246
For all the analyses above (Fst, PCA, and STRUCTURE), there was clear evidence of 247 population structure. In addition, the Procrustes transformation showed high correspondence 248 between PCA distance and geographic distance (i.e. isolation-by-distance, IBD). In order to 249 further confirm this, we used a Mantel test. This test showed a significant positive correlation 250 between Slatkin's linearized Fst (used as a measure of population structure) and geographic 251 overwater distance (r = 0.83, P = 0.001; Fig. 5 Here we have quantified population genetic structure in a direct-developing marine 265 invertebrate across a wide range of spatial scales. By employing a GBS approach we aimed 266 to increase our power to detect even low levels of population structure. test. Loci out of HWE were first identified using the R package pegas with Bonferroni 556 correction for multiple testing, and removed using the dartR package. 557

Supplemental Analyses 558
Analysis of Molecular Variance 559 We conducted a hierarchical Analysis of Molecular Variance (AMOVA) using the R packages  We conducted this analysis using 1000 randomly sampled loci, and used the following 573 parameters from BA3-Autotune: allele-frequency mixing of 0.2125, migration mixing of 574 0.075, and inbreeding mixing of 0.075. We ran the analysis for 10 million iterations, with a 575 burn-in of 5 million iterations. 576

Supplementary Results 577
Supplementary Table 2 SNP filtering results and the amount of SNPs remaining after each filter. Datasets were first filtered separately and shared SNPs between the datasets were retained for the combined dataset. This was performed to minimise 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. Observed Heterozygosity < 0.5 8,020 Loci out of Hardy Weinberg Auckland populations, and was used to identify fine grain structure.

589
In order to understand the levels at which genetic variation could be attributed, we performed 590 an AMOVA. This analysis indicated some degree of population structure (Supplementary 591 Table 3). The majority of genetic variation was within samples (56.2%), followed by regions 592 (28.9%). Between population within region variance was just below 5%. 593

Estimation of Migration Rates 602
We tested for migration events using BayesAss3. This analysis suggested low migration 603 between most populations. These three groups correspond to 1) the Auckland based 604 populations, 2) the East Coast of North Island populations and, 3) the Southern populations, 605 similar to what was found in other analyses. These estimates indicated no migration in or out 606 of the Southern population, but some negligible migration into Wellington from Kaikōura ( Fig  607   6 and Table 3).

616
Supplementary Figure 6 Chord diagram of migration estimates between each population. 617 These estimates represent the proportion of individuals in a population that are estimated to 618 be migrants. The largest arc for each population indicates within population movement, such 619 as a self-recruitment, while the smaller arcs are indicative of between population movement 620 or migration. As a result, wider arcs indicate higher migration rates. Displayed migration 621 rates are those for which the migration rate minus the standard deviation was greater than 622 0.005.

624
High migration within the Auckland populations was observed, with estimates of migration 625 rates approximately 0.25, this value represents the proportion of individuals within a 626 population that are estimated to be migrants from the source population. Further south, 627 among the populations along the East Coast of the North Island, migration estimates were 628 consistently higher among adjacent populations than between non-neighbouring 629 populations, suggestive of a stepping stone model of migration (Kimura & Weiss, 1964). 630 These migration rates were relatively low -generally between 0.01 and 0.02 (Supplementary 631