Genetic Structure of the Goniopora lobata and G. djiboutiensis Species Complex Is Better Explained by Oceanography Than by Morphological Characteristics

Species delimitation of closely related corals is often challenging due to high intraspecies morphological variation and phenotypic plasticity with a lack of characteristic features and scarcity of relevant molecular markers. Goniopora spp. are one such coralline group, and the species status of Goniopora lobata and Goniopora djiboutiensis, an Indian and Pacific Ocean hermatypic coral species complex, has been questioned on the basis of previous molecular and morphological analyses. To further examine the species boundaries between G. lobata and G. djiboutiensis in Japan, specimens collected from areas spanning from Ryukyu Island to temperate Japanese regions were morphologically identified based on traditional morphological descriptions. Then, the genetic structure of the G. lobata and G. djiboutiensis species complex was examined using six newly developed microsatellite markers. The majority of the collected specimens had intermediate morphologies, and a STRUCTURE analysis using the LOCPRIOR model based on typical G. lobata and G. djiboutiensis morphology indicated that there were no genetic differences between these morphologies. On the other hand, STRUCTURE analysis based on oceanographic regions revealed that there was a genetic break between the temperate and subtropical regions. This weak genetic break corresponded with the Kuroshio-associated barrier, which prevents larval transport between subtropical and temperate regions. This study confirms that the current morphological identification criteria for G. lobata and G. djiboutiensis do not match the existing genetic boundaries and thus the two should be regarded as a species complex. This study also highlighted the robustness of using a multi-locus population genetic approach, including a geographic context, to confirm the species boundaries of closely related species.


INTRODUCTION
Species are the fundamental units of biodiversity and are worth consideration in the conservation of ecosystems. However, the delimitation of some closely related taxa, including corals, is often challenging due to the fact that coral skeletons exhibit high phenotypic plasticity under different environmental conditions. The paucity of effective molecular tools for coral identification, together with frequent hybridization and possible reticulate evolution, has long hindered the classification of many reefbuilding corals (Todd, 2008). The delimitation of coral species is essential for conservation because one-third of all coral species are now facing the threat of extinction globally (Carpenter et al., 2008). In Japan, corals in subtropical areas have experienced mass bleaching because of elevated water temperatures, whereas some corals have reportedly migrated to and expanded in temperate areas (Yamano et al., 2011). Therefore, studies on the genetic connectivity between subtropical and temperate areas of coral species are important to understand such phenomena. However, no study has investigated the genetic connectivity of coral species between subtropical and temperate areas, except for Acropora (Nakabayashi et al., 2019). Moreover, no study has focused on other coral species, possibly because the identification of coral species in the field is challenging.
Corals belonging to the genus Goniopora are widely distributed across the Indian and Pacific Oceans, including subtropical and temperate areas in Japan, and have few skeletal characteristics in addition to highly variable skeletal and polyp morphologies. There are 15 nominal species in the genus Goniopora found in Japan (Nishihira and Veron, 1995). Veron (1993) reported that one of the relatively common species, Goniopora lobata, is distributed in both subtropical and temperate areas (up to Tatayama) and is very difficult to identify accurately, especially in temperate areas where skeletal variations are primarily correlated with non-reef environments. Kitano et al. (2013) reported the existence of many intermediate morphologies between closely related G. lobata and G. djiboutiensis species in Japan. G. djiboutiensis is less common than G. lobata, and its distribution overlaps with that of G. lobata, especially in subtropical areas (Veron, 1993;Nishihira and Veron, 1995). The first comprehensive genetic approach to delineate the species covering this whole genus was conducted in the Red Sea (Terraneo et al., 2016). The authors of the aforementioned study found discordance among genetic, morphological, and micromorphological data related to G. lobata and G. djiboutiensis, wherein G. lobata and G. djiboutiensis shared common alleles and thus could not be delineated using genetic markers including a mitochondrial region, the internal transcribed spacer (ITS) region, and the calmodulin (CalM) and ATPaseβ genes (ATPsβ). Based on this molecular and morphological evidence, the aforementioned authors questioned the species status of G. lobata and G. djiboutiensis, although the study was ultimately inconclusive.
Although species delimitation analysis is often attempted via a phylogenetic approach, population genetic approaches considering the geographic contexts of processes driving population divergence and speciation have proven to be more effective in identification of species boundaries in some phylogenetically complex groups (Medrano et al., 2014;Nakabayashi et al., 2019). Thus, further genetic studies using a multi-locus population genetic approach may help to confirm the species boundary between G. lobata and G. djiboutiensis. In this study, G. lobata-like and G. djiboutiensis-like specimens were collected from a wide range of coral habitats in Japan, including different oceanographic regions (subtropical and temperate regions) separated by an ocean current system. Six novel polymorphic microsatellite markers for Goniopora spp. were developed to analyze the genetic structures of the G. lobata and G. djiboutiensis species complex. Then, the specimens were comprehensively studied based on traditional morphological descriptions, ITS sequences, and microsatellites to reveal the species boundaries between G. lobata and G. djiboutiensis.

Sample Collection
Owing to the fact that species identification of living colonies is difficult in the field, G. lobata-and G. djiboutiensis-like samples were collected from 15 areas spanning from Ryukyu Island to temperate regions of Japan (Table 1 and Figure 1). Samples (100-500 cm 3 ) were collected using a hammer and chisel. A small portion of each coral sample (<1 cm 3 ) was preserved in CHAOS solution (Fukami et al., 2004) to dissolve proteins for DNA analyses, and the remnants of the samples were bleached for morphological analyses. Total genomic DNA was extracted from colonies of the samples using the method outlined by Fukami et al. (2004). The skeletal specimens are retained either in University of Miyazaki (MUFS) or in Seto Marine Biological Laboratory (SMBL).

Morphological Identification
Morphological identification was based on original species descriptions and related literature (Milne-Edwards and Haime, 1851;Bedot, 1907;Vaughan, 1907;Veron and Pichon, 1982;Terraneo et al., 2016). Each five mature corallites per colony were measured to determine the range of calice diameter using a vernier caliper. Seven skeletal characteristics were also examined (calice depth, width of calice wall, morphology of columella, columella divided to six parts or not, tips of columellar threads, existence of pali or paliform lobes, and septa under binocular stereomicroscope to classify the species. If specimens had morphological characteristics typical of both G. lobata and G. djiboutiensis in different categories, they were classified as "intermediate."

ITS Analysis
In order to compare the microsatellite results of the present study with those obtained in previous studies, the ITS region (including partial 18S, ITS-1, 5.8S, ITS-2, and partial 28S) was sequenced and compared with the sequences reported by Terraneo et al. (2016). Polymerase chain reaction (PCR) and sequencing conditions were exactly the same as in protocol (Kitano et al., 2013). Previously published sequences (G. lobata and G. djiboutiensis in the Red Sea; Supplementary Table 2) were also included in the following analysis to identify the conformity of the resultant genetic clades with those reported by Terraneo et al. (2016). The ITS sequences with a maximum length of 926 bp were aligned using MAFFT online (Katoh et al., 2017), and all gap regions were trimmed using the "-nogaps" option in trimAl v. 1.4 (Capella-Gutiérrez et al., 2009), which removes all gap sites in the alignment, resulting in 899 bp in length for analysis. To identify the genetic clades defined by Terraneo et al. (2016), a phylogenetic tree was constructed using two different algorithms: maximum likelihood (ML) using IQ-TREE2 v.2.0.6 (Minh et al., 2020) and Neighbor-Joining (NJ) using MEGA X (Kumar et al., 2018). In total, 135 sequences across 11 Goniopora spp. were used for the analysis, including 44 reference sequences reported by Terraneo et al. (2016). For the ML tree, the best-fitting nucleotide substitution model (K2P + R2) was estimated based on the Bayesian information criterion (BIC), using the "-m MFP" option, which performs extended model selection followed by tree inference. Confidence values for the trees were inferred using the ultrafast bootstrap method (Hoang et al., 2018) with 1,000 replicates. For the NJ tree, the Kimura two-parameter (Kimura, 1980)

Microsatellite Marker Development Using High-Throughput Sequencing
The total genomic DNA for microsatellite marker isolation was obtained from the samples listed in Table 1. The 24 genomic DNA isolates from CHAOS solution samples were mixed and were further purified using the DNeasy Blood and Tissue kits (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The purified genomic DNA was fragmented into c.a. 300-800 bp, and the adapter sequences were ligated to each end of the fragments. The resulting DNA library was subjected to shotgun sequencing on GS FLX Titanium picotiter plates using the Roche GS FLX 454 system. A total of 603,951 sequences were obtained with a total of 306 Mb. These sequencing reads were assembled using Newbler Assembler software ver. 2.6 (Roche Applied Sciences, Indianapolis, IN, United States) with default settings. Subsequent contigs and singletons were screened using a Perl pipeline coupled with Tandem Repeats Finder ver. 4.0.4 (Benson, 1999) and PRIMER 3 ver. 2.2.2 beta (Rozen and Skaletsky, 2000), which were eventually used to design primers for the potential microsatellite loci (Nakamura et al., 2013). A total of 1,136 pairs of primers with di-, tri-, and tetra-nucleotide motif loci were designed, from which 240 primer sets were selected for initial amplification trials. After removing the primers that did not amplify correctly (multiple bands or no amplification), coral microsatellite markers were further discriminated from those of their symbiont algae (family: Symbiodiniaceae) using the same primers. Owing to the fact that Goniopora spp. is known to host Symbiodiniaceae species of Cladocopium goreaui (formerly clade C type C1) in Taiwan (Chen et al., 2005) and Hong Kong (Wong et al., 2016) and Durusdinium trenchii (formerly clade D type D1-4) in Indian Ocean (LaJeunesse et al., 2010), the following four Symbiodiniaceae culture strains were tested: AJIS2-C2 (reported by Yamashita and Koike, 2013, Symbiodinium microadriaticum former clade A type A1), CCMP2466 (C. goreaui), CCMP2556 (D. trenchii), and AJIS5-A10 (D. trenchii, previously isolated by one of the authors H. Yamashita). The culture strains CCMP2466 and CCMP2556 were purchased from the Provasoli-Guillard National Center for Culture of Marine Algae and Microbiota (ME, United States). As a result, total 39 primer pairs of potential microsatellite markers including trinucleotide motif loci for Goniopora spp. were selected. To examine the effectiveness of PCR amplification of the 39 fluorescent-labeled primer pairs that were successfully identified as Goniopora microsatellite regions, we performed PCR in a thermal cycler (PC-808, ASTEC) in a reaction mixture (10 µL) containing 1 µL (ca. 10 ng) template DNA, 5 µL of 2× Go Taq Colorless Master mix (Promega Corporation, United States, containing Taq DNA polymerase, dNTPs, MgCl 2 , and reaction buffers at optimal concentrations), 0.35 µM of each designed primer pair, with one fluorescently labeled primer, and 3.76 µL of ultra-pure quality water. The PCR cycling conditions were as follows: 2 min at 95 • C, 30 cycles at 95 • C for 30 s, 50 • C-55 • C (primer-specific annealing temperature, Table 3) for 30 s, 72 • C for 30 s, and a final elongation at 72 • C for 5 min.
Genetic Structure of the G. lobata and G. djiboutiensis Species Complex To examine the species boundaries and spatial genetic structure of the G. lobata and G. djiboutiensis complex, 76 samples were genotyped (Table 1 and Figure 1). One microliter of the 80times diluted PCR product was added to 10 µL Hi-Di formamide containing 0.15 µL GeneScan 600 LIZ size standard (Applied Biosystems), followed by denaturation at 95 • C for 3 min. The samples were run on an ABI 3130 xl genetic analyzer (Applied Biosystems) and were analyzed using GeneMapper ver. 5.0 software (Applied Biosystems). Using Amakusa, Shirahama, and Miyazaki samples from which more than 10 samples could be collected (Table 1), the observed and expected heterozygosity and the Hardy-Weinberg equilibrium (HWE) were calculated using GenAlEx ver. 6.5 (Peakall and Smouse, 2012). Linkage disequilibria between the loci were tested in Genepop ver.4.0 (Raymond and Rousset, 1995) with 5,000 dememorizations, 500 batches, and 5,000 iterations per batch. Micro-Checker version 2.2.3 (Van Oosterhout, 2004) was used to check the presence of scoring error due to stuttering, large allele dropout, and null alleles for each locus.
A total of 10 independent runs of Structure ver. 2.3.4 (Pritchard et al., 2000) were calculated assuming K = 1-6 with 500,000 Markov chain Monte Carlo (MCMC) iterations after a 500,000 burn-in using admixture and an allele frequency correlated model to examine genetic structure. To examine the genetic structure assuming different hypotheses, three different priors were used as follows: (1) no prior, (2) morphological prior (two categories: typical G. lobata-morph and typical G. djiboutiensis-morph) using LOCPRIOR (Hubisz et al., 2009), and (3) oceanographic regions prior (two categories: temperate and subtropical regions). While it was assumed that the largest mean LnP(K) with a biologically meaningful low number of K was the optimal number of clusters [as suggested by Pritchard et al., 2000], delta K (Evanno et al., 2005) was also estimated using STRUCTURE HARVESTER (Earl and von Holdt, 2012). Then, CLUMPAK (Kopelman et al., 2015) was used to visualize the STRUCTURE results.
Additionally, the population-based genetic structure was visualized by principal coordinate analysis (PCoA) using GenAlEx ver. 6.5 (Peakall and Smouse, 2012). Furthermore, Arlequin ver. 3.5 (Excoffier and Lischer, 2010) was used to examine the genetic structure (F CT ) between temperate (Ooura and Iriomote) and subtropical (the others) regions and the calculated global F ST . With regard to these population-based analyses, the Kikai population was excluded due to an extremely small (n = 3) sample size.

RESULTS
As mentioned by Kitano et al. (2013), most of the collected samples (64 out of 76) had intermediate morphologies, thereby indicating that it was impossible to distinguish between the two species groups using the current morphological criteria (Figure 2 and Supplementary Table 1). None of the eight morphological features previously used were effectively able to discriminate between G. lobata and G. djiboutiensis (Supplementary Table 1).
The ML and NJ trees using the ITS region were consistent; therefore, only the ML tree has been shown in Figure 3 (raw data available in Supplementary Data 1). The reconstructed ITS tree suggested that all the samples used in this study belonged to clade II in Terraneo et al. (2016), which was previously found to be a G. lobata and G. djiboutiensis complex clade (Figure 3). Consistent with the findings reported by Terraneo et al. (2016), there was no subclade corresponding to geographic areas or specific morphological features within clade II.
After PCR testing, six primer pairs of Goniopora spp. that could amplify microsatellite regions were identified, including the tri-nucleotide motif ( Table 2). No significant linkage disequilibriums were detected among the loci after Bonferroni correction (α > 0.05). Only the G329 locus in the Miyazaki population was found to be significantly deviated from the HWE locus (P = 0.001, Table 3). Additionally, Micro-Checker analysis detected that the G329 locus might include null alleles. Although analysis was completed with and without G329 and it was confirmed that there was no discrepancy in the results, only the results excluding G329 have been presented this paper. All the used genotyping data is available in Supplementary Table 2. STRUCTURE analysis with a uniform prior showed admixed patterns, and the maximum mean LnP(K) was observed at K = 1, thereby indicating that only one cluster was estimated (Figure 4  Table 3). A similar result was obtained using a morphological prior, and while K = 2 was selected as the optimal number of clusters based on the largest Mean LnP(K) value, no genetic structuring between the typical G. djiboutiensis-morph and G. lobata-morph was found (Figure 4). On the other hand, STRUCTURE analysis using oceanographic regions suggested that K = 2 was the best model based on Mean LnP(K) as well as delta K estimation (Supplementary Table 3), and a clear genetic break between subtropical and temperate regions was detected with a low mean r value (0.365), thereby indicating that the geographic prior was able to detect a genetic structure within the data (Figure 4 and Supplementary Table 4).
Principle coordinate analysis indicated that the subtropical populations (Oura and Iriomote) were separated from the six temperate populations by coordinate 1, which explained 65.12% of the variation in the data (Figure 5). A significant genetic structure was observed according to the global F ST value (F ST = 0.033, P = 0.027). Genetic structuring among predefined oceanographic groups (subtropical vs. temperate) was also significant (F CT = 0.095, P = 0.037).

Microsatellite Markers as a Tool to Examine Species Boundaries
Recent advances in high-throughput sequencing technology have made microsatellite markers rather outdated. In general, a higher number of loci can resolve finer scale genetic structures (Jeffries et al., 2016). However, in some circumstances, using microsatellite markers can be advantageous compared to newly emerged genomic methods such as RAD-seq (Wagner et al., 2013). Hodel et al. (2016) reviewed and compared microsatellite markers with more recently developed genomic methods and concluded that microsatellite markers remained useful in terms of both effectiveness and cost-efficiency (Hodel et al., 2016).
Microsatellite markers are useful in many molecular ecological applications, including parentage analysis, conservation genetics, phylogeography, and population genetics (Hodel et al., 2016). Among the many molecular markers, microsatellite markers are more useful in detection of species boundaries between closely related coral species. Previous population genetic studies using microsatellite markers in corals with appropriate sampling design including geographical contexts have revealed underestimated species diversity even when phylogenetic trees cannot detect recent speciation owing to incomplete lineage sorting (Bongaerts et al., 2011;Nakajima et al., 2012;Yasuda et al., 2015).
Microsatellite analyses using sufficient sample sizes and including a geographic context enable us to detect species boundaries by accounting for incomplete lineage sorting and/or introgression (Nakajima et al., 2012;Nakabayashi et al., 2019). However, caution is required when planning the sampling design in terms of both the number of samples and the sampling strategy itself because individual-based Bayesian clustering methods, including STRUCTURE, are known to be sensitive to sampling design (Schwartz and McKelvey, 2008). In addition, because the molecular evolutionary patterns of microsatellite regions are not as well-known and the "phylo-" part of the analysis FIGURE 3 | Unrooted Maximum likelihood phylogenetic tree using internal transcribed spacer (ITS) regions including 44 reference sequences (Goniopora somaliensis, G. savignyi, G. djiboutiensis, G. lobata, G. stokesi, G. albiconus, G. tenuides, G. pedunculata (G. minor in Terraneo et al., 2016), G. gracilis, G. pendulus, and G. norfolkensisis) taken from Terraneo et al. (2016). Clade names also correspond with those reported by Terraneo et al. (2016). Nodes in red colors represent the Japanese samples used in the present study. The values on each node represent bootstrap values for the ML and NJ methods. Bootstrap values from the ML method are shown in boldface. Scale bar shows 0.005 substitutions par site.   will be lost in population-based analysis, integrated approaches taking phylogenetics, morphology, and ecology into account are nonetheless important. This study confirms the absence of species boundaries between G. lobata and G. djiboutiensis morphologies using microsatellite markers. Such a multi-locus population genetic approach, including a geographic context, would be also useful for delineating the species boundaries of other closely related species including corals.

Species Status of G. lobata and G. djiboutiensis
Although the concept of "species" is a fundamental unit for biodiversity conservation, the species delimitation of closely related coral species has been morphologically and genetically challenging due to high phenotypic plasticity and the paucity of suitable molecular markers. There are many proposed and widely used species concepts, including morphological species concepts that require morphological differentiation, the genealogical species concept that requires allelic coalescence (Baum and Donoghue, 1995), and the biological species concept that requires reproductive isolation (Mayr, 1942). However, whichever species concept is used, a species represents an independent metapopulation lineage through time (de Queiroz, 2005).
In this study, the species boundary between G. lobata and G. djiboutiensis was examined based on previously used morphological characteristics as well as newly developed microsatellite markers, including a geographic context for testing the hypothesis of "independent metapopulation lineages through time" for each species. The majority (84.2%) of the samples collected in Japan had intermediate forms, suggesting that the eight morphological characteristics used in the present study could not distinguish G. lobata from G. djiboutiensis. This fact corresponds with a previous study that reported an abundance of morphotypes in Goniopora spp. including G. lobata and G. djiboutiensis in the region (Kitano et al., 2013) and Veron (1993) who noted that G. lobata in Japan are "always difficult to identify with certainty in higher latitudes, where skeletal variations occur that are primarily correlated with nonreefal environments." Additionally, microsatellite STRUCTURE analysis using either uniform prior or morphological prior indicated genetic homogeneity and could not detect any genetic boundaries between typical G. lobata and G. djiboutiensis morphologies, which was consistent with previous studies using other molecular markers (Kitano et al., 2013). Lastly, a genetic break (0.030 of allele frequency divergence between pops, F CTlike value estimated by STRUCTURE) was discovered using a geographic prior, which could be explained by the presence of population isolation caused by oceanographic conditions rather than by morphological characteristics. This suggests that the current morphological characteristics defining G. lobata and G. djiboutiensis do not represent independent metapopulation lineages over time. Goniopora lobata and G. djiboutiensis could not be morphologically and genetically delimited as species. It is likely that the high morphological diversity in the G. lobata species complex may lead to overestimation of morphological species. The present study underlines the usefulness of a population genetics approach including a geographic context in examining the boundaries between closely related species. Given that the nomenclature of Goniopora lobata Milne-Edwards and Haime, 1851, has been used more than that of Goniopora djiboutiensis Vaughan, 1907, it would be appropriate to use "G. lobata species complex" rather than "G. djiboutiensis species complex" to describe the samples in clade II.

The Kuroshio-Associated Oceanographic Barrier
Remarkably, a clear genetic break was discovered between subtropical and temperate regions, irrespective of morphological characteristics, suggesting that genetic exchange via larval dispersal was somewhat limited between subtropical and temperate regions in the G. lobata species complex. This genetic break corresponds not only to climatic regions (subtropical vs. temperate) but also to the Kuroshio-associated larval dispersal barrier reported by Nakabayashi et al. (2019). Thus, it is possible that this genetic barrier could be caused by oceanographic barriers and/or local adaptation owing to environmental factors, such as water temperature.
According to larval dispersal oceanographic simulations, larval transport between subtropical and temperate regions in Japan via the Kuroshio Current and its associated currents is limited to a few generations, and this phenomenon is referred to as the Kuroshio-associated barrier (Nakabayashi et al., 2019). Goniopora spp. are known to be gonochoric and are broadcast spawners (Richmond and Hunter, 1990) with at least 5 days of pelagic larval duration (Babcock and Heyward, 1986). Although the dispersal potential of the G. lobata species complex is not FIGURE 4 | STRUCTURE bar-plot results using no prior (upper) at K = 3, geographic prior at K = 2 (middle), and morphological prior at K = 2 (lower). The x-axis indicates each individual, and the y-axis indicates the proportion of assigned clusters shown in different colors.
FIGURE 5 | Population-based principal coordinate analysis. Coordinate 1 explained 65.12% of the variation in the data and coordinate 2 explained 21.45% of the variation. Blue rhomboids indicate temperate populations and red rhomboids indicate subtropical populations. The name of the sites corresponding to each abbreviation is found in Table 1. precisely known, considering the markedly lower population density of the G. lobata species complex than that of Acropora spp. especially in subtropical areas, the number of larvae that can cross the Kuroshio-associated barrier would be fewer than that of Acropora spp., which have similar genetic breaks.
Alternatively, it is possible that this genetic break is caused by environmental factors. Among the environmental factors that we summarized for each sampling region [the mean surface water temperature of the coldest month, yearly average of mean surface water temperature, chlorophyll-a concentration, and particulate inorganic carbon (PIC) concentration], (Yara et al., 2012;NASA, 2018) the range of chlorophyll-a and PIC concentration overlapped between temperate and subtropical regions (Supplementary Table 4). The ranges of water temperature between the two regions did not overlap, although the difference between the two regions was minor (Supplementary Table 4). In addition, the following gradual changes along latitude were observed: 20.82-22.87 • C (coldest month average) and 24.63-26.17 • C (year average) in the subtropical region and 14.23-19.57 • C (coldest month) and 20.27-24.22 • C (year average) in the temperate region (Yara et al., 2012;Supplementary Table 4). If there were no clear physical barrier between the temperate and subtropical regions, one would expect a gradual change in different genotypes along the Kuroshio Current. The clear genetic break between temperate and subtropical regions is, therefore, impossible to explain only by the differences in environmental factors (i.e., water temperature, PIC concentration, and chlorophyll-a concentration) alone, although we cannot deny the possibility that these differences in environmental factors further facilitate local adaptation and genetic divergence.
In the present study, we used putatively neutral microsatellite markers. Thus, further genomic studies that can detect loci under natural selection, such as restriction site-associated DNA sequencing (RAD-seq, Wagner et al., 2013), as well as physiological studies of larval and adult corals between temperate and subtropical regions are needed to test this hypothesis. Similar genetic evidence of peripheral temperate isolation has been reported in other coral species, such as Acropora hyacinthus (Nakabayashi et al., 2019), and other marine species, including sea slugs (Hoang et al., 2018) and groupers. Further genetic connectivity studies of coral and other marine species with different dispersal potential would provide insights into the formation of genetic breaks associated with the Kuroshioassociated barrier.

Limitations of This Study
While this study indicated a clear genetic break between temperate and subtropical regions irrespective of morphological features, there are some limitations in this study. The number of microsatellite loci analyzed in this study is limited. We only succeeded in isolating six independent microsatellite loci, which may be insufficient to accurately identify asexually produced clones within some populations. Fortunately, to our knowledge there is no report on asexually produced clones in the G. lobata species complex. We have also confirmed that there are no possible clones (having exactly the same genotype based on six loci) within our data, indicating all the samples used in this study were sexually reproduced genetically unique individuals.
A second limitation is the number of samples we analyzed in this study. Although this is characteristic of our target region, 84.2% of studied specimens had intermediate morphologies, and only 13 specimens could be identified as morphotypes of G. lobata (n = 7) or G. djiboutiensis (n = 6). While Terraneo et al. (2016) found major discordance between morphology and molecular data between G. lobata and G. djiboutiensis morphotypes, concordant with our study, they could at least identify G. lobata and G. djiboutiensis morphotypes in the Red Sea (a type locality of G. lobata). Thus, genetic analysis using more G. lobata and G. djiboutiensis morphotypes from other regions including type localities would be necessary to clarify the species boundaries of the G. lobata species complex.
Lastly, this study lacks a morphometric analysis, and individual-based environmental data. Further study, providing more detailed morphometric analysis, together with environmental data for each specimen would offer insights into the morphological variations and phenotypic plasticity of this group.
Using the present data as a baseline, genome-wide genetic analyses such as RAD-seq or MIG-seq using a larger number of collected morphotype specimens, including geographic context, would be useful to further clarify species boundaries and possible ecological speciation in relation to morphological features and evolutionary history of G. lobata species complex.

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: DNA Data Bank of Japan (accession numbers listed in Supplementary

AUTHOR CONTRIBUTIONS
NY and YK conceived the study. YK, TM, and HY collected the samples. YK conducted the morphological identification and Sanger sequencing. NY, SN, HY, and HT developed microsatellite markers. NY conducted the genetic analysis and drafted the manuscript. All the authors have checked and agreed to the final manuscript. governments and worked in agreement with local fishing corporations. We would like to thank Editage (www.editage.com) for English language editing. We appreciate the reviewers and editors who assisted in improving this manuscript.