Genetic structure of the grey side-gilled sea slug (Pleurobranchaea maculata) in coastal waters of New Zealand

Pleurobranchaea maculata is a rarely studied species of the Heterobranchia found throughout the south and western Pacific–and recently recorded in Argentina–whose population genetic structure is unknown. Interest in the species was sparked in New Zealand following a series of dog deaths caused by ingestions of slugs containing high levels of the neurotoxin tetrodotoxin. Here we describe the genetic structure and demographic history of P. maculata populations from five principle locations in New Zealand based on extensive analyses of 12 microsatellite loci and the COI and CytB regions of mitochondrial DNA (mtDNA). Microsatellite data showed significant differentiation between northern and southern populations with population structure being associated with previously described regional variations in tetrodotoxin concentrations. However, mtDNA sequence data did not support such structure, revealing a star-shaped haplotype network with estimates of expansion time suggesting a population expansion in the Pleistocene era. Inclusion of publicly available mtDNA sequence sea slugs from Argentina did not alter the star-shaped network. We interpret our data as indicative of a single founding population that fragmented following geographical changes that brought about the present day north-south divide in New Zealand waters. Lack of evidence of cryptic species supports data indicating that differences in toxicity of individuals among regions are a consequence of differences in diet.


Introduction
The grey side-gilled sea slug (Pleurobranchaea maculata) is an opportunistic carnivore that feeds on invertebrates including sea anemones, marine worms and other molluscs [1] but also on algae [2]. It is native to New Zealand (NZ), southeastern Australia, China, Sri Lanka and Japan where it is found in habitats ranging from sandy sediments to rocky reefs, and from shallow sub-tidal flats to depths of 300 m [1,3]. Little is known about the life history of the species but studies of comparative development report the production of planktotrophic veligers that hatch within eight days and remain planktonic for three weeks before juveniles settle [1,3,4]. In late 2009 this otherwise little-known sea slug attracted attention after it was implicated in dog deaths on beaches in Auckland [5]. Analyses of vomit and gastrointestinal contents revealed that deaths were a consequence of tetrodotoxin (TTX) poisoning associated with ingestion of P. maculata [5]. This was the first time that TTX had been reported in NZ and in a species of the taxonomic clade Heterobranchia [5]. P. maculata that have recently invaded coastal waters of Argentina also contain TTX [6,7].
TTX is a potent neurotoxin found in numerous terrestrial and marine organisms, but neither the origin of TTX nor the causes of variation in TTX levels among species are understood. The structure of TTX suggests a microbial origin [8] and while certain microbes have been implicated in TTX production (reviewed in [9]), many of the claims have been refuted [9,10]. Nonetheless, while not excluding a microbial origin, there is recognition that TTX in animals is often acquired via diet. For example, variability in TTX levels found in puffer fish has been attributed to exposure to toxic food sources (reviewed in [8]). For P. maculata, there is mounting evidence that toxin accumulation occurs through feeding [11][12][13][14]. An alternate possibility is that TTX arises from commensal or symbiotic microorganisms that are associated with P. maculata [14], but no TTX-producing bacteria have been found [15,16]. Observation of a significant number of egg masses during the period when toxin levels peak in TTX-associated P. maculata populations [3], and vertical transfer of TTX from adults to [13] suggests that TTX may serve a defensive function.
Studies of individual and temporal differences in TTX concentration have established that P. maculata populations from northern regions of the North Island (Whangarei, Auckland, Tauranga) have high TTX concentrations (the highest average being 368.7 mg kg -1 per individual), while populations from the South Island (Nelson and Kaikoura) have trace amounts of TTX (<0.5 mg kg -1 ) or none at all [3,5,11,14]. A recent study reported TTX concentrations as high as 487 mg kg -1 [14]. Significant individual and seasonal variations have also been observed [3]. A single individual obtained from Wellington in the south of the North Island was found to have a low concentration of TTX (2.2 mg kg -1 ) supporting the notion of a geographical cline [3].
The genetic structure of P. maculata is unknown, but variation in the established differences in toxicity between northern and southern populations suggests that geographic variability in TTX concentration correlates with genetic structure-even the possibility that northern and southern populations define separate species. Here we test this hypothesis using a combination of microsatellite and mitochondrial DNA (mtDNA) sequence markers.

Sampling, DNA extraction and tetrodotoxin assay
A total of 156 samples were collected from nine regions around New Zealand between 2009 and 2013 (Fig 1 and Table A in S1 File). DNA was extracted as described in Yıldırım et al. [17]. The Tauranga (TR) population included samples from Tauranga Harbour whereas the Auckland (AKL) population included samples from Tamaki Strait and Waitemata Harbour. Some samples were from the studies of Wood et al. [3] and Khor et al. [11].
At the outset of this study there was limited knowledge of the toxicity of P. maculata individuals from Wellington (WL) as only one individual had been previously tested [3]. To obtain a better understanding, the TTX concentration of eight (of eighteen) individuals collected from WL in October 2012 was determined as described in Khor et al. [11]. The TTX assay was performed at the Cawthron Institute (Nelson) using a liquid chromatography-mass spectrophotometry method that is described in McNabb et al. [5].
Population-level analyses were performed only for five populations which are Ti Point (TP), AKL, TR, WL and Nelson (NL) due to the small sample sizes of the other locations. TP, AKL and TR, which included highly toxic individuals [11] were designated as the "northern cluster", whereas the WL and NL population, which contained slightly toxic and non-toxic individuals [3,11,12] were designated as the "southern cluster".
A 1060 bp and 1153 bp region of mitochondrial cytochrome B (CytB) and cytochrome oxidase I (COI) genes, respectively, were amplified and sequenced in all 156 P. maculata individuals. For details regarding the primer pairs and amplification see Methodology and Table C in S1 File. Geneious Pro 6.1.6 (Biomatters, New Zealand) was used to trim, assemble, align and concatenate the resulting DNA sequences.

Statistical analysis
Genetic diversity. Microsatellite genotyping data were tested for scoring errors due to stuttering, null alleles, and large allele dropout using MICRO-CHECKER v.2.2 [18]. Departures from Hardy-Weinberg equilibrium (HWE) were estimated using inbreeding coefficient F IS [19] with 9,999 permutations in GenoDive v2.0b25 [20]. Correction for multiple testing was performed using the false discovery rate (FDR) method [21]. Linkage disequilibrium (LD) between pairs of loci was estimated in FSTAT v2. 9 Patterns of differentiation were also analysed using a multivariate approach. For microsatellite data, the Manhattan distance (DM), which calculates the mean character differences between individuals, and clonal distances (DCL), which assumes a stepwise mutational model (SMM) [45], were used. For mtDNA data, a distance matrix (D SEQ ) was calculated as a standardized bp difference between every pair of individuals in Mothur v1.33.3 [48]. Statistical analyses on resulting distance matrices were done using PRIMER v6 [49] with PERMANOVA + [50]. Patterns of inter-sample distances were visualized using non-metric multi-dimensional scaling ordination (MDS) [51]. Permutational multivariate analysis of variance (PERMA-NOVA) [52,53] was used to formally test for differences in genetic structures among different locations and canonical analysis of principal coordinates (CAP) [50,54] was used to discriminate among specific populations and identify their distinctiveness, using leave-one-out allocation success. PERMDISP was used to test the null hypothesis of homogeneity of within-group dispersions among populations [55]. All permutation tests used 10,000 permutations. Multivariate analyses were used as an alternative approach because they do not require strong assumptions about an underlying genetic model of population structure. The method was also used to investigate the relevance of north-south disjunction for both microsatellite and mtDNA data.
Migration. The microsatellite data were analysed with GeneClass2 [57] to identify the first-generation migrants using the Bayesian criterion of Rannala and Mountain [58] and the L home /L max likelihood by introducing northern and southern clusters as populations, with a threshold p-value of 0.01 and a Monte-Carlo resampling algorithm [59].
Neutrality tests and demographic analyses. BOTTLENECK v1.2 [60] was used to test the possibility of recent population reduction for microsatellite data assuming SMM and twophase models (TPM) with default settings using a Wilcoxon signed rank test [60]. A possible sign of a recent bottleneck was investigated also by a mode-shift analysis [61].
We used isolation-with-migration models [62] implemented in the program IMa2, [63] to evaluate historical demographic parameters of the two main P. maculata populations (as determined by the phylogenetic reconstructions), i.e., northern and southern populations, using the microsatellite allele data (twelve loci). IMa2 uses Bayesian inference and MCMC simulations of genealogies to estimate several demographic parameters, including population size (Θ) of the extant (Θ North and Θ South ) and ancestral populations (Θ Ancestral ), the split time (t) for the branching event, and asymmetric migration rates between the extant populations (m North!South , m South!North ). A step-wise mutation model was applied. Upper bounds of prior distributions of parameter values were evaluated in several trial runs. When the highest posterior probability peaks of all parameters fell well within the prior boundaries in these test runs, we ran five IMa2 runs with prior settings chosen according to these trial runs (three with wider and two with narrower prior boundary settings), all with different random seed numbers. The runs began with a burn-in period of 10 6 steps followed by 10 8 steps where every 10 3 genealogy was sampled. We achieved adequate convergence and mixing of the Markov chains as indicated by visual inspection of trend line plots, sufficient effective sample size values and similar posterior probability distributions in the five runs. We present parameter estimates corresponding to the average highest peak of the posterior probability distribution of the five runs. The parameter estimates are scaled to the mutation rate (μ) and generation time (g) and to convert them to biologically interpretable demographic units, we calculated the population sizes (N) as Θ/4μ, split time (T) in years as tg/μ, and population migration rates per generation (2NM) as Θm/2. The mutation rate is uncertain for P. maculata microsatellites and was set to 10 −4 per generation, a commonly assumed mutations rate for microsatellites [64,65], and the generation time was set to one year. Note that the direction of gene flow for 2NM is forward in time (e.g., 2Θ North m South!North indicates the number of migrants per generation that the northern population receives from the southern population).
Deviations from neutrality and demographic changes within and across the populations were calculated with Tajima's D [66], Fu's Fs [67] and mismatch distribution analysis in ARLE-QUIN for the mtDNA COI and Cytb sequences. The null hypothesis of expansion was statistically tested with the sum of squared deviations (SSD) from the expected values [68] and Harpending's raggedness index [69]. McDonald and Kreitman's [70] neutrality test was performed pooling all P. maculata COI sequences (1153 bp) in DnaSP using P. meckeli COI sequences as an outgroup species. Fisher's exact test (two-tailed) was used to identify significant deviations from neutrality.
Population size reconstruction based on the COI and CytB sequence data using Bayesian skyline plot analyses [71] was performed in BEAST 2.4.3 [72]. Since there was no evidence of mtDNA population structure (see Results) a single population was modelled. The input was prepared in BEAUti (part of BEAST suite). We used the HKY substitution model with five gamma categories, an exponential gamma shape prior of 1.0 and a log-normal kappa of 10.0, and a strict clock model with a uniform clock rate prior (parameters estimated empirically). The substitution rate was set to 2.0%/My (the average substitution rate used in several recent studies of marine invertebrates [65,73]. In the MCMC, a chain length of 1.1 × 10 8 and a preburn-in of 10 7 , with sampling every 10 4 , were used. Results were inspected and the Bayesian skyline plot analysed and reconstructed in Tracer 1.7 [74].

Tetrodotoxin levels in P. maculata from Wellington
Previous analyses have established that northern WH, AKL, TR, and Coromandel populations have high levels of TTX, marking these populations as "toxic", while southern populations from NL and Kaikoura (KK) are recorded as containing either trace, or no TTX [3,5,11,14]. For WL populations, previous measurements existed for only one individual documented as having a low level of TTX (2.2 mg/kg) [3]. For this study, 18 slugs were obtained from WL of which eight randomly chosen individuals were subject to TTX assay. Three contained extremely low concentrations (0.12, 0.16 and 0.5 mg/kg) of TTX. No TTX was detected in the remaining five individuals. Accordingly, the WL, NL and KK samples (the southern cluster) were classified as "non-toxic".

Genetic diversity
Microsatellite analyses. All loci were highly polymorphic with between five and 23 alleles for each locus (diversity statistics in Table 1 and Table D in S1 File). H e across populations ranged from 0.407 to 0.843, with an average of 0.655. Rarefaction curves for A R across each locus levelled off for each sampling location indicating that a reasonable portion of the available allelic diversity was sampled at each location (Figure A and in Table E in S1 File). Populations did not exhibit significant differences in genetic diversity for either A R (P!0.370) or H e (P!0.504). No significant linkage disequilibrium was found after Bonferroni correction (P>0.05) (

Genetic structure
Microsatellite analyses. Global genetic differentiation, estimated using F ST and Dest, was low to moderate: 0.064 and 0.122, respectively, and highly significant (P 0.0001). Differentiation was significant for most loci (Table 1). Pairwise comparisons with all estimators (Genic: χ2 = infinity, d.f. = 24; F ST = 0.074-0.122; and Dest = 0.153-0.246) showed that southern populations (WL and NL) were significantly differentiated from northern populations (TP, AKL and TR) (P<0.0001, Table 3 and Table G in S1 File). Most loci supported this pairwise differentiation pattern among populations (Dest values in Table H in S1 File). None of the estimators suggested significant differentiation among the northern populations of TP, AKL and TR. A weak but significant differentiation between the WL and NL populations (P<0.046) was found based on some estimators (Table G in S1 File). Analysis of statistical power by POWSIM showed a 100% probability of detecting population differentiation at an F ST value of 0.01. The probability of α error was~0.05 (P = 0.04 and 0.057 based on chi-square and Fisher methods, respectively), suggesting a low risk for Type I error. The differentiation pattern was therefore considered real and reinforced by F ST values between the significantly differentiated populations being for the most part greater than 0.01.
A north-south differentiation is evident from MDS ordinations of the DM (Fig 2A) and DCL (not shown). Note that although stress is relatively high (0.23) for the 2-dimensional MDS ordinations, the north-south distinction was also very clearly apparent in the corresponding 3-dimensional MDS ordinations (not shown here), which had lower stress (0.18). Additionally, PERMDISP analysis showed no significant differences in dispersion for either  Figure B in S1 File) the two-group north-south split (P<0.001), with the CAP model showing a leave-one-out allocation success of 97.26% (142/ 146 for DCL and also for DM), while there was no discrimination among the three northern (TP, AKL and TR) and two southern populations (WL and NL) (CAP, P>0.76 for all matrices), justifying their pooling into a single group. Bayesian clustering of individuals based on allele frequencies as implemented by STRUCTURE showed a ΔK value and mean log probabilities of data (LnP (x/K)) that were maximal at K = 2 (Fig 2B), again supporting the same two distinct north-south clusters (Fig 2C and 2D). This finding was not affected by inclusion of sampling locations as priors (data not shown). Mitochondrial DNA analyses. The sequence data for 156 individuals obtained from COI and CytB were submitted to GenBank (accession numbers: KY094153-KY094309 for COI and KY094310-KY094466 for CytB). MJN analysis of CytB (Fig 3) and COI (Figure C in S1 File) sequences resulted in highly similar patterns. Thus, for simplicity, the CytB network was used to infer evolutionary relationships among individuals. The network shows a lack of noticeable geographical structure; common haplotypes are shared across populations. The two most common CytB haplotypes are shared by 25 (16.0% of the total dataset) and 22 (14.1% of the total dataset) individuals (frequencies in Table I in S1 File). The network is complex and reticulated with a star-like topology; many private haplotypes descend from central shared  Figure C in S1 File). Samples from Argentina are closely related to NZ samples with just a single base pair distance from a commonly shared haplotype. The index of substitution saturation (ISS) was used to test for homoplasy due to multiple substitutions [42]. For both symmetrical and asymmetrical tree topology models and for both genes, the observed ISS values are significantly larger than the critical ISS (ISS.c) values (Table J in S1 File), which indicate that the paired partitions are not saturated, and the degree of homoplasy is low.
Analysis based on F ST values showed a moderate but highly significant global genetic differentiation for both COI (F ST = 0.022, P = 0.0001) and CytB (F ST = 0.057, P = 0.0001). However  Table 3). Differences in the distribution of allele frequencies were observable for both CytB (Fig 4A) and COI (Figure D in S1 File) and pairwise Θ ST values were low for all the comparisons for both COI and CytB (0.007-0.029); a weak significant differentiation was observed only between TR and NL for CytB (Θ ST = 0.026, P = 0.0226, Table 3). MDS ordination of the D SEQ matrices revealed no observable structure for either gene (Fig 4B). However, PERMANOVA analysis suggests that sampling location has a significant effect on the population structure for COI (F 4, 145 = 1.8791, P = 0.0173), but not for CytB (F 4, 145 = 1.3386, P = 0.1806), which was consistent with the results of F ST analysis. Pairwise PERMANOVA tests showed significant but weak differentiation between the AKL and TP (P = 0.0276), and AKL and TR populations (P = 0.0218) (Table G in S1 File). PERMDISP revealed no significant differences in dispersion among populations for either COI (F 4,141 = 0.3793, P = 0.847) or CytB (F 4, 141 = 0.358, P = 0.833). The CAP results were consistent with the structures revealed by PERMANOVA: among the five populations, AKL is weakly distant from TP and TR based on COI data. PERMANOVA analysis was repeated by contrasting the TP, AKL and TR populations (the northern population cluster) against WL and NL (the southern population cluster) to test whether the north-south differentiation identified with microsatellite data could be observed with the mtDNA data. The results indicate no significant differentiation between the clusters (COI: F 1, 145 = 0.819, P = 0.5268; CytB: F 1, 145 = 1.2325, P = 0.2982, Figure E in S1 File) and do not support the north-south differentiation identified by microsatellite data.
A phylogenetic analysis of P. maculata shows that P. maculata samples from NZ form a single clade ( Figure F in S1 File). Inclusion of the sequences from Argentina did not change the topology of the tree. Phylogenetic relationships between P. maculata individuals were

Migration and demographic changes
Microsatellite analyses. Migration analysis with GeneClass2 detected four first generation migrants (P = 0.01): two individuals sampled from the northern cluster (TR) were migrants from the southern cluster, and two individuals sampled from the southern cluster (WL) are migrants from the northern cluster (Table K in S1 File). Individuals from each cluster were more likely to belong to populations from the same cluster. However, presence of some first generation migrants between the clusters shows that these clusters are genetically connected. This was compatible with the results of CAP and STRUCTURE showing low connectivity between the clusters and admixed individuals in both northern and southern clusters (Fig 2C and 2D). The highest misclassifications between the clusters detected by CAP analysis and highest admixture proportions detected by STRUCTURE were noted in TR and WL. This suggests that the TR and WL populations could be bridges between the northern and southern clusters. Admixture in the WL population may also explain weak differentiation between the WL and NL populations that was found with F ST and PERMANOVA tests.
The Wilcoxon test did not detect recent bottlenecks in any population under either TPM or SMM models (Table 4). In addition, analysis of mode-shift in the distribution of allele frequencies suggests that all the populations exhibit a normal L-shaped pattern indicating no modeshift in the frequency distribution of alleles. Taken together these data suggest that none of the populations experienced a recent or sudden bottleneck.
Isolation-by-migration models suggest that the northern and southern clusters split approx. 8.0 kya. The population sizes were estimated as 8.2 k in north, 11.4 k in south and 38.1 k in the ancestral population. The migration rate was slightly higher in the south (3.2 migrants per generation) than in the north (1.2 migrants per generation). Note, however, that the estimated 95% highest posterior density (HPD) intervals were wide for all parameters (e.g. the 95% HPD  (Table 4) suggesting a recent population expansion or purifying selection. This was also suggested by the uni-modal mismatch distributions of pairwise base pair differences for COI and CytB haplotypes (Fig 4C and Figure G in S1 File, respectively), non-significant SSD and raggedness patterns (Table 4). Furthermore, the McDonald-Kreitman test found no evidence of positive selection: the ratio of nonsynonymous to synonymous substitutions within P. maculata (Pn/Ps = 5/136) and between species (Dn/Ds = 5/100) was statistically similar (neutrality index = 0.400, P = 0.1104). Reconstruction of the demographic history with Bayesian skyline plots showed signs of an expansion at around 5 kya for COI and CytB, but the 95% highest posterior density intervals are wide (Fig 5).

Discussion
This study marks the first attempt to describe and account for patterns of genetic diversity in P. maculata. Overall, we observed marked signals of population structure; however, the population structure suggested by microsatellite versus mtDNA data differed. The nuclear data exhibit patterns of diversity indicative of a north-south disjunction. The northern samples formed one group and southern (WL and NL) samples formed another (with few examples of migration). However, this disjunction was not supported by F ST analysis of mtDNA data, which indicated divergence among all populations.
Discordance between results obtained from nuclear versus mitochondrial markers is not uncommon, with explanations ranging from variation in lineage sorting to differences in rates of nuclear versus mitochondrial evolution. However, data from microsatellite markers that come from multiple unlinked nuclear loci are expected to provide a more accurate representation of population structure that mtDNA that are linked at a single locus [76]. Taken together, we interpret our data as indicative of a single founding population that subsequently became fragmented following geographical and oceanographical changes that led to the present northsouth divide in NZ waters.
According to this scenario, P. maculata inhabited NZ waters before the end of the last glacial maximum (LGM) when sea levels were low and NZ was a single land mass [77,78]. At some later time, most likely following the LGM (~22,000 years ago), large benthic habitats became available [79,80], and this may have facilitated population expansion possibly from north to south and fragmentation aided by warming temperatures and rising sea levels [81].
In support of this model is the haplotype network based on mtDNA data showing star-like structure and high haplotype diversity, indicative of a population expansion arising from a small initial population [75]. Additional support comes from the unimodal mismatch distribution pattern of pairwise base pair differences for COI and CytB haplotypes (Fig 4C and Figure G in S1 File, respectively), non-significant SSD and raggedness patterns (Table 4), significant negative values for Tajima's D and Fu's Fs, and Bayesian Skyline plots (Fig 5). The approximate date of population expansion was estimated to be 5 kya. This estimate assumes a 2% divergence rate for the COI gene, which while an estimate, nonetheless suggests a recent expansion that could reasonably have followed the cyclic climatic oscillations defined by the late Pleistocene era (~110-15 kya) [82]. Glaciation has been suggested as the possible cause of demographic changes in other NZ marine organisms, including triplefin species [83], whelk species Cominella virgata and C. maculosa [84] and red alga Bostrychia intricata [85].
This scenario is also supported by the analysis of more rapidly evolving microsatellites, which are useful for uncovering recent barriers to gene flow [86][87][88]. As sea levels rose at the end of the last glacial cycle (~13,000 years ago; [~13,000 years ago; 78]), geographical factors and associated oceanography established barriers to gene flow for marine organisms. In particular, confluence of the East Cape current with the Wairapapa Eddy off the east coast of the North Island (between 37-39˚S) created an oceanographic barrier [89]. A barrier was also formed by waters separating North and South Islands (the Cook Strait) [89,90]. The split time (ca. 8 kya) estimated for the north and south P. maculata populations suggests that the genetic differentiation between the populations might have happened recently following creation of the present coastline after the last glacial maximum [77,89]. The north-south disjunction identified for P. maculata can be explained by oceanographic barriers specific to NZ, but also with an isolation-by-distance model. The geographical gap between the sampling locations made it difficult to draw firm conclusions as to the origin of the disjunction.
A north-south genetic differentiation has been observed in other marine organisms from NZ (reviewed in [91,92]). Confluence of the East Cape current with the Wairarapa Eddy is thought to be responsible for population differentiation in organisms such as the amphipods Paracorophium excavatum and P. lucasi [89] and the gastropod Diloma subrostrata [93]. The Cook Strait barrier is thought to have underpinned north-south differentiation in organisms such as the green shell mussel (Perna canaliculus) [94], blackfoot paua (Haliotis iris) [95] and the Ornate limpet (Cellana ornate) [96]. It may also explain the weak differentiation between WL and NL with connectivity between them being explained by the D'Urville Current that flows from the west into Cook Strait [97].
One additional factor that has likely promoted recent population subdivision in P. maculata is the distribution of invasive species [98] that constitute a food source for P. maculata. The Asian date mussel, Arcuatula senhousia has been established in the Auckland region since the 1970s, forming large transient beds in sub-and inter-tidal areas of the Hauraki Gulf, Manukau Harbour and Whangarei Harbour [99]. Expansion of A. senhousia beds in the Hauraki Gulf appears to have preceded increases in the density of P. maculata populations. Interestingly, subsequent decline of near-shore beds of A. senhousia post 2010 was followed by rapid decline in the density of P. maculata populations [100]. Further evidence that range expansion of P. maculata may have been facilitated by availability of prey species comes from off-shore mussel farms in Tasman Bay (Nelson, NZ), where culture of the green shell and blue mussels have created new habitats for P. maculata, with high-density populations being found beneath mussel farms [100]. Recently, P. maculata was identified in Argentinean waters with the species rapidly spreading along the Atlantic coast [6,7]. The minor difference between mtDNA sequence in slugs from Argentina versus NZ raises the possibility that NZ maybe the source of the recently discovered population in Argentina.
Life history traits such as the nature of egg and larval stages are of understandable importance in shaping population structure of the species. Species with benthic eggs tend to have more structured populations than ones with pelagic eggs [101,102] and an inverse relationship between pelagic larval duration (PLD) and genetic structure has been found [92,103]: in a comparative analysis of NZ pelagic marine species, Ross et al. [92] showed a significant negative correlation between PLD and genetic differentiation. However, in the same meta-study, when NZ-wide sampling regimes were considered, NZ organisms with PLD durations similar to P. maculata (2-4 weeks) exhibit structural patterns of diverse types ranging from no structure, to a north-south disjunction, IBD and differentiation within and between sampling locations [92]. The three-week pelagic larval stage of P. maculata [3,4,91] likely confers a high dispersal capacity on the species. Migration analysis identified first generation migrants between the two clusters. However, the north-south disjunction still shows that dispersal is limited. Beyond barriers formed via ocean currents, density-dependent processes acting at regional scales may act to limit invasion by new types [104].
Our prediction that the previously recorded cline in TTX might be explained by genetic structure holds only for microsatellite markers. Had this also held for mtDNA markers a case may have been made that P. maculata is a complex of two cryptic species, but no such evidence exists, at least for the samples included into this study. Our phylogenetic analysis indicates that all P. maculata populations-including samples from Argentina waters [6]-are conspecific. Short branches with no or low bootstrap support is also indicative of lack of genetic differentiation among P. maculata. Similar lack of differentiation between toxic and non-toxic populations has been shown for Taricha granulosa newts from various localities in western North America [105,106] and the red-spotted newt Notophthalmus viridescens [107].
Having called into question substantive genetic differences between north and south populations, differences in TTX levels are thus likely attributable to exogenous factors, such as differences in associated bacteria, exposure, or diet. Work to date is strongly suggestive of diet as the major source of TTX, with P. maculata accumulating TTX via feeding [11], while offspring from TTX positive individuals raised in a TTX-free environment become free of TTX [13]. Recent work studying cultured bacteria from P. maculata found no evidence for a bacterial origin of the toxin [16], but TTX has been reported in certain prey, including a Platyhelminthes Stylochoplana species that co-occurs with TTX-containing P. maculata [12].