Integrated Population Genomic Analysis and Numerical Simulation to Estimate Larval Dispersal of Acanthaster cf. solaris Between Ogasawara and Other Japanese Regions

The estimation of larval dispersal on an ecological timescale is significant for conservation of marine species. In 2018, a semi-population outbreak of crown-of-thorns sea star, Acanthaster cf. solaris, was observed on a relatively isolated oceanic island, Ogasawara. The aim of this study was to assess whether this population outbreak was caused by large-scale larval recruitment (termed secondary outbreak) from the Kuroshio region. We estimated larval dispersal of the coral predator A. cf. solaris between the Kuroshio and Ogasawara regions using both population genomic analysis and simulation of oceanographic dispersal. Population genomic analysis revealed overall genetically homogenized patterns among Ogasawara and other Japanese populations, suggesting that the origin of the populations in the two regions is the same. In contrast, a simulation of 26-year oceanographic dispersal indicated that larvae are mostly self-seeded in Ogasawara populations and have difficulty reaching Ogasawara from the Kuroshio region within one generation. However, a connectivity matrix produced by the larval dispersal simulation assuming a Markov chain indicated gradual larval dispersal migration from the Kuroshio region to Ogasawara in a stepping-stone manner over multiple years. These results suggest that the 2018 outbreak was likely the result of self-seeding, including possible inbreeding (as evidenced by clonemate analysis), as large-scale larval dispersal from the Kurishio population to the Ogasawara population within one generation is unlikely. Instead, the population in Ogasawara is basically sustained by self-seedings, and the outbreak in 2018 was also most likely caused by successful self-seedings including possible inbreeding, as evidenced by clonemate analysis. This study also highlighted the importance of using both genomic and oceanographic methods to estimate larval dispersal, which provides significant insight into larval dispersal that occurs on ecological and evolutionary timescales.

The estimation of larval dispersal on an ecological timescale is significant for conservation of marine species. In 2018, a semi-population outbreak of crown-of-thorns sea star, Acanthaster cf. solaris, was observed on a relatively isolated oceanic island, Ogasawara. The aim of this study was to assess whether this population outbreak was caused by large-scale larval recruitment (termed secondary outbreak) from the Kuroshio region. We estimated larval dispersal of the coral predator A. cf. solaris between the Kuroshio and Ogasawara regions using both population genomic analysis and simulation of oceanographic dispersal. Population genomic analysis revealed overall genetically homogenized patterns among Ogasawara and other Japanese populations, suggesting that the origin of the populations in the two regions is the same. In contrast, a simulation of 26-year oceanographic dispersal indicated that larvae are mostly selfseeded in Ogasawara populations and have difficulty reaching Ogasawara from the Kuroshio region within one generation. However, a connectivity matrix produced by the larval dispersal simulation assuming a Markov chain indicated gradual larval dispersal migration from the Kuroshio region to Ogasawara in a stepping-stone manner over multiple years. These results suggest that the 2018 outbreak was likely the result of self-seeding, including possible inbreeding (as evidenced by clonemate analysis), as large-scale larval dispersal from the Kurishio population to the Ogasawara population within one generation is unlikely. Instead, the population in Ogasawara is basically sustained by self-seedings, and the outbreak in 2018 was also most likely caused

INTRODUCTION
Many benthic marine invertebrates in coral reef ecosystems exhibit planktonic larval dispersal during their early life histories. Such larval dispersal connects different populations and forms a meta-population structure (Shanks, 2009). Thus, for effective conservation and management of benthic marine invertebrates, assessment of larval dispersal is important (Almany et al., 2009;Botsford et al., 2009). Because it is challenging to directly track the movement of numerous larvae that are smaller than 1 mm in the field, several different indirect methods, including population genetic analysis (e.g., Arndt and Smith, 1998;Yasuda et al., 2009), oceanographic simulation (e.g., Miyake et al., 2009;Storlazzi et al., 2017), towed nets (e.g., Uthicke et al., 2015;Yasuda et al., 2015a;Suzuki et al., 2016), and drifting buoys (e.g., Fukuda and Hanamura, 1996), have been used to estimate larval dispersal. Although each method has advantages and disadvantages, the first two methods can be applied to relatively large spatial scales and thus have been used more frequently than the other methods. On the one hand, the data obtained by population genetic analysis reflect the biological processes of larval dispersal and recruitment. However, these results not only include larval dispersal occurring on an ecological timescale but also reflect historical gene flow and genetic breaks caused by past climate change and plate tectonics (Benzie, 1999a;Crandall et al., 2014). Oceanographic simulation, on the other hand, can capture snapshots of larval dispersal, although even a sophisticated biophysical model cannot accurately simulate larval behavior, recruitment, and survival in the ocean (White et al., 2019). It is therefore important to integrate different methods to estimate larval dispersal (Marko and Hart, 2017), although few studies have attempted to do so (but see Alberto et al., 2011;Schunter et al., 2011;Nakabayashi et al., 2019;Taninaka et al., 2019).
In this study, we examined larval dispersal of a species of crown-of-thorns sea star, Acanthaster cf. solaris, which is a notorious predator of reef-building corals in the Pacific. Recent chronic outbreaks of crown-of-thorns sea stars have destroyed coral reef ecosystems in the Indian and Pacific oceans (Baird et al., 2013). For example, more than 90% of coral reefs over 38 km 2 in area were lost due to the feeding damage of crown-ofthorns sea star in Guam (Chesher, 1969). On Australia's Great Barrier Reef (GBR), more than half of the coral reefs have been reduced by 50% in size in the last 27 years and 42% of the decline of corals was due to population outbreaks of crown-of-thorns starfish (De'ath et al., 2012). The crown-of-thorns starfish is a species complex consisting of four genetically different lineages (Haszprunar et al., 2017); Pacific crown-of-thorns sea star (which is considered to be Acanthaster cf. solaris), and three other Indian ocean crown-of-thorns sea stars (Acanthaster planci in the north Indian Ocean and two more Acanthaster spp. in Red Sea and south Indian Ocean). Population outbreaks of crown-of-thorns sea star in the Indian and Pacific Oceans are largely categorized into two types. One is called "primary outbreak, " which is defined as an independent, sudden increase of crown-of-thorns sea star population at a location. The other is called "secondary outbreak, " which is defined as an outbreak caused by the spread of larvae released by another outbreak in neighboring regions. While the causes of primary outbreaks are still controversial, the mechanisms of secondary outbreaks are more obvious. With high fecundity and a larval stage that lasts a few to several weeks (Lucas, 1973;Yamaguchi, 1977), once a population outbreak takes place at a population, the large number of larvae produced and dispersed from the population can cause successive population outbreaks at neighboring sites, especially in the western Pacific regions that are impacted by strong warm currents (reviewed in Birkeland and Lucas, 1990;Pratchett et al., 2014;Yasuda, 2018). Given that early monitoring and control are key for management of A. cf. solaris (Birkeland and Lucas, 1990), assessment of the spatial range of secondary outbreaks of A. cf. solaris is necessary for predicting successive outbreaks and thus important for coral reef conservation.
Many population genetic and phylogeographic studies of A. cf. solaris have been conducted using various molecular markers (e.g., Benzie, 1992Benzie, , 1999aYasuda et al., 2009Yasuda et al., , 2015bTimmers et al., 2011;Vogler et al., 2013). These studies have indicated strong gene flow over a long distance along western boundary currents such as the Kuroshio Current region (genetic homogeneity from across the Ryukyu Islands to the temperate Pacific coast of mainland Japan) and the Eastern Australian Current region (the GBR, Nash et al., 1988;Benzie, 1999b;Vogler et al., 2013;Harrison et al., 2017), while relatively limited gene flow and genetic differentiation were observed among distant Pacific islands (Yasuda et al., 2009;Timmers et al., 2012;Vogler et al., 2013). Overall, the studies found strong gene flow and genetic homogeneity in regions with strong currents. Genetic data, however, strongly reflect historical gene flow that occurred during past climate change, and the spatial extent of secondary outbreaks could be overestimated (Benzie, 1999b;Yasuda et al., 2009, Yasuda et al., 2015b. Larval dispersal simulation studies in the GBR have indicated southward larval dispersal along the Eastern Australian Current, which is consistent with historically observed patterns of secondary outbreaks of A. cf. solaris (Dight et al., 1990a,b;James and Scandol, 1992;Scandol and James, 1992). Net surveys also have shown evidence of large-scale larval dispersal of A. cf. solaris along the GBR (Uthicke et al., 2015). While intensive studies of connectivity for A. cf. solaris have been conducted at the Indo-Pacific scale and in some localized regions where secondary outbreaks were suspected (e.g., the Kuroshio region and the GBR), connectivity between such regions and some relatively isolated oceanic island regions has yet to be examined.
In 2018, a population outbreak of A. cf. solaris on Ogasawara Island was reported (Biodiversity Center, Natural Environment Bureau, Ministry of the Environment, 2021), which was the second population outbreak on Ogasawara since 1979 (Kurata, 1984). Ogasawara Island is an oceanic island situated ca. 1,000 km south of the Tokyo urban center, and it is out of the mainstream of the Kuroshio Current. It is therefore still unknown whether this population outbreak is caused by mass larval recruitment produced by the outbreaks in the Kuroshio region.
The aim of this study was to estimate larval dispersal between the Kuroshio region and Ogasawara to determine whether the population outbreak in Ogasawara was successively caused by outbreaks in the Kuroshio region. To meet this aim, we used two different methods to estimate larval dispersal; population genomic analysis called multiplexed ISSR genotyping by sequencing (MIG-seq) and oceanographic larval dispersal simulation based on the Global HYbrid Coordinate Ocean Model (HYCOM).

Sampling and DNA Extraction
A total of 187 starfish samples were collected from six Kuroshio regions (Tatsukushi, Miyazaki, Sakurajima, Onna Village, Miyako, and Sekisei Lagoon) and Ogasawara via scuba diving ( Table 1 and Figure 1). The tube feet were preserved in 99.5% ethanol prior to DNA extraction. Genomic DNA extraction was carried out using a modified heat alkaline method (Meeker et al., 2007;Nakabayashi et al., 2019).

Whole-Genome Sequencing and Assembly of a Reference Genome
We obtained a new high-quality whole-genome sequence for mapping. Whole-genome sequencing of the sample collected from Miyazaki (BioSample ID, SAMD00056692) was performed using Sequel, which is a single-molecule real-time sequencer from Pacific Biosciences (PacBio; Menlo Park, CA, United States), according to the manufacturer's protocol. Previously published whole-genome short-read sequencing data (Wada et al., 2020) were downloaded from the BioProject PRJDB4009 in the ENA/Genbank/DDBJ database. Briefly, we obtained two pairedend (PE) libraries (insert sizes: 300 and 500 bp) and six mate-pair (MP) libraries (insert sizes: 3k, 5k, 8k, 10k, 12k, and 15k).
For the in silico procedures in this section, default parameters were used except where otherwise noted. Illumina PE and MP reads were filtered and trimmed using Platanus_trim v1.0.7. 1 We assembled trimmed PE and MP reads using Platanus v1.2.5 (Kajitani et al., 2014) with the commands "assemble, " "scaffold, " and "gap_close, " resulting in the scaffolds.
Gaps in the short read-based scaffolds (Platanus result) were filled with polished long read-based contigs using an in-house program. The procedure is as follows (Supplementary Figure 1): (1) Fixed-length regions flanking the gaps in scaffolds were extracted.
(2) The flanking sequences were aligned (queried) to the long read-based contigs using Minimap2 v2.17 (Li, 2018) with the options of "-c -k 19." (3) The alignment results were filtered out if sequence identity < 90%, query coverage < 25%, or multiple bestscoring hits were identified. (4) Each gap was filled with a region between the alignments in the contig if the following conditions were satisfied: (i) The flanking region pair was aligned to the same contig. (ii) The distance between alignments < 50 kbp.  (iii) The strands of the alignments were consistent.
(iv) There was no "N" in the corresponding regions in the contig.
Finally, to eliminate contamination, we performed homology searches of the assembly against the NCBI genome database (human, bacteria, and virus), the mitochondrial genome of A. planci (accession no. NC007788), and the chromosome sequence of the COTS symbiont (COTS27, accession no. AP019861) using BLASTN v2.7.1 (Altschul et al., 1990). The assembled sequences (scaffolds or contigs) were eliminated as contamination if sequence identity ≥ 90% and query coverage ≥ 50% (Supplementary Table 1

Multiplexed ISSR Genotyping by Sequencing Analysis
MIG-seq is a population genomic method that can easily detect a moderate number of neutral single nucleotide polymorphisms (SNPs) (Suyama and Matsuki, 2015). Previous studies using MIG-seq on marine species successfully discovered species boundaries and genetic structures that were undetectable using traditional genetic markers such as mitochondrial DNA and nuclear ribosomal internal transcribed spacer 2 (Nakabayashi et al., 2019;Takata et al., 2019).
We used eight pairs of multiplex ISSR primers (MIG-seq primer set 1; Suyama and Matsuki, 2015) with the Multiplex PCR Assay Kit Ver. 2 (TaKaRa Bio Inc., Shiga, Japan) for the first polymerase chain reaction (PCR). We slightly modified the protocol by Suyama and Matsuki (2015) and used a total reaction volume of 7 µL in a T100 TM Thermal Cycler (Bio-Rad, Hercules, CA, United States) under the following conditions: initial denaturation at 94 • C for 1 min; followed by 25 cycles of 94 • C for 30 s, 38 • C for 1 min, and 72 • C for 1 min; and final extension at 72 • C for 10 min. The samples were indexed in the second PCR, which was performed with PrimeSTAR GXL DNA polymerase (TaKaRa Bio Inc.) using a total reaction volume of 12 µL in a thermal cycler under the following conditions: 21 cycles of 98 • C for 10 s, 54 • C for 15 s, and 68 • C for 1 min. PCR products were run on an agarose gel, and 350-800 bp PCR products were extracted from the gel. The gel-extracted DNA was pooled and sequenced using MiSeq Control Software v2.0.12 (Illumina) with the MiSeq Reagent v3 150 cycle kit (Illumina). Image analysis and base calling were performed using Real-Time Analysis Software v1.17.21 (Illumina).
To eliminate low-quality reads and primer sequences from the raw data, we used the FASTX-Toolkit v 0.0.14 (Gordon and Hannon, 2012) 4 with a fastq-quality-filter setting of -q 30 -p 40. We removed adapter sequences for the Illumina MiSeq run from both the 5 end (GTCA GATCGGAAGAGCACACGTCTGAACTCCAGTCAC) and the 3 end (CAGAGATCGGAAGAGCGTCGTGTAGGGAAAGAC) using Cutadapt v1.13 (Martin, 2011), and then we excluded short reads less than 80 bp. We used Stacks v2.2 (Catchen et al., 2013;Rochette and Catchen, 2017) to stack the reads and extract the SNPs. We used the newly sequenced A. cf. solaris genome as a reference, as described below. We used the Gstacks option (-rm, -pcr, -duplicates) in Stacks v2.2 to remove PCR duplicates by randomly discarding all but one pair of each set of reads. We used population software in Stacks to prepare the dataset. We used the minimum percentage of individuals required to process a locus across all data (r) equal to 095. After excluding individuals missing more than 10% of SNP data, we used a single SNP option to avoid strong linkage between loci for Structure analysis. We used all SNPs for the remaining genetic analyses.

Population Genetics
BayeScan v2.1 (Narum and Hess, 2011)  GenoDive v. 3.0 (Meirmans, 2020) was used to estimate possible clones within each population. Because the maximum pairwise genetic distance across all the samples ignoring missing data was 62 and there was a clear disjunction between a few possible clonemates and normal distribution histogram of nonclonemates that start from pairwise genetic distance of 13 (Supplementary Figure 2), we heuristically set a limit of possible clonemate as genetic distance of six to identify clones. Probability of Identify (Table 1) using all loci for each population was also estimated using GenAlEx 6.5 (Peakall and Smouse, 2012).
To test whether the Ogasawara population is genetically different from other populations, a hierarchical analysis of molecular variance was conducted with prior grouping (Ogasawara vs. other populations) using Arlequin v3.5 (Excoffier and Lischer, 2010).
Bayesian clustering analysis was conducted using Structure software v2.3.4 (Pritchard et al., 2000) with 500,000 burnin followed by 500,000 steps of Markov chain Monte Carlo simulations and 10 iterations. We used admixture-and allele frequency-correlated models that were expected to best explain the data set. The LOCPRIOR model (Hubisz et al., 2009) was used with Ogasawara vs. other populations and without prior regional grouping, and the number of putative clusters (K) was set from 1 to 7. The online software Structure Harvester (Earl and von Holdt, 2012) was used to summarize the likelihood value. Structure plots were summarized using the online software CLUMPAK (Kopelman et al., 2015). To further visualize the genetic relationships among different populations, principal coordinate analysis was conducted using GenAlEx v6.502 based on pairwise genetic distances among populations. PGDspider v2.0.8.3 (Lischer and Excoffier, 2012) was used to convert the data formats.

Simulation of Larval Dispersal
We conducted a Lagrangian larval dispersal simulation to elucidate the oceanographic connectivity between the regions. For this simulation, we used Connectivity Modeling System v2.0 (Paris et al., 2013) with the ocean analysis/reanalysis products produced by the Global Ocean Forecasting System 3.0 of Global HYCOM (Chassignet et al., 2007), which are 1/12 • resolution (∼9 km) daily products archived from 1993 to 2018. 5 We extracted the surface water current and sea surface temperature (SST) datasets from the Global HYCOM products ranged from 115 • E-155 • E longitude and 15 • N-40 • N latitude. Thirty source and sink sites were chosen to create source sink polygons that fit purpose of this study (see Figure 1). The particles regarded as COTS larvae were numerically released from ocean areas inside the polygons and tracked. Once a particle passed more than 14 days after release entered a polygon, the particle was removed and judged as successful settlement, and the connection between the source polygon and sink polygon was counted. If a particle passed 35 days without entering a polygon, it was removed and scored as failure to recruit. The timing of particle release (start time of COTS mass spawning) and duration of spawning at each site were determined by the following rules using a history of SST at the site provided by Global HYCOM analysis/reanalysis products based on Yasuda et al. (2010): (1) If SST at the site rose above 28 • C, the particles were released until SST dropped below 26 • C or 30 days passed; (2) If SST at the site did not reach 28 • C, the particles were released while SST was higher than 26 • C. A total of 100 particles per day were released from each source polygon under the above conditions. Finally, we obtained a connectivity matrix that was created from larval dispersal simulations over 26 years .
Because the connectivity based on simulated larval dispersal is a one-generation direct linkage, indicating that stepping-stone migration had not been considered, we conducted back-tracing to the origin through a discrete-time Markov chain using the connectivity matrix. To this end, we first set the following form of the connectivity matrix N(m × m): The matrix elements were normalized by the total number of sink particles as follows: Where The normalized connectivity matrix can be regarded as the transition matrix of the Markov chain. The n-th power of the transition matrix P n represents the n time-step back-trace to the origin source.

Sequencing and Single Nucleotide Polymorphism Metrics
In total, 11,461,866 raw reads ranging from 9,190 to 281,474 reads per sample were obtained for 187 A. cf. solaris samples. After filtering out the low-quality reads, 6,125,996 reads ranging from 5,996 to 107,468 reads per sample remained. On average, 189,536 reads per individual remained after two-step filtering. We excluded seven samples missing more than 10% data and used a total of 182 samples for analysis. Using all SNP options, we obtained a total of 464 SNPs, which were used for all subsequent analyses except for Structure analysis. For Structure analysis, 115 SNPs using a single SNP option were used. BayeScan indicated that all loci were neutral (q-values > 0.05).

Genetic Diversity and Genetic Differentiation
Genetic diversity (H O and H E ) of the studied populations, including that from Ogasawara, were comparable and ranged from 0.039 to 0.044 and 0.039-0.045 for H O and H E , respectively ( Table 1). In total five pairs of possible clonemates were found in Tatsukushi (TTK, one pair out of 12), Onna Village (ONN; one pair out of 41), and Ogasawara (OGS; three pairs out of 30) populations (Table 1). Global F ST values were significant (P < 0.001) but small both for groups including (F ST = 0.005) and excluding the Ogasawara population (F ST = 0.006). Likewise, analysis of molecular variance with prior grouping (Ogasawara vs. other populations) indicated that F CT was low and not significant (0.00167, P = 0.219). Pairwise F ST values across the studied populations were also small, and only two pairs [Sekisei Lagoon (SKS) and Tatsukushi (TTK), Miyazaki (MYZ), and Ogasawara (OGS)] became significant after sequential Bonferroni correction (α < 0.05) ( Table 2).
Mean log likelihood values [mean lnP(D)] were highest at K = 1 without prior grouping and comparable at K from 1 to 7 with prior grouping, with the lowest standard deviation value of mean lnP(D) found at K = 1 (Supplementary Table 2). Structure bar plots with and without prior geographic grouping suggested genetic homogeneity across the populations (Figure 2A and Supplementary Figure 4). Principal coordinate analysis indicated that, along the x-axis, SKS, MYK, and Sakurajima (SKR) were divided from the others, although it was not related to geographic location or the year of sampling ( Figure 2B).

The Larval Dispersal Simulation
The results of the larval dispersal simulations are shown as a connectivity matrix (Supplementary Figure 3). The diagonal elements of the matrix indicate self-seeding. Basically, site numbers from 1 to 20 were assigned from south to north along the Kuroshio Current; site numbers from 21 to 25 were assigned from Nii-Jima and Shikine islands, which are closer to Japan's main island, to the Ogasawara and Iwo-Jima Islands, which are further; and site numbers from 26 to 30 were assigned along the Tsushima Warm Current. The linkage of the site numbers for numerical simulation to site names for genomic analysis was as follows: 5 = Sekisei Lagoon (SKS), 6 = Miyako (MYK), 7 = Onna Village (ONN), 11 = Miyazaki (MYZ), 13 = Tatsukushi (TTK), 25 = Ogasawara (OGS), and 27 = Sakurajima (SKR). The results of this analysis confirmed that the larvae were transported from south to north by following the Kuroshio Current. However, they  also confirmed that self-seeding dominated in the Ogasawara area, and it is very rare for larvae from the other sites to reach the Ogasawara area directly or for larvae from the Ogasawara area to reach the other sites directly. Figure 3 shows the back-tracing results obtained using Markov chains. The bar of each site indicates the ratio of the number of particles released from the source site (laws of transition matrix P n ). The results indicate that the origin sites gradually converged with some of the upstream sites of the Kuroshio Current. The population in the Ogasawara area is not well mixed with that of the Kuroshio region, but the 200-500 time-steps, which may be regarded as representing over a few hundred generations, confirmed that the populations in the Kuroshio region and the Ogasawara region are almost mixed.

DISCUSSION
This study examined larval dispersal of the coral predator A. cf. solaris between the Kuroshio region and Ogasawara using population genomic analysis and oceanographic simulation to test the secondary outbreak hypothesis. Contrasting results were obtained using different methods. Population genomic analysis indicated genetic homogeneity between the populations of the two regions, and oceanographic numerical simulation indicated that multiple-generation stepping-stone migration was required to disperse larvae from the Kurshio region to Ogasawara. These results indicate that the Ogasawara population of A. cf. solaris has the same genetic origin as that of the Kuroshio region, but a large amount of one-generation larval dispersal, such as a secondary outbreak from the Kuroshio region to Ogasawara, is unlikely. Rather, the data suggest that the population outbreak in Ogasawara might have been caused by successful selfseeding. This study highlights the importance of an integrated approach for estimating larval dispersal, including its timescale, for conservation purposes.

Potential of Secondary Outbreak Between Ogasawara and Kuroshio Region
In this study, all population genomic results indicated genetic similarity and no genetic structure between A. cf. solaris populations in the Ogasawara and Kuroshio regions. The analysis of molecular variance and Structure results indicate that the origins of A. cf. solaris populations in the Ogasawara and Kuroshio regions are genetically similar. Low F ST /F CT values between the Ogasawara and Kuroshio regions indicate (1) a contemporary, large amount of larval dispersal, (2) accumulation of stepping-stone migrants over multiple years, or (3) recent separation of a population with no or limited ongoing gene flow. Based on the oceanographic numerical simulation analysis results and historical records of A. cf. solaris in Ogasawara, the latter two are more likely than contemporary, large amounts of larval dispersal.
The result of the 26-year oceanographic simulation suggests that the chance of direct one-generation larval dispersal between the Kuroshio region and Ogasawara is very low, although stepping-stone migration (e.g., via Hachijyo-jima in Izu islands) between the two regions over multiple generations is physically possible. The Markov chain simulation, which simulated stepping-stone migration over multiple generations, indicated that the Ogasawara population is dominated by migrants from the Kuroshio region after 200-500 generations. This implies that the accumulation of larvae from low levels of migration between the Kuroshio and Ogasawara regions over multiple years would result in genetic homogeneity. Given that there are no observed reports of population outbreaks of A. cf. solaris in stepping-stone regions (e.g., Izu islands) in the last 40 years (Yasuda, 2018), genetic homogeneity between the Ogasawara and Kuroshio regions observed by F-statistics and Structure analysis could be attributed to stepping-stone larval migration over multiple generations. Crandall et al. (2018) demonstrated observed results of the marine species along Hawaii archipelagos where F ST is not significantly different from zero were not actually panmictic but having weak hierarchical structure of isolation-by-distance caused by stepping-stone migration. Coalescent simulation including A. cf. solaris samples from Izu islands would confirm this hypothesis.
Alternatively, it is possible that the observed strong gene flow resulted from recent colonization. According to a literature survey, there was no record of A. cf. solaris in the Ogasawara region before 1945, and the first record of this species in Ogasawara was in 1968 (Kurata, 1984;Yasuda, 2018). Since then, the density of A. cf. solaris in Ogasawara has been low (Yasuda, 2018). However, a lack of a record of A. cf. solaris in Ogasawara before 1945 does not necessarily mean a complete absence of this species in the region, given that almost no studies or surveys were conducted before 1945. In addition, crown-of-thorns sea stars are cryptic, and it would have been difficult to detect a few individuals that have been colonized for the first time. Recent migration with a small population size would predict low genetic diversity due to the founder effect, although we observed almost the same genetic diversity in the Ogasawara and Kuroshio regions. Therefore, high genetic diversity in Ogasawara and strong gene flow between the Ogasawara and Kuroshio regions is more likely caused by stepping-stone larval migration over multiple generations, as suggested by the oceanographic simulation.
Whichever may be the case, the oceanographic simulation suggests that contemporary, one-generation larval dispersal between the Ogasawara and Kuroshio regions is limited, and thus, direct secondary population outbreaks that require a large amount of larval dispersal within a generation from the Kuroshio region to Ogasawara are unlikely. Kurata (1984) speculated that A. cf. solaris in Ogasawara might have colonized either from the Kuroshio region or from Mariana regions such as Guam. According to a population genetics analysis using microsatellite markers, genetic differentiation between the Kuroshio and Guam regions has been identified (Tusso et al., 2016), implying that colonization and/or secondary population outbreaks from Guam to the Ogasawara region are less likely than those from the Kuroshio region.
Only one previous study examined the genetic connectivity of coral reef organisms between the Kuroshio and Ogasawara regions. A reef-building coral species, Acropora sp., which has a shorter (average 3-4 days) larval stage than that of A. cf. solaris, showed genetic differentiation between the Ogasawara and Kuroshio regions (Nakajima et al., 2012). This implies that larval dispersal from the Kuroshio region to the Ogasawara region of coral reef organisms with short larval stages is quite limited.

Possible Causes of Population Outbreak in Ogasawara
Oceanographic simulation suggests that many larvae are selfrecruited in Ogasawara compared with other A. cf. solaris populations in the Kuroshio region (Supplementary Table 3). This fact indicates that the observed population outbreak in Ogasawara in 2018 was caused by a sudden increase in regional self-seeding (namely "primary outbreak" as opposed to "secondary outbreak"). While several hypotheses have been proposed to explain the primary outbreak of A. cf solaris, the causes of primary outbreaks are still unclear (Pratchett et al., 2014). The causes of outbreak would also depend on several regional environmental conditions being met simultaneously Okaji et al., 2019): elevated levels of larval nutrition and subsequent higher survival rates supported by favorable physical, chemical, and biological conditions such as lower predation on larvae and juveniles; and physical processes acting on larvae that increase their survival (e.g., weather conditions and current patterns) (Keesing and Halford, 1992;Okaji et al., 2019).
The population outbreak in Ogasawara in 2018 could be attributed to successful local larval recruitment in 2016, given that it takes 2 years for A. cf. solaris larvae to become adults (Birkeland and Lucas, 1990). Oceanographic simulation also suggests that the number of self-recruited larvae in Ogasawara was larger (397) in 2016, which exceeded the average (252) over 26 years (Supplementary Table 4). The large number of selfrecruited larvae might be partly due to the smaller number of typhoons during the larval period (July to October, estimated by the spawning period associated with water temperature; Yasuda et al., 2010); in 2016, three typhoons passed through the Ogasawara region, which is below the 70-year average (4.47) (data from Japan Meteorological Agency, Supplementary Table 5). A previous field survey discovered dense clouds of A. cf. planci larvae before one typhoon, but they were scattered and disappeared from the coral reef area after the typhoon (Suzuki et al., 2016), implying that typhoons prevent mass settlement of larvae. However, in some years with a higher number of self-recruits and fewer typhoons, no population outbreak was observed in Ogasawara.
Recently, self-fertilization and/or clonal propagation hypothesis has been proposed by Guerra et al. (2020) to explain the cause of primary outbreak as evidenced by the existence of hermaphrodites and the finding of possible larval cloning of A. cf. solaris (Allen et al., 2019). In the present study, we found three pairs of individuals in Ogasawara population were possible clonal or genetically very close individual (Table 1 and Supplementary  Figure 2). It is also notable that we also found one pair of possible clonemates in Tatsukushi population where previous study found hermaphrodites (Guerra et al., 2020). Although we need more precise analysis to identify true clonemates, it is possible that inbreeding including self-fertilization or larval cloning also partly facilitated the initiation of population outbreak in the oceanographically isolated Ogasawara populations.
Other factors such as increased success of cross-fertilization (Rogers et al., 2017), lower salinity, and higher phytoplankton biomass during the larval stage (Lucas, 1973), and reduced predation during the post-settlement stage (McCallum, 1987;Keesing and Halford, 1992) also may be associated with the demographic fluctuation of A. cf. solaris in the Ogasawara area. An Allee threshold of three starfish ha −1 (below which reproductive capacity is greatly reduced regardless of the aggregation level) is proposed and aggregation during spawning period is key in causing population outbreaks (Rogers et al., 2017). Laboratory experiments indicated that lower salinity and increased phytoplankton biomass that are often caused by territorial run-off enhances the survival rate of A. cf. solaris (Lucas, 1973;Birkeland, 1982;Fabricius et al., 2010). Keesing and Halford (1992) also suggested that juveniles of A. cf. solaris would be strongly influenced by predation because of their slow growth and movement, underlining the importance of posttreatment processes for population outbreak. All these factors may influence the population dynamics of A. cf. solaris in Ogasawara.

Multiple Approaches to Estimate Larval Dispersal in Marine Organisms
This study estimated larval dispersal of A. cf. solaris between Kuroshio and Ogasawara using two different approaches: population genomic analysis and larval dispersal simulation. Gene flow estimation by the population genomic approach often reflects both contemporary and evolutionary processes, rendering the interpretation of the amount of contemporary dispersal difficult because of the long evolutionary timescale. Genetic homogeneity can be observed due to a large amount of contemporary larval dispersal (secondary population outbreak) or stepping-stone migrants over multiple generations or recent divergence with no contemporary dispersal. Our study is unique in that it included an oceanographic model with Markov chain simulation. Such a simulation was useful for testing these alternative hypotheses. A simple oceanographic simulation averaged over 26 years first rejected the secondary population outbreak hypothesis by indicating limited physical larval transport from the Kuroshio to the Ogasawara region within a generation. This result supports the idea that steppingstone migrants over multiple generations is responsible for the genetic homogeneity observed in the F ST and Structure analysis between the Kuroshio and Ogasawara populations.
This may give the impression that oceanographic simulation alone is enough for testing the secondary population hypothesis; however, it is generally imperfect for estimating larval dispersal. Our model regarded larvae as neutral passive particles without any of the biological characteristics of larvae that may influence actual larval dispersal, such as vertical movement, natural selection, food availability, and mortality. In addition, the model only simulates the dispersal process, which does not include post-settlement survival. In some marine species, genetic connectivity is quite limited despite having long pelagic larval duration possibly caused by several biological features of larvae and juveniles (e.g., Barber et al., 2000Barber et al., , 2002. Thus, a comparison with genetic connectivity that reflects all the biological characteristics of the species is also important. While several previous studies have reported good agreement between oceanographic models and genetic analysis and have provided robust and straightforward interpretations of larval dispersal occurring on an ecological timescale (e.g., Galindo et al., 2006;White et al., 2010;Sunday et al., 2014;Taninaka et al., 2019), only a few have discussed the usefulness of oceanographic modeling for identifying alternative hypotheses or provided additional insights into the causes of observed population genetic structure (Galindo et al., 2010;Crandall et al., 2014). By applying Markov chain simulation in the numerical simulation, we could also include the effect of migration over multiple generations, which is easier to compare with genetic results than a single generation simulation result in terms of timescale. This method underlines the usefulness and importance of using both genetic and oceanographic methods to estimate larval dispersal and provides significant insight into larval dispersal that occurs on ecological and evolutionary timescales.

Conservation Implications
This study indicated that population outbreaks of A. cf. solaris in the Ogasawara Islands are independent of successive outbreaks in the Kuroshio region and that outbreaks are more likely caused by local self-seedings. Considering it takes at least 2 years for A. cf. solaris to reach adulthood (Birkeland and Lucas, 1990), caution and early monitoring should be taken 2 years following the observation of local ocean current system retention and other preferable conditions for the survival during early life stage of A. cf. solaris. For example, a large amount of terrestrial run-off will increase the survival rate of A. cf. solaris by increasing food available for larvae. In addition, high levels of larvae aggregation were observed only before a typhoon (Suzuki et al., 2016), and thus, a fewer number of typhoons during spawning periods may increase the chance of high levels of larval aggregation and settlement. Local monitoring of juvenile starfish 2 years after such preferable conditions would be useful for early detection and predicting population outbreaks and to avoid being too late for coral reef management (Okaji et al., 2019).

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: DDBJ (accession: DRA011826).

AUTHOR CONTRIBUTIONS
NY conceived the study. MH, NY, RK, HY, and TY drafted the manuscript. RK, HY, TI, and AT sequenced, assembled, and reconstructed the reference genome. NY, YA, and TS collected the Acanthaster samples. MH, HT, HY, TK, and NY conducted the MIG-seq analysis. TN and TY conducted oceanographic simulations. All authors edited the draft and approved it for publication.

ACKNOWLEDGMENTS
We express our gratitude to Eriko Nagahiro for her help in the molecular experiment.