Temporal and spatial variation in population structure among brooding sea stars in the genus Leptasterias

Abstract Temporal genetic studies of low‐dispersing organisms are rare. Marine invertebrates lacking a planktonic larval stage are expected to have lower dispersal, low gene flow, and a higher potential for local adaptation than organisms with planktonic dispersal. Leptasterias is a genus of brooding sea stars containing several cryptic species complexes. Population genetic methods were used to resolve patterns of fine‐scale population structure in central California Leptasterias species using three loci from nuclear and mitochondrial genomes. Historic samples (collected between 1897 and 1998) were compared to contemporary samples (collected between 2008 and 2014) to delineate changes in species distributions in space and time. Phylogenetic analysis of contemporary samples confirmed the presence of a bay‐localized clade and revealed the presence of an additional bay‐localized and previously undescribed clade of Leptasterias. Analysis of contemporary and historic samples indicates two clades are experiencing a constriction in their southern range limit and suggests a decrease in clade‐specific abundance at sites at which they were once prevalent. Historic sampling revealed a dramatically different distribution of diversity along the California coastline compared to contemporary sampling and illustrates the importance of temporal genetic sampling in phylogeographic studies. These samples were collected prior to significant impacts of Sea Star Wasting Disease (SSWD) and represent an in‐depth analysis of genetic structure over 117 years prior to the SSWD‐associated mass die‐off of Leptasterias.

While studies of low-dispersers commonly find high spatial genetic structure among populations (Collin, 2001;Hellberg, 1996;Hunt, 1993), studies assessing temporal genetic structure among organisms with low-dispersal life histories are less common. Je Lee and Boulding (2009) compared the genetic structure of four littorinid gastropods, two of which were planktonic-dispersers and two of which were brooders, over 10 years and four sites leading to predictions of variability based upon inferred dispersal ability. On a temporal scale, planktonic-dispersers are predicted to have high levels of genetic turnover between generations due to stochastic processes affecting larval mortality and high gene flow between populations (Eldon et al., 2016;Je Lee & Boulding, 2009). According to the sweepstakes hypothesis, relatively few individuals will contribute to the recruits of the next generation due to high variance in success at varied life-history stages, leading to high temporal genetic variability (Hedgecock, 1994;Johnson & Black, 1982). Conversely, low-dispersers with significant parental care are expected to experience low juvenile mortality and have low levels of temporal variability across generations (Je Lee & Boulding, 2009;Puritz et al., 2017).
In brooders, low gene flow from other populations also contributes to high population genetic stability over time. Je Lee and Boulding (2009) provide predictions of limited temporal genetic structure for low-dispersers over 10 years. Additional studies are needed to test this pattern in other taxa over variable lengths of time. Assessing temporal genetic stability in brooding species will provide insight into long-term environmental and demographic processes contributing to population structure and inform population responses to large-scale environmental changes.
Leptasterias is a genus of small-bodied lecithotrophic sea stars ranging from Alaska to central California composed of several cryptic species complexes. Leptasterias occur in rocky intertidal and subtidal habitats, typically measure less than 6 cm from ray tip to ray tip (Chia, 1966;Fisher, 1930;Niesen, 1973), and mature around 2 years of age (Menge, 1974). Leptasterias are lecithotrophic and females brood their young underneath their rays until the fully developed juveniles crawl away to disperse (Barreto & Bauer, 2019;Chia, 1966;Menge, 1975). Due to their brooding life history and small size, these sea stars have limited dispersal to new sites. Dispersal likely occurs by individuals rafting on macroalgae or other floating substrate; long-distance dispersal is possible though likely infrequent (Highsmith, 1985;Parker & Tunnicliffe, 1994). Due to the limited vagility of these sea stars and high susceptibility to local selection pressures like algal blooms and disease outbreaks, Leptasterias can be a coastal-indicator species reflecting local environmental health; however, proper species identification is necessary to assess changing distributions, abundances, and population health. Assessing cryptic diversity is also important for providing baseline data for monitoring ecological effects of mortality events and other environmental perturbations, especially as these types of events are predicted to increase with global climate change (Harvell et al., 2004;Jurgens et al., 2015). For example, Sea Star Wasting Disease (SSWD) is a syndrome that resulted in mass mortalities in many Pacific Coast sea star populations and is multifactorial in cause (Bates et al., 2009;Eisenlord et al., 2016;Hewson et al., 2014Hewson et al., , 2018Kohl et al., 2016;Menge et al., 2016). SSWD impacts on sea star genera, including Pisaster and Pycnopodia, were first noted in central California in 2013 (Eisenlord et al., 2016). Collections for this study were completed in 2014, before major impacts of SSWD were evident in Leptasterias (Eberl et al., 2017;Eisenlord et al., 2016;Jaffe et al., 2019;MARINe, 2015).
Many lineages within the Leptasterias genus have an unresolved taxonomic status. In several broad-scale analyses of the Leptasterias genus, cryptic lineages were identified within both L. hexactis and L. aequalis using mitochondrial molecular data (Foltz et al., 2008;Hrincevich et al., 2000). Leptasterias hexactis is comprised of two distinct clades: L. hexactis C found in Washington and L. hexactis G found in Alaska. Leptasterias aequalis is comprised of four allopatric and sympatric clades: L. aequalis B in Washington, L. aequalis A ranging from Washington to north of San Francisco Bay, L. aequalis D ranging from Washington to south of Monterey Bay, and L. aequalis K ranging from Cape Mendocino to south of Monterey Bay.
Morphological characters are challenging for identification within the Leptasterias genus due to high morphological variability within and among clades (Foltz et al., 1996) and potential hybridization (Foltz, 1997). Taxonomic uncertainty exists for L. aequalis D, which is also referred to as Leptasterias pusilla in some literature (Foltz et al., 2008). Fine-scale genetic analysis of Leptasterias will contribute to taxonomic revision and resolution in this genus.
Past studies on the broad-scale distributions of Leptasterias spp.
do not account for the fine-scale cryptic diversity within the genus and previous studies addressed the need for fine-scale analysis (Foltz, Nguyen, Nguyen, & Kiger, 2007. Indeed, a recent study used one mitochondrial locus to reveal the presence of a previously undescribed clade, Clade Y, localized around the San Francisco Bay outflow . Here, we investigate temporal and spatial population structure using nuclear and mitochondrial sequence data with widespread contemporary sampling and historic museum sampling. Previous range estimates of Leptasterias might be incorrect as museum samples were historically identified based on unreliable morphological characters; molecular identification in this study could contradict early classifications. Multilocus sequence data will be used to (a) clarify the phylogenetic relationship of California Leptasterias lineages, and (b) assess temporal and spatial patterns of population structure. We predict high population structure and high temporal stability in structure over time due to the low-dispersal potential of Leptasterias.  (Table 1). Several tube feet from each sea star were transferred from ethanol into milli-Q water and left on a shaker for 2 days to remove excess ethanol prior to extraction. Sampling was meant to be non-destructive with minimal tissue removal and whole stars were placed back into the invertebrate collection upon tube feet removal. All DNA extractions were carried out using NucleoSpin Tissue Columns (Macherey-Nagel Inc), except for samples collected in 2008 from Marshall Gulch, Bodega Bay, and Mussel Rock, which were extracted with a phenol-chloroform extraction.

| COI amplification
Primers were designed for this study from published mitochondrial cytochrome oxidase subunit I (COI) sequences of L. aequalis and L. hexactis (Foltz et al., 2008;Hrincevich et al., 2000). Forward primer COILF, 5′ GCA-GGA-TTT-ACC-CAC-TGA-TTT-C 3′ and reverse primer COILR, 5′CCT-GGC-TTC-ACA-GGC-AGA-T 3′ amplified 378 bp of COI, 68 bp of tRNA-Arg, and 90 bp of ND4L genes (henceforth referred to as COI for simplicity). PCR reactions were carried out using the same reaction concentrations and volumes as in D-Loop amplification. Thermal cycling conditions were: initial denaturation at 96°C for 120 s, 35 cycles of denaturation at 94°C for 30 s, annealing at 46°C for 30 s, extension at 72°C for 60 s, and a final extension at 72°C for 300 s. For all historic samples, thermocycling conditions were run for 40 cycles.

| Intron amplification
Five Exon Primed Intron Crossing (EPIC) loci (Chenuil et al., 2010;Gérard et al., 2013) were screened for amplification in Leptasterias based on successful amplification in other echinoderm taxa: i1, i9, i39, i43, and i51. One EPIC locus which offered the highest resolution among sites and clades was chosen and optimized for population genetic analyses. The i51 primer pair amplified a region in the gene group UDP-N-acetylglucosaminyl-transferase (Chenuil et al., 2010). Primers were redesigned for i51 to decrease primer degeneracy. Forward primer i51LF GAT-CGA-CCC-AGC-CAC-ATT and reverse primer i51LR TTG-AAG-CAA-CAG-GGG-AGA-AG were exclusively used to amplify a 277 base-pair intronic region. PCR reactions were the same as D-Loop and COI, but used 0.1 µM of each primer. Thermal cycling conditions were: initial denaturation at 96°C for 60 s, 35 cycles of denaturation at 94°C for 40 s, annealing at 45°C for 30 s, extension at 72°C for 40 s, and a final extension at 72°C for 120 s.

| Sequencing reactions
Amplification of PCR templates was assessed with gel electrophoresis using a 1.5% agarose gel stained with ethidium bromide. PCR products were cleaned using a SAP/EXO reaction following manufacturer's instructions (Affymetrix). Cycle sequencing reactions were carried out in the reverse direction using the 1/8 reaction BigDyeTerminator v3.1 (Applied Biosystems, ABI). Products were sequenced using an ABI 3130 genetic analyzer. Cloning was used to resolve and confirm a subset of alleles for i51 in heterozygous individuals using Vector System II, pGEM-T (Promega).

| Phylogenetic analysis
COI, D-Loop, and i51 sequences were edited by eye in Geneious v7.1.7 (Kearse et al., 2012) (Maddison & Maddison, 2007) to ensure the correct reading frame was used and to determine base saturation and nucleotide position changes. D-Loop and COI haplotypes were aligned to published Leptasterias sequences for phylogenetic analyses (Table S1).

Leptasterias camtschatica was chosen as an outgroup based on previ-
ous phylogenetic analysis identifying it as a sister group to L. hexactis and L. aequalis (Foltz et al., 2008). Maximum likelihood ( Rozas et al., 2003) was used to resolve the allelic phase for single nucleotide polymorphisms in i51 sequences. Allelic phase that could not be resolved with greater than 60% confidence were not used in the analysis (only one individual was excluded for less than 60% confidence). Only three individuals were assigned a phase with confidence less than 98%. The i51 phased haplotype alignment was imported into Seqstate v.1.4.1 (Müller, 2005) to code indels as simple characters (Simmons & Ochoterena, 2000) and complex characters (Müller, 2006). A phylogenetic tree was estimated in MrBayes for i51 using indels as both missing characters and coded as informative characters.
A log-likelihood value was calculated to determine whether a molecular clock was appropriate for the data using BEAST v1.8.2 , and the p-value was non-significant indicating the application of a molecular clock to be appropriate.
Divergence times were measured between concatenated historic and contemporary mtDNA haplotypes to estimate the time of differentiation between previously undescribed clades within L. aequalis in BEAST. Nuclear data were omitted due to low resolution. The TN93 + G substitution model was employed in BEAST. An uncorrelated lognormal relaxed clock with an estimated substitution rate was used with all tips set to zero. Leptasterias aequalis was constrained as monophyletic, but taxon sets within L. aequalis were not constrained as monophyletic. Foltz et al., (2008) used a molecular clock calibrated with the crossover of Leptasterias muelleri through the Bering Strait as the prior probability. The estimated divergence between L. hexactis and L. aequalis using the putative mitochondrial control region and COI was 2.56 Mya and 3.08 Mya, respectively.
Both estimates were averaged for the combined mtDNA divergence time and set as the normal distribution calibration point in this study.
Starting trees were randomly generated and the tree prior assumed a Yule speciation process. For each BEAST analysis, the MCMC was performed for 10 7 generations, sampling every 1,000 generations with a burn-in of 10%. Summary statistics were generated and visualized in Tracer v1.6.0 . Maximum clade credibility trees with median node heights of >50% posterior probabilities were calculated with TreeAnnotator v1.8.2  and were drawn in FigTree v1.4.0 (Rambaut, 2009). The analysis was run five times to confirm convergence and the combined results are reported. were constructed using indels as both informative and uninformative characters.

| RE SULTS
Historic samples amplified preferentially at COI but were often unsuccessful for amplification of D-Loop or i51 (Table S2)

| Phylogenetic analysis of contemporary and historical samples
The Genetic distances were calculated for concatenated D-Loop and COI using the TN93 + G model ( Reference sequences for L. aequalis clades D, K, A, and B and L. hexactis were obtained from GenBank ranged from 0.10%-1.30%, with a mean genetic distance of 0.70%. Leptasterias aequalis K had the highest within-clade genetic distance (1.30%). Mean inter-clade genetic distance was higher (1.97%) with a range of 0.9%-2.9%.
Bayesian reconstruction in BEAST resulted in a tree with overall similar topology to that of Figure 1. The estimated mean rate of evolution for mtDNA was 0.0135 ± 0.0001 substitutions/site/My. Divergence times were estimated between clades and time to most recent common ancestor (TMRCA) was estimated within each clade ( Figure 2). All effective sample sizes were over 1,000 with the exception of the split between L. aequalis A and Clade Z (n ESS = 949), TMRCA for Clade Y (n ESS = 517) and TMRCA for Clade Z (n ESS = 987).
The two main groupings of lineages (A and B) were estimated to have diverged between 1.74 and 1.25 Mya.
ML and Bayesian analysis of i51 alleles produced trees with the same topology. Samples obtained from Alaska populations, alleles CC, G, and SS, were used as outgroup sequences, as these samples were identified as L. hexactis through D-Loop and COI barcoding. Alleles within L. aequalis were not phylogenetically resolved ( Figure S1). The topology of the i51 phylogeny was unchanged when the 13 base-pair repeat indel was used as informative. The i51 phylogenetic tree showed low resolution and did not further resolve the relationships within lineages of the putative L. aequalis complex.  Figure S2). The haplotype map for nuclear DNA revealed patterns in which high-frequency haplotypes were abundant at central sites and present in almost all sites, revealing shared alleles between populations ( Figure S2).

| Population genetic analysis
The mtDNA minimum spanning network revealed many lowfrequency haplotypes separated by a large number of mutations ( Figure 4). Haplotypes in Clade Y showed a more typical pattern of one high-frequency haplotype with many low-frequency haplotypes separated by one or two mutations. Leptasterias aequalis K was comprised of haplotypes found both north and south of San Francisco Bay, while Clade Y was comprised almost completely of haplotypes that were bay-proximal. The i51 minimum spanning network ( Figure 4) revealed four common alleles with other low-frequency alleles separated by one or two mutations. Indels were included in the minimum spanning network and resulted in many mutations separating each cluster of alleles. When indels were excluded from the minimum spanning tree, alleles were separated by fewer mutations; however, a geographic pattern was still not evident ( Figure S3).
Haplotype diversity for all populations ranged from 0.0 to 0.95 for mtDNA and 0.10 to 0.87 for i51 (Table 3). The lowest haplotype diversity values for both mitochondrial and nuclear loci were measured at sites with predominant Clade Y abundance, whereas, the highest haplotype diversities occurred at sites with high abundances of L. aequalis K. Nucleotide diversity at all sites ranged from 0.00 to 0.025 for mtDNA and 0.0004 to 0.010 for i51 (Table 3).  (Table 3).

Mismatch distributions of mtDNA haplotypes for L. aequalis K and
Clade Y did not differ significantly from the unimodal curves expected for a sudden demographic expansion or for a rapid spatial expansion ( Figure 5). The raggedness index values were not significant for either clade and a hypothesis of sudden expansion could not be rejected (Table S5). Clade Y showed a steep peak in the distribution, consistent with a recent bottleneck or expansion event. The L. aequalis K distribution was slightly more ragged, but still consistent with a demographic expansion. While the unimodal distributions of both clades indicate recent population expansion, both expansion events and selective processes can result in a distribution of low diversity.

| Historic samples
Historic samples successfully amplified at COI, however, successful amplification of D-Loop or i51 was variable. COI haplotypes (from historic and contemporary samples) were used to build ML and Bayesian phylogenetic trees (not shown; the topology of supported branches did not differ from the concatenated mtDNA tree). Two dominant haplotypes resolved into Clades Z and Y, respectively, which were historically widespread (Figure 3), in contrast to contemporary samples, where these two dominant haplotypes were only found in bay-proximal populations (with the exception of one individual at Pigeon Point). Historical sampling revealed two clades, Clades Z and Y, as more widespread and abundant than indicated by contemporary sampling. Historic samples showed Clade Z was once found across 800 km of coastline from Crescent City, CA to Diablo Canyon, CA and Clade Y was found across 750 km of coastline from Crescent City, CA to Piedras Blancas, CA ( Figure 6). In contemporary samples, Clade Y was found at sites across 100 km of coastline with one individual found at Pigeon Point as the southern range limit.

| Divergence times and phylogenetics
COI coupled with D-Loop offered higher resolution of lineages within the Leptasterias genus than analysis using D-Loop alone  and revealed the presence of two potential species complexes historically grouped into one. Previously, L. aequalis D has been interchangeably referred to as L. pusilla (Foltz et al., 2008), implying the phylogenetic grouping found in this study of L. aequalis D, L. aequalis A, Clade Y, and Clade Z makes up a nominal species complex of L. pusilla that is separate from the L. aequalis complex.
Genetic distances between the two potential species complexes F I G U R E 2 BEAST consensus tree for Leptasterias concatenated D-Loop and COI mtDNA haplotypes constructed with Bayesian MCMC analysis . Leptasterias camtshatica was included as the outgroup sister taxon. The L. hexactis and L. aequalis split, marked with a red asterisk, was estimated by Foltz et al., (2008) and used for calibration. Reference sequences were included from GenBank (1.1%-4.3%, Hart et al., 2003). While genetic distances do not alone support the separation of clades into distinct species, additional loci and behavioral, morphological, or physiological analyses will help to resolve the relationship of lineages within the Leptasterias genus (e.g., Shaw & Cohen, 2015;Gong et al., 2019;Jaffe, 2020;Johnson et al., 2018;Johnson, 2020;Rupert, 2020).
The nuclear intron locus showed lower levels of variation than mitochondrial loci and did not resolve the putative L. aequalis into monophyletic lineages. Introns tend to accumulate mutations at a higher rate than exons and have a slower time to coalescence than mitochondrial loci (Hung et al., 2016) (Foltz et al., 2008). Following isolation in the Pleistocene, range expansion in the Holocene likely occurred, as seen in other taxa (Dawson et al., 2011;Ellingson & Krug, 2006;Hellberg et al., 2001;Jacobs et al., 2004;Marko, 1998 Figure 1. Letters represent sample site, see Table 1, (n = sample size), numbers above the circles represent the collection year, and circle size is representative of sample size driven speciation in Leptasterias, additional neutral and adaptive processes of divergence likely contributed to the maintenance of divergence and population structure.

| Patterns of population structure
A previous study by Je Lee and Boulding (2009)

| Colonization events during the San Francisco bay formation
The historically widespread and abundant Clade Y or Clade Z haplotypes could represent a source of founders colonizing the bay area.

It is possible Clade Y individuals inhabited the area that eventually
formed San Francisco Bay and colonized coastal areas around the bay following formation, resulting in their current localized distribution. Colonization events are supported by negative neutrality statistics and low haplotype diversities at Clade Y sites, as seen in other taxa experiencing genetic bottlenecks and reduced diversity associated with colonization (Marko, 1998;Hess et al., 2011) Concurrently, the unimodal mismatch curve for Clade Y suggested a recent population expansion and could reflect expansion around San Francisco Bay.

| Ocean circulation and transport processes
Leptasterias lack a planktonic dispersal stage, however, they can be considered epi-planktonic through long-distance dispersal on algal rafts (Highsmith, 1985). While these events are likely infrequent,

F I G U R E 4
Haplotype network for concatenated D-Loop and COI mtDNA haplotypes (middle) and nuclear i51 haplotypes with indels as informative characters (right when they do occur, rafters could be affected by ocean circulation processes much as planktonic dispersers are. The central CA region of coastline in this study has two potential physical barriers to dispersal for rafting organisms: San Francisco Bay and Point Reyes. Estuarine outflow from San Francisco Bay is a potential barrier to dispersal due to uninhabitable effluent conditions. Young or adults on rafts might experience mortality due to conditions associated with San Francisco Bay effluent: low salinity, high temperatures, pollutants from wastewater run-off (Conomos et al., 1985;Luoma & Cloern, 1982;Nichols et al., 1986), or even facilitated offshore transport. Puritz and Toonen (2011) found reduced genetic diversity and connectivity in the planktonic disperser Patiria miniata across areas of high human impact and pollutant run-off in the Southern California Bight attributed to larval mortality. San Francisco Bay effluent could act as a similar barrier to dispersal in Leptasterias, though additional genetic assays of other intertidal organisms around the bay outflow would help elucidate this theory.
The Point Reyes peninsula is a prominent geographic feature in the range of Leptasterias. Several studies indicate Point Reyes is a barrier to dispersal for other taxa: in the low-dispersing genera, Alderia (Ellingson & Krug, 2006), Tigriopus (Edmands, 2001), and Nucella (Marko, 1998), and in the high dispersing species, Mesocentrotus franciscanus (Moberg & Burton, 2000). Retention embayments occur both north and south of Point Reyes, which can retain nearshore waters and entrain non-local propagules (Morgan et al., 2009;Wing et al., 1998). These retention zones could effectively limit connectivity between southern and bay-proximal Leptasterias populations and populations north of Point Reyes.
The Monterey Bay region also has documented retention zones (Graham & Largier, 1997;Vander Woude et al., 2006) which could facilitate connectivity between northern populations and southern populations. The California current is southward driven during upwelling months (Checkley & Barth, 2009;Huyer, 1983;Largier et al., 1993) and water entrained in the Point Reyes eddy will eventually move offshore or south to Monterey (Rosenfeld et al., 1994;Steger et al., 2000). Oceanic current conditions along the coastline provide a potential for water transport north to south, which could connect northern and southern populations while reducing movement of water toward the San Francisco Bay gateway.

| Local adaptation of Clade Y to bay effluent conditions
Rather than divergence due to neutral processes, Leptasterias divergence could be the result of adaptive processes. Interestingly, Leptasterias patterns of clade distributions appear to coincide with regions of upwelling exposure in central California. Clade Y might be locally adapted to warm, low salinity conditions from San Francisco effluent affecting local coastal areas Smith & Cohen, 2013).
Intense upwelling zones span from Point Arena to Cape Mendocino (Bakun, 1990;Huyer & Kosro, 1987), and occur near Año Nuevo (Rosenfeld et al., 1994). Leptasterias aequalis K and L. aequalis B occur at upwelling exposed regions north of Point Reyes and south of Half-Moon Bay ( Figure 6). Bay-proximal populations are exposed to the warm, low salinity effluent from San Francisco Bay. Low haplotype diversities and negative neutrality statistics at the mitochondrial and nuclear loci used in this study could reflect selection upon other genes favoring Clade Y individuals at bay-proximal sites. Adaptive divergence is consistent with expectations of brooders that lack a highly dispersive life stage (reviewed by Sanford & Kelly, 2011;Sotka, 2012;Strathmann, 1986). While the genetic break of Leptasterias clades around the bay area appears to be upwelling associated, other factors associated with estuarine effluent may be causing further differentiation of Clade Y. Behavioral assays assessing the tolerance of clades to variable temperature and salinity conditions are an area of ongoing investigation (Contreras & Cohen, 2014;Shaw & Cohen, 2015, Braun et al., 2016Rupert, 2020), though are made difficult due to population declines from SSWD, discussed below.   Table 1) and colors represent clade delineation (see Figure 1) 1897-1974 1975-2000 2001-2014 (Jurgens et al., 2015).

| Conclusions
Beginning in 2013, SSWD, putatively attributed to a densovirus, caused massive population declines in Leptasterias and many other sea star genera (Harvell et al., 2019;Hewson et al., 2014, but see Hewson et al., 2018;Hewson et al., 2019;MARINe, 2015;Miner et al., 2018;Menge et al., 2016;Jaffe et al., 2019) and allele frequency shifts in Pisaster ochraceus (Schiebelhut et al., 2018). The frequency of these documented events over 10 years suggests local extinction events might have been common in the evolutionary history of Leptasterias as they have been for other echinoderms (Uthicke et al., 2009). The decline in abundance and range of Clades Y and Z illustrate brooding species' susceptibility to variation in population density. This study provides a record of population structure in Leptasterias sea stars over 117 years along the California coastline and can be used to understand changing population dynamics caused by large-scale mortality events.
Range shifts of native species have the potential to heavily impact community structure and function in regions experiencing expansions (Sorte et al., 2010) and poleward range shifts have been documented in many invertebrate species (Lonhart et al., 2019;Sanford et al., 2019;Sorte et al., 2010). We cannot rule out climate change as a potential factor in the range shifts found in this study. Leptasterias function as important predators in the intertidal environment by preying upon snails, limpets, barnacles, and juvenile mussels. Changes in population abundance of Leptasterias could impact standing algal stocks by indirectly affecting grazing by herbivores (Gravem & Morgan, 2019). Long-term changes in upwelling processes associated with climate change such as stronger upwelling-favorable winds, colder water, and a higher frequency of upwelling occurrences (García-Reyes & Largier, 2010) have the potential to impact selective and demographic forces that can lead to further shifts in population dynamics.
Historic genetic sampling has important implications in conservation management practices through monitoring population genetic diversity and interpreting environmental influences on diversity (Fenderson et al., 2020;Nielsen & Hansen, 2008). We found dramatic change in a genus of sea stars over a relatively short time span on the California coastline and we suggest several mechanisms for how the environmental landscape has shaped the recent evolutionary history of a low-dispersing sea star. We recommend further studies to understand the species delineation within this genus through morphological and physiological analysis. Additional monitoring of genetic diversity over time following sea star wasting disease, paired with this dataset, would be a valuable look at changing genetic diversity caused by mass mortality events.

ACK N OWLED G M ENTS
The authors dedicate this manuscript to Dr. John Pearse who provided inspiration and knowledge to so many. He encouraged us to tackle the Leptasterias species complex in a multifaceted way sharing his insights on morphology, behavior, and distributional variation on the shores of Monterey Bay. We would like to thank members of the Cohen Lab, Jason Helvey, and Shana Gallagher for help with sample collections and lab work. We would also like to thank John

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

DATA AVA I L A B I L I T Y S TAT E M E N T
Leptasterias spp. haplotype GenBank identification numbers can be found in Supplementary File 2. Leptasterias sample information, including IDs for museum samples, is found in Table 1.