Higher Genetic Diversity of the Common Sea Cucumber Holothuria (Halodeima) atra in Marine Protected Areas of the Central and Southern Ryukyu Islands

Recently, sea cucumbers (Echinodermata: Holothuroidea) have been over-exploited in many areas of the world, including in the Ryukyu Islands, southern Japan, due to increases in their economic importance. Nevertheless, management and protection of sea cucumbers are insufficient worldwide. The black sea cucumber Holothuria (Halodeima) atra Jaeger, 1833, inhabits a large range across the Indo-West Pacific Ocean and is a widely harvested species. Here we conducted population genetic analyses on H. atra using partial mitochondrial DNA sequences of cytochrome c oxidase subunit I (COI) and 16S ribosomal RNA (16S) to examine 11 different populations around three island groups in the middle Ryukyus; Okinawajima Island, the Kerama Islands, and the Sakishima Islands, all within Okinawa Prefecture. We found 27 haplotypes for COI and 16 haplotypes for 16S. Locations within national and quasi-national parks (Zamami Island, Keramas, and Manza, Okinawajima; managed by the national Ministry of Environment and Okinawa Prefecture, respectively) had the highest number of haplotypes, whereas locations with less management and more anthropogenic pressure had lower numbers The mean of all samples' genetic diversity indices was moderate with regards to both haplotype and nucleotide diversity. According to our results, Zamami Ama was the most genetically diverse location based on both markers used, likely because it is located within Kerama-Shoto National Park with comparatively stricter regulations than most other locations. Based on our COI sequences, three-quarters of the locations with the highest haplotype diversity were found to be distant from Okinawajima Island, indicating that the genetic diversity of H. atra was reduced around Okinawajima Island. Our results possibly reflect negative impacts from anthropogenic pressures such as over-harvesting and coastal development, although future comprehensive research including sequences of nuclear loci is needed to confirm this hypothesis.


INTRODUCTION
Sea cucumbers (Echinodermata, Holothuroidea) are distributed across the world's oceans, from shallow coastal areas to the deep sea (Tyler et al., 1992;Purcell et al., 2016). Roughly 20 species of sea cucumbers are known to perform transverse fission (Dolmatov, 2014), and in some species, a large fraction of individuals produce asexual propagules (e.g., up to 76% of Holothuria atra experience fission each year at Fantome Island on the Great Barrier Reef; Uthicke, 1997). Among sea cucumbers, the order Holothuriida is widely distributed in tropical and subtropical waters (Purcell et al., 2012;Miller et al., 2017). Holothuriida species generally ingest small-sized organic matter and microalgae in sediment using their tentacles (Yingst, 1982;Mercier et al., 1999;Mfilinge and Tsuchiya, 2016;Purcell et al., 2016). It is well known that Holothuriida can enhance seagrass meadow productivity (Wolkenhauer et al., 2010), microphytobenthos productivity (Uthicke and Klumpp, 1998), calcium carbonate dissolution (Schneider et al., 2011(Schneider et al., , 2013, and can also buffer ocean acidification by discharging ammonia (Uthicke, 2001;Purcell et al., 2016). For all of these reasons Holothuriida species are considered ecosystem engineers that play crucial roles in diverse marine ecosystems, including coral reefs (Purcell et al., 2012).
Some Holothuriida species are prized as delicacies in Asia (Conand, 2004) and are therefore important commercially harvested species. However, the protection and management of these organisms have generally been insufficient. In recent years, overfishing of sea cucumbers has become a serious problem worldwide due to rising market demand (Uthicke and Conand, 2005;Toral-Granda et al., 2008;Purcell et al., 2013). Sea cucumbers in Japanese coastal waters have not been exempt from such problems, and some species have been over-harvested along the Japanese coastline for consumption, both domestically and for export to other Asian countries (Akamine, 2005).
Holothuria (Halodeima) atra (Jaeger, 1833), is one of the most common sea cucumber species along Indo-West Pacific coastlines and inhabits coral rubble, sandy tidal flats, and reef slopes (Purcell et al., 2012). The pelagic larval duration (PLD) of H. atra is estimated to be around 20 days (Ramofafia et al., 1995), and this time span is thought to be long enough for larvae to cross geographical regions. This species is known to have small and large morphotypes, but it is not yet clear whether they correspond with genetic variation . Although this species is not as commercially important as some other species such as Holothuria (Metriatyla) scabra Jaeger, 1833 or Thelenota ananas Jaeger, 1833, demand for H. atra has increased due to the depletion of other, more highly-valued species in the western central Pacific region (Kinch et al., 2008). Currently, H. atra is harvested in Micronesia, Polynesia, Melanesia, Australia, and New Zealand (Kinch et al., 2008), as well as in many Asian countries (Choo, 2008). Commercial harvesting of this species has also been reported in Okinawa Prefecture, southern Japan, although information on the local stock conditions is sparse (Purcell et al., 2014). Since H. atra is common in Okinawa and the Ryukyu Islands, investigating this species' population health is crucial for sustainable fishery management as well as for better preservation of local ecosystems.
In order to assess whether changes occur in the population's genetic composition of ecologically important organisms, molecular methodologies have been applied to marine organisms. Such methods have also been applied to holothurian species, and some studies, in various regions, have been performed examining their genetic diversity and population structure, mostly based on COI sequences. In Fiji, despite the number of specimens examined being low at some locations (three of the four populations considered had no more than 10 specimens each), only eight haplotypes were observed among 40 specimens. Such low genetic diversity has been proposed to be due to asexual recruitment processes such as fission (Eastwood et al., 2016). In contrast, Skillings et al. (2014) found that the haplotype diversity of 252 individuals of H. atra in Hawai'i was relatively high compared to marine population genetic surveys of other species in the same area.
The Kuroshio Current, which is a warm tropical current that runs north along the Ryukyu Islands in southern Japan, is known as one of the most important oceanographic features across this region (Barkley, 1970), and is thought to be important in structuring the population of coral reef species in this region (Yasuda et al., 2014: Zayasu et al., 2016. In the central Ryukyus around Okinawajima Island, population genetic research has been conducted on local scale on the sea cucumber species Holothuria (Halodeima) edulis Lesson, 1830(Soliman et al., 2016a and Stichopus chloronotus Brandt, 1835(Soliman et al., 2016b. These studies have indicated that even within small geographic areas around Okinawajima Island, genetically different populations can be distinguished in these species. Regarding H. atra, a previous study explored the genetic diversity of 44 specimens from two locations in southern Japan. However, only one population has been investigated from Okinawa (Skillings et al., 2011), and the local population connectivity of H. atra in Okinawa has still not been well studied. As well, research on sea cucumber population genetics on a larger scale to examine the effects of the Kuroshio Current have not been conducted. The knowledge of local genetic population structures is critical to ensure effective conservation measures (e.g., Fiji, López et al., 2017), such as in establishing conservation areas, such as Management Units (MUs). Since the global demand of holothurians has increased rapidly, deposition of such information is becoming increasingly important for better management. Thus, using partial sequences of cytochrome c oxidase subunit I (COI) and 16S ribosomal DNA (16S), we investigated the population genetic diversity and structure of H. atra around Okinawajima Island and nearby outlying islands in the central and southern Ryukyus. Finally, by investigating the genetic population diversity among locations and comparing population haplotype diversity of each island and site in this research, we discuss MUs. MUs are defined as "populations with significant divergence of allele frequencies at nuclear or mitochondrial loci, regardless of the phylogenetic distinctiveness of the alleles" (Moritz, 1994). Such populations are connected with weak gene flow, but may be possibly independent of each other, and thus the consideration of such units is quite important in the maintenance of potential Evolutionally Significant Units (ESU) (Moritz, 1994). Here, we define potential MUs for H. atra inhabiting Okinawa.

Sampling
Holothuria atra specimens were collected by reef walking and/or snorkeling between April to September 2019 at 11 locations across the central and southern Ryukyu Islands, from three island groups, Okinawajima Island (Chatan, Kayou, Kin, Manza, Odo, Sesoko and Uruma), the Kerama Islands (Zamami Ama and zamami Port), and the Sakishima Islands (Ishigaki St. 20 and Ishigaki St. 27) (Figures 1-3; Table 1). Zamami Island and Ishigaki Island are included in Kerama-Shoto and Iriomote-Ishigaki National Parks, respectively, whereas Manza is in Okinawa-Kaigan Quasi-National Park (Figure 1). Individual H. atra specimens were collected by hand, and a small amount of tissue was cut from around the mouth of each individual. All collections were non-lethal, and the individuals were returned to the original location from where they were taken. This nonlethal sampling method was confirmed to not require permits with the relevant prefectural administrator. In between sampling individuals, scissors were cleaned with 99% EtOH in order to avoid contamination. Tissues obtained were preserved in 99% EtOH until further investigation.
PCR was conducted in a 20 µl reaction volume, with 7 µl of purified water, 10 µl of HotStarTaq R Plus Master Mix (Qiagen, Tokyo), 1 µl each of forward and reverse primers, and 1 µl of DNA. Amplification was performed under the following conditions: 95 • C for 15 min, followed by 35 cycles of 30 s at 94 • C, 45 s at 50 • C and 1 min at 72 • C, with a final extension of 10 min at 72 • C. PCR clean-up was done with the Shrimp Alkaloid Phosphate (SAP) method following the manufacturer's protocol. For SAP, 0.3 µl of SAP solution, 0.15 µl of Exonuclease I solution and 2.55 µl of purified water were added, and then samples were incubated for 20 min at 37 • C followed by 30 min at 83 • C. Sequencing reactions were performed using a BigDye TM Terminator v3.1 Cycle Sequencing Kit (ThermoFisher). The reaction mixtures contained a final volume of 10 µl and included 2.5 µl of BigDye 3.1 Master Mix and buffer, 0.5 µl of 10 µM primer, 1.5 µl of PCR product, and 5.5 µl of water. The following reaction conditions were applied: initial denaturation at 96 • C for 1 min followed by 25 cycles of 96 • C for 10 s, annealing at 50 • C for 5 s and extension at 60 • C for 4 min. Sequencing products were purified using Mag-Bind SeqDTR magnetic beads (Omega Bio-tek) following the manufacturer's instructions. The resulting purified products were analyzed on an ABI3130 Genetic Analyzer (Applied Biosystems, ThermoFisher) at the Genomic Unit, Scientific and Technological Support Center for Research (CACTI), University of Vigo, Spain. The obtained sequences were visually checked, assembled and aligned with Geneious v.8.1.9 (Kearse et al., 2012). Subsequently, COI and 16S sequences were compared with previously reported H. atra sequences using the Basic Local Alignment Search Tool (BLAST) on GenBank (www.ncbi.nlm.nih.gov/GenBank). The obtained sequences were then analyzed with those available from GenBank (clade 1-EU848217, EU848265, LC217308; clade 2-EU848222, EU848283, LC217309; clade 3-EU848266; also see Supplementary Table 2). After including these reference sequences, a COI alignment consisting of 260 sequences and 559 nucleotide positions was constructed. The partial 16S sequence of H. edulis (MZ081847) was added to the 16S alignment. Thus, 16S alignment consisted of 235 sequences and 349 nucleotide positions. For phylogenetic analyses, the best-fit substitution models (K2P+G for COI, T92+G for 16S) were selected with MEGA v.7 (Kumar et al., 2016). The Maximum Likelihood (ML) phylogenetic trees were built for both datasets using MEGA. Support for branches was assessed using a rapid bootstrap analysis with 1,000 pseudo-replicates. The Bayesian tree was also inferred with MrBayes v.3.2.7 (Ronquist et al., 2012), under K2P+G model. The MCMC ran for 1,000,000 generations, sampling every 1,000 steps, with 25% of the generated trees discarded as burn-in. Bootstrap values >50 and posterior probabilities values >0.95 are shown in the ML tree (Figure 4).

Genetic Diversity Indices and Neutrality Tests
Genetic diversity indices, including haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), and polymorphic sites (S), were calculated with DNAsp v.5.10.01 (Librado and Rozas, 2009) under default options for both datasets. Based on the results of the haplotype estimation, a median-joining network was generated for visualization of connectivity with PopART v.1.7 (Leigh and Bryant, 2015). Pairwise distance analysis was

Local Fishing Pressure Analyses
We calculated the number of fishers per total coast length at our locations by acquiring data on fisher populations, and numbers of fishing boats in publicly available data in the "Basic plan for coastal conservation of Ryukyu Islands" (Okinawa Prefecture Office, 2008), "Document of the Remote Islands No. 2, Industry" (Okinawa Prefecture Office, 2019) and the "47th Statistical Annual Report of Agriculture, Forestry and Fisheries of Okinawa Prefecture" (Okinawa General Bureau, 2019). The total population of fishers for each area was divided by the coastline length of each area (Zamami Village, Ishigaki City, Okinawajima Island). We also calculated the amount of the catch of "others" per coastline length. The category "others" includes all collected marine animals aside from fishes, mammals, and mollusks for Zamami Village and Ishigaki City, and fishes, mammals, mollusks, cephalopods, crustaceans, and echinoids for Okinawajima Island (Okinawa General Bureau, 2019). The total "others" catch (Zamami Village and Ishigaki City−2016 report; Okinawajima−2017 report) of each area was divided by the coastline length of each area (Zamami Village, Ishigaki City and Okinawajima Island).

Management Units Investigation
To consider potential Management Units (MUs) of H. atra at the locations investigated, we explored the genetic diversity based on COI for each island, also considering the locations excluded from all other population genetic analyses due to small sample sizes (e.g., Miyako Island <10 specimens for each population; Supplementary Table 1). For this step, the same as in the separated locations' analyses, DNAsp v.5.10.01 (Librado and Rozas, 2009) under default options was employed to calculate genetic diversity indices including haplotypes (H), haplotype diversity (Hd), nucleotide diversity (π), and polymorphic sites (S). Based on the results obtained for haplotype estimation, a median-joining network was generated for visualization of connectivity with PopART v.1.7 (Leigh and Bryant, 2015).

Sampling Locations and Numbers of Specimens
From a total of 295 specimens, 253 COI and 234 16S sequences were successfully sequenced and used in further analyses ( Table 1, Supplementary Tables 1, 2). Although sampling was performed in other areas, including Oura Bay and several locations at Miyako Island, specimens from these locations were excluded from most analyses because of small sample sizes (<10).

Genetic Distances and Phylogenetic Tree Reconstruction
ML trees reconstructed based on both COI and 16S sequences did not show apparent correspondence with either location of specimens or body size of each specimen. The Bayesian tree reconstructed based on COI sequences showed a different topological structure compared to the ML tree, but with some nodes showing high posterior probabilities (Figure 4). Although the phylogenetic tree based on 16S showed no apparent clades Bayesian support values are shown on the right. Each color indicates clades as reported in a previously published phylogeny . Note several haplotypes includes multiple samples. Please see also Supplementary Table 2 to match samples to each haplotype. Figures 2, 3), the phylogenetic tree based on COI sequences was composed of two major groups, one consisting of nine haplotypes (clade 2, bootstrap = 53%) and another including 17 haplotypes (clades 1, bootstrap = 62%) (Figure 4, also see Supplementary Figure 1). Haplotype 20 was sister to EU848266 into clade 3 (bootstrap = 99%). The greatest genetic distances for COI sequences were between haplotype 20 and haplotypes 6, 11, and 27 (4.1%; 23/559 bp), whereas for 16S, the greatest distances were between haplotype 8 and haplotypes 10 and 16 (2.3%; 8/349 bp). Haplotype 20 was found only at Zamami Ama, whereas the haplotypes in clade 2 appeared at all locations, and those within clade 1 at all locations except Ishigaki St. 27 (Figure 5).

Population Genetic Structure
Pairwise Φ st estimates based on COI sequences showed Chatan and both Ishigaki H. atra populations to be distinct from the other locations examined ( Table 2). Zamami Ama was also significantly different from all other locations, including from Zamami Port despite being only ∼1 km apart. Uruma showed significant differences for most comparisons except with Kayou and Odo. All Okinawajima Island populations, excluding Chatan and Uruma, did not show significant differences in any comparison.
When considering 16S sequences, Zamami Ama was significantly different from all other populations, including from the other population of Zamami Island (Zamami Port). Among Okinawajima populations, Chatan showed significant differences when compared to all other locations except with Uruma, and also the Uruma population had significant differences from Kin, Manza, and Sesoko.

Population Dynamics
Neutrality tests based on COI sequences indicated that only Chatan showed a significant negative p-value in its Tajima's D test (p < 0.001), while Zamami Ama had a marginally not significant value (p = 0.068) ( Table 4). All Fu's FS values were positive except for that of the Zamami Ama population, yet none were statistically significant ( Table 4).
The results based on 16S sequences were not fully consistent with the results obtained with COI sequences. Tajima's D of Sesoko and Chatan populations showed negative values (p = 0.044 and 0.032 respectively), whereas only Zamami Ama showed negative Fu's FS value (p = 0.001) ( Table 4).

Potential MUs Establishment
Among the islands investigated, Miyako Island showed the highest haplotype diversity (Hd = 0.98718) followed by Zamami Island (Hd = 0.94168), Ishigaki Island (Hd = 0.69737) and Okinawajima Island (Hd = 0.63980). Surprisingly, almost all the specimens from Miyako Island had different haplotypes (12 haplotypes within 13 individuals), and thus potential haplotype diversity at the area was assumed to be very high (Supplementary Table 3). Accordingly, each location possessed many unique haplotypes (Zamami Island: nine unique haplotypes, Okinawajima Island: eight unique haplotypes, and Miyako Island: five unique haplotypes), although Ishigaki had only one unique haplotype (Supplementary Figure 4).
Pairwise comparisons of genetic distance for all islands showed significant differences except for between Zamami and Miyako Islands, probably due to small sample sizes (Supplementary Table 4).

DISCUSSION
Based on COI sequences, Uthicke et al. (2010) showed that H. atra includes three clades separated by a relatively large genetic differentiation. Following this, in Okinawa, Takano et al. (2017) confirmed the local presence of clades 1 and 2. In this study, the most common H. atra clade observed in Okinawa was clade 2 ( Figure 5; Uthicke et al., 2010), and we confirmed in our study that this clade is common in the central and southern Ryukyus.
On the other hand, in the current study, clade 3 was very rare, with only a single specimen observed from Zamami Ama. The maximum genetic distance between H. atra specimens for COI sequences obtained in this study (= 4.1%) was close to the maximum distance previously reported between New Caledonia and Great Barrier Reef individuals (4.5%, Uthicke et al., 2010). The phylogenetic tree based on 16S showed no obvious clades. Although 16S has been little used in genetic investigations of H. atra compared to COI, the maximum distance observed in our study (= 2.3%) does not support the existence of cryptic species within H. atra considering the minimum genetic distance (= 4.1%) between two different species in Holothuroidea as observed in a previous study (Kerr et al., 2005). However, further investigations are needed in order to establish if cryptic species truly exist within H. atra. Even though there have been no previous focused studies on the genetic population structure of H. atra inhabiting the Ryukyus, our observed genetic diversities were relatively high compared to those of previous works conducted on the closely related species H. edulis and S. chloronotus across the same region (eight and three 16S haplotypes observed; Soliman et al., 2016a,b). Additionally, despite Skillings et al. (2011)'s results showed that the Okinawan H. atra population had the lowest genetic diversity in a study area spanning central to western Pacific Ocean with only three COI haplotypes, in the current study the actual diversity of this species seems to be much higher, comparable to the diversity observed in Hawai'i (Hd = 0.88 in Hawai'i in Skillings et al., 2011;Hd = 0.77 in Okinawa).
Potential reported causes of genetic diversity decline have included over-harvesting , coastal development (Soliman et al., 2016a), and coastal alteration (Nehemia and Kochzius, 2017). Our study partially supports these previous observations as Chatan and Uruma, which have comparatively well-developed urban coastlines (Masucci and Reimer, 2019), had the smallest number of COI haplotypes. However, in the 16S dataset, Ishigaki St. 27 had the lowest number of haplotypes (one haplotype) followed by Ishigaki St. 20, Chatan, Uruma, Manza and Kayou (two haplotypes). Among these locations, only Chatan and Uruma have developed or artificial coastlines.
In the current study, Zamami Ama had the highest genetic diversity with many unique haplotypes in both the COI and 16S datasets. Similar results have been reported for the scleractinian coral Acropora digitifera Dana, 1846 in Okinawa, in which Zamami showed the highest number of unique haplotypes, particularly when compared with Okinawajima (Nakajima et al., 2010). The Kerama Islands, including Zamami Island, have been a quasi-national park from 1978 and have been included within Kerama-Shoto National Park since 2014. Similarly, Manza, where the highest COI genetic diversity around Okinawajima Island was observed in our study, has been within a quasi-national park since 1972. However, these results were not confirmed by the 16S analyses, in which Manza had only two haplotypes. In Japan, quasi-national and national parks have restrictions put on activities that may potentially disturb the natural environment such as coastal landfilling (land reclamation), sewage discharge, and construction (Ministry of Health Labor Welfare, 1957). High genetic diversity in protected areas have been reported around the world in other sea cucumber species, such as for example Holothuria polii Della Chiaje, 1824 and H. tubulosa Gmelin, 1791 in Turkey . Therefore, we speculate that these environmental protections, as imperfect as they may be (Shinbo, 2016), could be crucial for the local conservation of H. atra's genetic diversity. Although Zamami Island waters are not a complete no-take zone, fisheries rules are much stricter than in completely non-protected areas.
To investigate this further, we compared the geographical regions in our study in terms of fishing pressure. Our results showed that the number of fishers and fishing boats per total coast length was much less in Zamami Village than in both Ishigaki City and Okinawajima Island. Likewise, the annual catch of "others" (including sea cucumbers) per total coastal length was also lower for Zamami Village than Ishigaki City, although the catch of "others" was at the same approximate level for Zamami Village and Okinawajima Island. It should be noted that the definition of "others" is different between Zamami Village and Okinawajima Island, and this may obscure the true amounts of sea cucumber harvests in each region. Another possibility is hand-to-hand trade ("hamauri" in Japanese), which has been reported in Okinawa (Okinawa General Bureau, 2017). Hamauri is an unregulated way of trading that makes it difficult to track the actual amount of harvested and traded organisms. Based on past research and the data available from Okinawa Prefecture regarding fishing pressure, we conclude that the high genetic diversity we observed in H. atra at Zamami is most likely due to the positive impacts of being within a national park, including benefits from reduced fishing pressures. To better confirm this, Okinawa Prefecture should work toward obtaining more accurate sea cucumber catch information.
In pairwise differentiation analyses based on COI sequences, Chatan and Zamami Ama were unique populations. As well, both Ishigaki populations were similar to each other but significantly different from all other populations. Chatan's COI results showed low genetic diversity with only three haplotypes, one of which was dominant (28/30 individuals). Surprisingly, the two Zamami populations (Ama and Port) showed a significant genetic break despite being only one km apart. Likewise, Kin and Uruma, separated by 13 km on Okinawajima Island's east coast, were significantly different. Similar results were previously observed for H. edulis in Okinawa (Soliman et al., 2016a), H. atra in Hawai'i (Skillings et al., 2011), H. scabra in Australia, (Uthicke and Benzie, 2001), and also for H. mammata (Grube, 1840) and H. polii (Borrero-Pérez et al., 2011;Valente et al., 2015) in the Mediterranean Sea. These results indicate that even geographically close sea cucumber populations may retain unique genetic diversity, and therefore each population may be needed to be protected separately.
In addition, another explanation to the low gene flow between adjacent Zamami locations could be caused by high rates of sexual reproduction failure, as previously implied in other locations (Chao et al., 1994). This hypothesis is further supported by the fact that juveniles of Holothuriida species were not observed in previous studies (Chao et al., 1994;Benzie, 2001, 2002;Thorne et al., 2013). However, in the current study, the biological factor(s) determining genetic structure were not investigated in detail. Therefore, ecological surveys assessing the reproductive behavior and recruitment of H. atra at Zamami are needed in the future. At a same time, other DNA markers such as from the nuclear region may provide finer scale resolution of populations and should be investigated in the future.
From the AMOVA analyses, populations were categorized into six groups that were not island-dependent groupings, except for the grouping of the two Ishigaki locations together. Populations around Okinawajima Island were very disjunct, with three groupings (Chatan, Uruma, Sesoko + Kayou + Manza + Kin + Odo) that did not coincide to geographical groupings observed in other previous studies on other marine invertebrates (e.g., Nishikawa and Sakai, 2005). For example, previous studies conducted around Okinawajima Island on two different species of sea cucumbers (Soliman et al., 2016a,b), and on the amphipod Leucothoe vulgaris White and Reimer, 2012 (White et al., 2015), showed clear east and west coast genetic breaks; patterns that we did not observe in the current study. However, a complex disjunct pattern was previously observed in the tidal snail Batillaria flectosiphonata Ozawa, 1996(Kojima et al., 2003 in Okinawajima Island and the authors speculated that the varied genetic composition may be caused by the complex geographical features of Okinawajima Island and/or different dispersal behavior. Based on the results obtained here we found that H. atra populations from Okinawajima Island and those on more distant islands such as Zamami and Ishigaki Islands are genetically isolated from each other. This is consistent with many previous studies, which have confirmed isolation between Okinawajima-Kerama-Sakishima Islands (including Ishigaki) populations in the scleractinian corals Pocillopora damicornis Linnaeus, 1758, Goniastrea aspera Verrill, 1866, and A. digitifera, as well as between Okinawajima and Sakishima Islands in the tideland snail B. flectosiphonata (Adjeroud and Tsuchiya, 1999;Kojima et al., 2003;Nishikawa and Sakai, 2005;Shinzato et al., 2015). These results are somewhat surprising, given that it is often hypothesized that the Kuroshio Current should bring connectivity to coral reef organism populations around different islands, particularly those with planktonic larvae and comparatively high dispersal capabilities (Yorisue et al., 2020). It should be noted that it has been shown in a previous study that H. atra did not have significant genetic differentiation in pairwise comparisons between locations up to 2,000 km distant in Hawai'i (Skillings et al., 2011). Therefore, the main factor(s) that determine H. atra population connectivity in Okinawa and the Kuroshio region remain unknown.
In the COI haplotype network, clade 2 had a star-like shape (dominant central haplotype connected to other haplotypes with few mutations), indicating recent population expansion (Ferreri et al., 2011). A star-like shape haplotype network is one of the common features of many echinoderm species (e.g., the blue sea star Linckia laevigata Linnaeus, 1758 in Otwoma and Kochzius, 2016; tropical sea urchins Acanthaster planci Linnaeus, 1758 and Tripneustes gratilla Linnaeus, 1758in Liggins et al., 2014sea cucumbers H. atra in Skillings et al., 2011;H. polii in Valente et al., 2015) and it is thought to show population expansion after the Last Glacial Maximum, when sea levels were much lower than present day (Ni et al., 2014). However, to better understand the historical fluctuation of echinoderm population sizes, more analyses on other species, combined with wider geographical surveys, should help confirm this theory in the future.
In this study, the Chatan population had a significantly negative Tajima's D value (Tajima's D = −2.23, P-value = 0.001) and Zamami Ama had a significantly negative Fu's FS value (Fu's FS = −5.74, P-value = 0.001). In the case of Zamami Ama, it is thought that this population has been expanding due to the presence of the highest number of unique haplotypes observed at any location, in both COI and 16S datasets. However, from previous studies, it has been inferred that smaller sea cucumber individuals tend to reproduce asexually rather than sexually (Bonham and Held, 1963;Chao et al., 1993Chao et al., , 1994Uthicke, 1997). Another study indicated that nearshore reefs had more asexually reproduced individuals than mid-shelf reefs, likely induced by the high levels of environmental disturbances near shore . Chatan is also a location comparatively influenced by anthropogenic disturbances such as nutrient input from terrestrial runoff (Yang et al., 2013;Mukai et al., 2020). Therefore, Tajima's D test's significantly negative result in the Chatan population might reflect population expansion following a shift to asexual reproduction. However, to validate this hypothesis, further periodic investigations focusing on individual size fluctuations and seasonal changes of the occurrence of fission are needed.
From a conservation point of view, as each island examined here possessed many unique haplotypes as well as clear population genetic structure, we suggest considering each main islands' population as a different MU. Further investigations, including finer-scale genetic analyses, a larger number of loci, and long-term monitoring or observations are required in order to better manage and preserve ecologically important species such as H. atra.

CONCLUSIONS
In this study, we genetically assessed the diversity of H. atra populations at 11 different locations in the Ryukyu Islands, southern Japan. We found that some populations are isolated, despite being geographically very close to each other. From these results, H. atra population diversity in southern Japan needs to be considered at the local scale through the development, in conjunction with local fisheries communities, of effective management and conservations plans. Our results indicate that comparatively better conserved H. atra populations in national and quasi-national parks show higher genetic diversity compared to locations where coastlines have been altered and where less harvesting regulations exist. This implies that even partial restrictions on fishing may help even partially protect the genetic diversity of H. atra.

AUTHOR CONTRIBUTIONS
Research plan was conceived by KH, TS, IF-S, and JDR. All field work performed by KH with help from AP. KH conducted all laboratory experiments and IF-S contributed to sequencing. Data analyses were performed by KH, TS, and JDR. Manuscript was written by KH, TS, AP, IF-S, and JDR. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Drs. Hin Boo Wee and Giovanni Masucci for advice on statistical analyses and Dr. Yee Wah Lau (all MISE, University of the Ryukyus) for checking English. As well, we thank Federico Clementoni (MISE), who helped with sampling. Zamami and Ishigaki Islands' expeditions were conducted with the Scripps Oceanographic Institute's 100 Island Challenge in summer 2019 and we thank Dr. Brian Zgliczynski, Lindsey Bonito, Christopher Sullivan, Sho Kodera, and Gabe Turner. We thank boat captains Takaya Naka at Ishigaki, Kazuya Hayashi at Okinawajima, Masaaki Isa at Miyako and Hisao Miyahira at Zamami. We thank Kazumi Inoha from the Okinawa Graduate School of Science and Technology (OIST) for logistical arrangements during the 100 Island Challenge expedition. We thank two reviewers for their constructive comments.