Life‐history predicts past and present population connectivity in two sympatric sea stars

Abstract Life‐history traits, especially the mode and duration of larval development, are expected to strongly influence the population connectivity and phylogeography of marine species. Comparative analysis of sympatric, closely related species with differing life histories provides the opportunity to specifically investigate these mechanisms of evolution but have been equivocal in this regard. Here, we sample two sympatric sea stars across the same geographic range in temperate waters of Australia. Using a combination of mitochondrial DNA sequences, nuclear DNA sequences, and microsatellite genotypes, we show that the benthic‐developing sea star, Parvulastra exigua, has lower levels of within‐ and among‐population genetic diversity, more inferred genetic clusters, and higher levels of hierarchical and pairwise population structure than Meridiastra calcar, a species with planktonic development. While both species have populations that have diverged since the middle of the second glacial period of the Pleistocene, most P. exigua populations have origins after the last glacial maxima (LGM), whereas most M. calcar populations diverged long before the LGM. Our results indicate that phylogenetic patterns of these two species are consistent with predicted dispersal abilities; the benthic‐developing P. exigua shows a pattern of extirpation during the LGM with subsequent recolonization, whereas the planktonic‐developing M. calcar shows a pattern of persistence and isolation during the LGM with subsequent post‐Pleistocene introgression.

Studies comparing sympatric, closely related species that differ in the presence or absence of larval dispersal provide an especially informative window into the specific influence of life history on marine population connectivity by controlling for environmental factors that may also be impacting genetic structure (Dawson et al. 2014). Multiple studies have shown that the mode of larval dispersal can predict relative levels of genetic and phylogeographic structure in species comparisons (Dawson, Louie, Barlow, Jacobs, & Swift, 2002;Duffy, 1993;Haney, Dionne, Puritz, & Rand, 2008;Hellberg, 1996;Hoskin, 1997;Hunt, 1993;Lee & Boulding, 2009;Riginos & Victor, 2001;Sherman, Hunt, & Ayre, 2008). There are instances, however, where other ecological factors such as ocean currents, climate change, or ecological niche appear to override the influence of life history on genetic structure (Ayre, Minchinton, & Perrin, 2009;Iacchei et al., 2013;Kyle & Boulding, 2000;Marko, 2004;Selkoe et al., 2016).
Two intertidal sea stars found on the rocky shores of Australia have been the subject of multiple analyses of population genetic structure (Ayre et al., 2009;Hunt, 1993;Sherman et al., 2008). Meridiastra calcar and Parvulastra exigua both have extensive ranges that span most of southern and eastern Australia, with a large area of sympatric overlap where both species are abundant. While largely codistributed, these two species offer a stark contrast in modes of larval development (Byrne, 1992(Byrne, , 2006Lawson-Kerr & Anderson, 1978). Meridiastra calcar, the larger of the two species, is a gonochoric, broadcast spawning species with a lecithotrophic (nonfeeding) larval stage and a short larval duration of 1-2 weeks (Byrne, 1992(Byrne, , 2006. Parvulastra exigua lays externally fertilized, attached egg masses with benthic larval development (Byrne, 1995;Lawson-Kerr & Anderson, 1978), and P. exigua has protandry-like hermaphroditism with small individuals being predominantly male and large individuals being predominantly female with a few simultaneous hermaphrodites present (Barbosa, Klanten, Jones, & Byrne, 2012;Byrne, 1991). Some individuals self-fertilize and their embryos develop to the juvenile stage (Barbosa et al., 2012). Given these differences in reproductive and developmental modes, all else being equal, P. exigua should have higher levels of population differentiation across smaller spatial scales and lower levels of within-population genetic diversity than M. calcar, especially across parts of their respective ranges where they co-occur. As expected, a study using allozymes showed that P. exigua has lower levels of among-and within-population genetic diversity and higher levels of population subdivision than M. calcar across a spatial scale of <200 km (Hunt, 1993). A more extensive allozyme-based regional study spanning New South Wales to Tasmania reported a similar pattern (Sherman et al., 2008), as did Barbosa, Klanten, Puritz, Toonen, and Byrne (2013) using the putative mitochondrial control region between rRNAs (12S-16S), with the latter study showing population subdivision at a very small scale (50 m). In contrast, a recent study using mtDNA COI showed that M. calcar exhibited significant phylogeographic structure across a major biogeographic barrier (Bass Strait) in southeastern Australia, whereas P. exigua unexpectedly did not (Ayre et al., 2009). On the surface, these studies present differing patterns of population structure; however, the different marker types may actually be revealing different aspects of population structure.
Studies using different marker types have helped to clarify different aspects of population structure in a variety of systems, and highlight the need for a comprehensive study combining the results from multiple marker types for P. exigua and M. calcar. Multiple markers allow for more robust analyses of population structure and can be effectively used to infer critical features such as directional migration rates, effective population size, and population divergence time Lowe & Allendorf, 2010;Marko & Hart, 2011).
Here, we present results from a multilocus, comparative phylogeographic analysis of M. calcar and P. exigua to reconstruct their history.
We show that these two species differ profoundly in the relative times of population divergence across major biogeographic barriers and also in the magnitude of migration between biogeographic regions. Based on the combined data from mtDNA sequences, microsatellite loci, and nuclear intron sequences, we test how well each of these species conform to expectations based on life-history mode. We conclude that inferred gene flow and genetic diversity of these two species conform better to previous work using allozymes than to conclusions based solely on mitochondrial COI data.

| Study area
The southeastern portion of coastal Australia is divided into three biogeographic regions: Flindersian, Maugean and Peronian (Figure 1; Bennet & Pope, 1953). These regions are maintained by the contemporary interactions of the Eastern Australia Current (EAC) moving warm surface waters from north to south along the New South Wales (NSW) coastline and the Leeuwin Current transporting waters from west to east along portions of southern Australia (SA), as well as the Pleistocene land bridge between Tasmania (TAS) and mainland Australia (Reviewed in Ayre et al., 2009;Dawson, 2005;Sherman et al., 2008;Waters, 2008). We collected samples from 12 populations of P. exigua and 13 populations of M. calcar from five different geographical regions (Table 1;

| Mitochondrial DNA amplification
Following the protocols in Keever et al. (2009), we PCR-amplified and sequenced two portions of the mitochondrial (mtDNA) genome.
The tRNA/COI fragment was 703 bp for P. exigua and 714 bp for M. calcar. The putative mitochondrial control region was amplified following the protocols in Barbosa et al. (2013) and was 603 and 644 bp in each species, respectively. PCR conditions were 95°C for 3 min followed by 35 cycles of 94°C for 30 s, 57°C for 30 s, and 72°C for 60 s, and then a final 72°C extension for 5 min. For all analyses, loci were concatenated to make a single mtDNA locus of 1,350 bp (M. calcar) and 1,306 bp (P. exigua) respectively.

| Microsatellite genotyping
We genotyped samples of P. exigua at five microsatellite loci and samples of M. calcar at seven microsatellite loci, following the protocols in Keever, Sunday, Wood, Byrne, and Hart (2008), Keever et al. (2009).
Because of the work required to genotype every sample, we reduced the number of sampled localities used for genotyping by removing some regional replicate localities while increasing within-locality sample sizes (Table 1; Fig. 2).

| Next-generation sequencing
We sequenced approximately 20 individuals from 12 localities of P. exigua and 13 of M. calcar at each of five nuclear intron (nDNA) loci as described in Puritz, Addison, and Toonen (2012) and Puritz and Toonen (2013). In short, for each individual, raw sequencing reads were demultiplexed, trimmed according to base quality, and mapped to reference sequences for each locus in the program GeNeIOUS were realigned by eye. If heterozygous base pairs were found in the contig alignment, the alignment was then sorted by that base pair. The two most common haplotypes were then selected as the two allelic states for that locus.
F I G U R E 1 Map of Australia with major currents and biogeographic provinces. Note that Tasmania was only partially glaciated during the Pleistocene with glaciers restricted to high elevation plateaus

| Genetic structure analyses
MODelTeST 3.7 (Posada and Crandal 1998) selected the Tamura-Nei model as the best fit for genetic distance among mtDNA haplotypes.
We then calculated all AMOVA and pairwise measures of population structure for mtDNA in aRleQUIN 3.5 using this model, and 50,000 permutations for significance testing. We used GeNODIVe to calculate all AMOVA and pairwise measures of population structure for microsatellites using the standardized F ′ ST (AMOVA) and G ′′ ST (pairwise) statistics with 10,000 permutations for significance testing (Meirmans, 2006;Meirmans & Van Tienderen, 2004).
For P. exigua, clusters were estimated without admixture based on the results of previous studies (Hunt, 1993;Sherman et al., 2008) using nuclear genes indicating very low levels of connectivity at both regional and fine spatial scales. 99,999 steps of burn-in and 100,000 sweeps were used for 2-8 groups (the maximum range allowed by the program), with three replicates per group. Replicate runs in which K = 6 had the lowest deviance information criterion (DIC) score and results displayed are from the single lowest DIC run from the K = 6 replicates, and other replicates were consistent with the presented results. For M. calcar, clusters were estimated with the BYM admixture model (described in Durand et al., 2009), with spatial interaction optimization using a 999,999 step burn-in and 1,000,000 sweeps for 2-6 groups with three replicates per group. Although runs with K = 5 had the lowest DIC score, we selected the K = 4 run with the lowest DIC as the best model. The ΔDIC from 4 to 5 was only 150.6 for the top runs, and the fifth genetic grouping was only present in a few individuals at very low percentages across a few populations.
T A B L E 1 Locations, site codes, sample numbers, and molecular diversity of sampled populations. Diversity indices are n, number of samples; H E , expected heterozygosity; H O , observed heterozygosity; N A , number of alleles; HD, haplotype diversity, and π, average number of pairwise differences. For microsatellites, diversity measures are jackknifed across loci We also estimated genetic clusters with STRUcTURe 2.3.3 (Pritchard, Stephens, & Donnelly, 2000), using source population as a prior (Hubisz, Falush, Stephens, & Pritchard, 2009). For P. exigua, clusters, we used the "nonadmixture" mode (see above), with correlated allele frequencies, 100,000 burn-in steps, 100,000 sampling steps from 1 to 7 groups, and three replicate runs per group. STRUcTURe HaRVeSTeR 0.6.8 (Earl, 2012) was used to implement the ΔK method (Evanno, Regnaut, & Goudet, 2005), which identified K = 6 clusters as the best fit to the data. For M. calcar, clusters were estimated with admixture using correlated allele frequencies, 500,000 burn-in steps, and 500,000 sampling steps from 1 to 5 groups, using three replicate runs per group.
ΔK was highest for K = 2, but the data suggested that the relatively high ΔK for this small number of groups was artificially inflated by the relatively low negative log-likelihood values from runs with K = 1. The likelihood score approached an asymptote at K = 3 ( Figure S1), and removing K = 1 from the analysis resulted in K = 3 producing the highest the ΔK value. For these reasons, we chose K = 3 as the best-fit model for M. calcar. Results from both species were summarized across replicates using clUMPP (Jakobsson & Rosenberg, 2007) and visually displayed using DISTRUcT, grouped by population instead of individual (Rosenberg, 2003).

| Coalescent analysis
We included all six sequence loci (the exception is for locus TBP in P. exigua which was fixed across all populations) for coalescent analyses in IMa2p (Hey, 2010;Sethuraman & Hey, 2016 were combined and randomly subsampled down to 10 individuals with equal representation among localities.
20-30 preliminary runs for each species were conducted adjusting priors between runs to obtain the priors used in final analysis.
MCMCMC analyses were run for burn-in until flat trend lines appeared in all estimated parameters, and afterward 300,000 genealogies per locus were saved with 10 steps between genealogy saves.
Joint peak locations and posterior probabilities were estimated with a single "L-mode" run for each species. Initial runs of P. exigua returned abnormally large estimates of the mutation rate scalar (~24) for the mtDNA locus, likely because of the large difference in the number of haplotypes for the nDNA loci compared to mtDNA locus. For the final run, a geometric mean of all mutation rate scalars was used instead.
This scalar more closely matched the scalar estimated for M. calcar and produced results similar to using the mtDNA locus alone for estimates.

| Choosing a mutation rate for coalescent analysis
We used an estimated mutation rate for the concatenated mtDNA se-

| Genetic diversity
Across all sampled loci, M. calcar has substantially higher levels of genetic diversity than P. exigua, and higher levels of within-and amongpopulation heterozygosity, allelic and haplotype diversity, and average pairwise distances between alleles (Figure 3; Table 1). Additionally, genetic diversity of M. calcar was similar to other asterinid species that have longer-lived planktonic larvae (Keever et al., 2009;Puritz & Toonen, 2011

| Haplotype networks
The haplotype networks for both species revealed distinct regional partitioning of mtDNA genetic variation, but the details differ between the two.
For P. exigua, few haplotypes were shared between regions, and there was complete lineage sorting among SA, NSW, and WA. Alleles were shared between northern and central NSW as well as central NSW and southern NSW, but no haplotypes were present in all NSW populations ( Figure 2).
In contrast, several haplotypes from M. calcar occurred in all three NSW regions and in the TAS samples, whereas WA was completely distinct from all other regions. The haplotypes also showed a partitioning of network characteristics. For example, NSW and TAS showed more starburst patterns, with most haplotypic diversity contributed by singletons that differed by only one mutation from a far more common central haplotype (Figure 2). In contrast, many haplotypes within SA occurred at low frequencies throughout the region and tended to be several mutations different from each other.

| Bayesian population clustering
Both clustering methods produced similar results for P. exigua

| Pairwise and hierarchical genetic structure
Parvulastra exigua displayed genetic structure across all loci and levels of hierarchical sampling. Pairwise values of ϕ ST ranged from 0.018 to 1.00, with 64 of 66 comparisons significant at the p < .05 level for mtDNA loci ( Figure 5; Table S1). G ′′ ST values ranged from 0.084 to 0.982 for microsatellite loci with all 50 comparisons significant at the p < .001 level (Table S2). For nDNA loci, G ′′ ST values ranged from −0.069 to 0.990 with the majority of comparisons significant at the p < .05 level (Table S3). AMOVA for all three sets of loci confirm that the majority of the population structure in P. exigua is among populations within regions (78.36% mtDNA, 68.80% microsatellite, 71.6% nDNA; Table 2). Grouping populations by geographical region explains less variance than individual populations did, even though Bayesian clustering analyses largely confirmed grouping populations by geography.
In contrast, M. calcar has low to moderate levels of pairwise population structure, with mtDNA ϕ ST = 0.000-0.637 (Table S4), microsatellite G �� ST = 0.017 − 0.708 (Table S5), and nDNA G ′′ ST = −0.078 to 0.694 (Table S6). Results from the mtDNA data showed that populations of M. calcar exhibited more regional differentiation than P. exigua  (Table S5) with the nDNA loci indicating strong regional structure but low levels of within-region structure (Table S6). AMOVA results also show a distinct regional pattern in genetic structure with the majority of microsatellite variation in this species among individuals within populations ( Table 2). The three groupings defined by Bayesian clustering (NSW; eastern TAS; western TAS + SA) explained more molecular variance in the mtDNA loci, but explained less of the microsatellite variation.

| Coalescent analyses of migration rates, population sizes, and divergence times
Estimates of gene flow from IMa2p were consistent with the strong regional population structure in P. exigua. The highest estimated migration rate was 1.35 migrants per generation from population NSW to population SH1, with all others below one migrant per generation (Table 3). Estimates of effective population size were also relatively low (88-6,038 individuals) with the exception of the population from NSW which had an estimated population size of about 70,000 (Table 3).
In M. calcar, estimated migration rates and effective population sizes were one to five orders of magnitude higher than in P. exigua (Table 4). Migration rates between northern/central NSW and southern NSW were high in both directions (~40 migrants per generation) F I G U R E 3 Bayesian clustering results for Parvulastra exigua. Top six panels are spatial interpolations of the six genetic clusters detected by TESS. The bottom panel are the K = 6 results from STRUCTURE plotted at the population level and similarly high between eastern and western Tasmania (though only in the east to west direction). Large numbers of migrants were also being exchanged from SA to eastern Tasmanian (~60 migrants per generation), but interestingly not nearly as high between SA and western Tasmania (~4 migrants per generation). In general, migration rates were higher into Tasmania than from Tasmania, although many rates have overlapping confidence intervals. Effective population sizes estimates ranged from ~25,000 individuals in NSW to over 2 million individuals in southern NSW, although all had overlapping confidence intervals with upper bounds well above 500,000 individuals (Table 4).
Estimates of divergence times between regions differed strikingly between the two species. In P. exigua, populations from northern NSW and central NSW diverged only a few thousand years before present

| DISCUSSION
With extensive sympatric ranges that span most of southern and eastern Australia, M. calcar and P. exigua provide a robust comparative system to study the population genetic effects of their striking life-history differences. The larger of the two species, M. calcar, is a gonochoric, broadcast spawner with a relatively short (10-20 days) nonfeeding pelagic larval stage. In contrast, P. exigua has protandry-like hermaphroditism and benthic larvae, presumably with more limited potential for dispersal. Thus, all else being equal, these differences in reproductive and developmental modes make clear predictions about the degree of gene flow and scale of population genetic structure expected in comparisons of these two species, but results of such comparisons have been somewhat mixed to date (reviewed by Liggins et al., 2016;Mercier et al., 2013). Here, we use the most extensive set of loci for two cosampled marine species and show that our results support the expectation that benthic larval development strongly affects genetic and demographic characteristics of these two intertidal sea star species. The benthic developer, P. exigua, showed lower levels of genetic diversity across all markers, more inferred genetic clusters, and higher overall and pairwise population genetic structure compared to the planktonic disperser, M. calcar. Additionally, coalescent analyses

| The effects of larval dispersal on population connectivity
Results from this study, and several others on these two species (Ayre et al., 2009;Barbosa et al., 2013;Colgan, Byrne, Rickard, & Castro, 2004;Hunt, 1993;Sherman et al., 2008), have consistently demonstrated that life-history characteristics predict levels of genetic structure. We showed on several different hierarchical levels and across three different types of genetic markers that the benthic larval development of P. exigua leads to higher levels of genetic structure and consistently lower levels of genetic diversity (presumably due to smaller effective population size) when compared to M. calcar, a sympatric species with planktonic development. AMOVA and pairwise comparisons of all three marker types show consistently that the majority of genetic variation in P. exigua is partitioned among populations. This is not surprising, considering that populations of P. exigua can be strongly structured (F ST ~ 0.6) on spatial scales of tide pools only meters apart (Barbosa et al., 2013).
Meridiastra calcar also displayed relatively high levels of overall and regional genetic structure, especially compared to other asterinid species with planktonic development (Keever et al., 2009;Puritz & Toonen, 2011), but with little genetic differentiation among sites within regions, especially within NSW. Coalescent estimates of withinregion migration line up with recent models of nearshore offshore oceanographic currents (Teske, Sandoval-Castillo, van Sebille, Waters, & Beheregaray, 2015, 2016, such as the highly unbalanced gene flow seen from eastern to western Tasmania  Even with varying levels of genetic structure, M. calcar has unusually high levels of genetic diversity, especially for nuclear loci. This nDNA and microsatellite diversity, consistent with previous allozyme analyses (Hunt, 1993;Sherman et al., 2008), eclipses even the high levels of nDNA and microsatellite diversity found in Patiria miniata, an asterinid with a long-lived feeding planktonic larval stage (Keever et al., 2009), and is substantially higher than the diversity levels found in another asterinid species with planktonic lecithotrophic larval development, Cryptasterina pentagona (Puritz et al., 2012). In short, our results confirm that life-history characteristics, especially whether larval development is benthic or pelagic, and their relationship to genetic structure and diversity, are consistent with other previous comparative studies (Barbosa et al., 2013;Dawson et al., 2002;Duffy, 1993;Haney et al., 2008;Hellberg, 1996;Hoskin, 1997;Hunt, 1993;Lee & Boulding, 2009;Riginos & Victor, 2001;Selkoe & Toonen, 2011;Sherman et al., 2008).
The contrasting patterns of phylogeography inferred for these two species provide insight into their recent history and population biology, with the benthic developer showing shallower divergence times between populations than the more broadly dispersing planktonic developer. These results suggest an interesting trade-off for these life-history syndromes; despite lower rates of gene flow, the benthic hermaphrodite is actually a better colonizer, and that local extinctions are more easily replaced by new colonists on shorter timescales (Johannesson, 1988). In the planktonic disperser, the populations are colonists to arrive and persist as a population. Further, the mutational distance between low frequency haplotypes in SA for the planktonic disperser indicates that it may have once been a larger population that has recently declined, or that it was once subdivided and is now admixed (Avise, 2004).

| The effects of larval dispersal on phylogeography
The spatial phylogeographic patterns and regional genetic structuring of P. exigua and M. calcar are very similar, despite their dramatically different life histories. However, coalescent analyses estimated substantially different population divergence times for the two species.
Linking spatial genetic patterns with population divergence times provides a clearer view of the relationship between life-history characteristics and phylogeography. Climate and sea-level changes of the late Pleistocene were particularly powerful influences on the marine fauna of southern Australia (Ayre et al., 2009;Dawson, 2005). When sea levels fell to lower than 50 m below contemporary levels, a land bridge formed between Victoria and Tasmania in the present day area of Bass Strait. The emergence and subsequent submergence of the Tasmanian land bridge would have cut off circulation and gene flow between the Southern Ocean and the Tasman Sea, providing opportunities for reproductive isolation in refugia (Sinclair et al., 2016), differentiation via diversifying selection, and lineage extinctions . Several other species show strong genetic breaks or phylogenetic structure at this point (Ayre et al., 2009;Dawson, 2005;Sherman et al., 2008;Waters, 2008;Waters & Roy, 2003;York, Blacket, & Appleton, 2008 with peaks in relative glacial periods (Billups, 2004 (MPE ~ 2,000-13,000 ybp), well after the end of the last glaciation and reestablishment of contemporary sea surface temperatures and sea levels around 11,000 ybp (Bostock, Opdyke, Gagan, Kiss, & Fifield, 2006). Interestingly, these divergence estimates (with the exception of northern and central NSW) largely precede the reestablishment of the East Australian Current around 5,000 ybp (Bostock et al., 2006), but the confidence intervals of our estimate include this period as well. We hypothesize that during the last glacial maximum, most populations of P. exigua in southern Australia were extirpated leaving refugia in northern NSW and either the eastern coast of Tasmania or the southern coast of Victoria. Following the resubmergence of the Tasmanian land bridge, populations in southern NSW, SA, and TAS were then colonized by sporadic successful dispersal events, through rafting or sporadic displacement of the benthic larvae (Bryan et al., 2012;Byrne, 1995;Hart, Byrne, & Smith, 1997;Johannesson, 1988;Thiel & Gutow, 2005;Waters & Roy, 2004)  Our results with multiple genetic markers augment the recent study using mtDNA COI that showed that life-history characteristics failed to predict the effect of Bass Strait, a known phylogeographic barrier affecting the population structure of several intertidal species (Ayre et al., 2009). In this previous study, COI haplotypes of M. calcar were significantly differentiated across the Strait, whereas genetic structure in P. exigua was entirely among populations, such that grouping populations into east and west of Bass Strait in a hierarchical AMOVA explained no additional molecular variance over populations alone. Our combined mtDNA, microsatellite, and nDNA results show that virtually all genetic variation is partitioned at the population, and not regional, level in P. exigua, clarifying the findings based on COI haplotypes alone.
In addition, our coalescent analyses on the combined dataset help to explain previous findings because they show that populations east and west of Bass Strait all share a post-Pleistocene origin. In that sense, the Tasmanian land bridge was not a barrier for P. exigua. This example highlights the value of multiple marker studies incorporating coalescent analyses to further our understanding of population connectivity in marine species. Our findings indicate that phylogeographic patterns in P. exigua and M. calcar are consistent with differences in their life histories, and this result is consistent with a previous meta-analysis showing that the divergence between benthic and pelagic development was the primary correlate between life-history differences and population structure (Weersing & Toonen, 2009).

| CONCLUSION
This study shows that incorporating multiple marker types provides a more complete understanding of the contrasting local and regional population connectivity of these two sympatric species than singlemarker approaches, and together with explicit modeling of historical processes helps to makes sense of otherwise apparently contrasting results among studies. Additionally, we show that populations of P. exigua are each evolving essentially independently, with the majority of genetic diversity limited to the single population scale. Moreover, our phylogeographic results indicate that this species may be particularly sensitive to changes in sea levels, with strong subsequent effects on both patterns of extirpation and recolonization as well as the resultant population structure. Within the current context of global ocean and climate change, these results have implications for conservation of these species, and likely also for other species in the region that lack pelagic larvae.

CONFLICT OF INTEREST
The authors declare that they have no conflict of interest.