A total of 1287 Heliopora samples from a total of 73 sites across three western Pacific regions (Japan, Taiwan and Guam) and Singapore, and three Indian Ocean regions (northwestern Australia, Thailand, and the Maldives) were collected using SCUBA from 2008 to 2019 (Supplementary Table 1). Additionally, we used previously collected samples from Japan and northwestern Australia (Supplementary Table 1). The samples include those from the northernmost distribution site (Yaku Island in Japan; Nakabayashi et al., 2017), and the known deepest depth (from 50 m depth from Green Island, Taiwan). All coral fragments were preserved in 99.5% EtOH and kept at −20°C until DNA extraction. During sampling, colonies were photographed, and GPS information was acquired wherever possible. All overseas sample collection was conducted through legal procedures and collection cooperation with the support of local research collaborating institutions; Academia Sinica in Taiwan, University of Guam in Guam, National University of Singapore in Singapore, Western Australian Museum and Curtin University in Western Australia, Phuket Marine Biological Center in Thailand, University of Milano-Bicocca and Marine Research and High Education Center in Maldives (see also FIELD STUDY PERMISSION). Genomic DNA was extracted from the preserved samples using the hot alkaline solution method (Meeker et al., 2007) or QIAGEN DNeasy Blood & Tissue kit (QIAGEN, Hilden, Germany) on site. The first PCR procedure of MIG-seq analysis was carried out at each overseas laboratory except for Thailand samples from which ethanol preserved coral fragments were imported in 2011 (CITES number: AC.0510.2/407). We used them in this study under an additional permit from the Department of Marine and Coastal Resources in Thailand.

All the DNA extraction and first PCR procedure except for domestic samples were carried out locally, and the first PCR products were transported to Miyazaki, Japan. The first PCR step was performed by using 8 pairs of multiplex ISSR primers of Suyama and Matsuki (2015). The fragments were amplified with the Multiplex PCR Assay Kit Ver.2 using a total volume of 7 μL reaction in a thermal cycler with the following modified profile: 94°C for 1 min followed by 29 cycles at 94°C for 30 s, 38°C for 1 min, 72°C for 1 min, and a final extension at 72°C for 10 min. The 50-fold diluted first PCR product was used as the template DNA for the second, tailed-PCR to add the individual index and the Illumina adapter sequence to each sample. PrimeSTAR GXL DNA polymerase (TaKaRa Bio Inc., Otsu, Shiga, Japan) was used for second PCR with a total volume of 12 μL reaction in a thermal cycler with a following profile: 20 cycles at 98°C for 10 s, 54°C for 15 s, 68°C for 1 min. The second PCR products were pooled as a single mixture library by using 1μL of each product. The mixed sample was electrophoresed on 0.1% agarose gel to verify amplification. The size ranged from 350 to 800 bp and DNA was manually extracted from the gel on a transilluminator. After extracting genomic DNA from the gel using FastGene Gel/PCR Extraction Kit (Nippon Genetics, Tokyo, Japan), library concentration was measured with a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, United States). Finally, the DNA library was sequenced on a MiSeq Sequencer (sequencing control software v2.0.12, Illumina, San Diego, CA, United States) using a MiSeq Reagent Kit v3 (150 cycles; Illumina). Image analysis and base calling were performed using real-time analysis software v1.17.21 (Illumina). DarkCycle option was changed “Amplicon-dark 17-3” to “Amplicon-dark 17-17” on the “Chemistry” line (see also Suyama and Matsuki, 2015).

The FASTX-toolkit version 0.0.14 (Gordon and Hannon, 2012), with a fastq-quality-filter setting of –Q 33 –q 30 –p 40, was used to eliminate low-quality reads and primer sequences from the MIG-seq raw data. Adapter sequences for Illmina MiSeq were removed from both 5’ end and 3’ end using Cutadapt version 2.10 (Martin, 2011). Short reads less than 80 bp were excluded using an in-house python script. To further exclude contamination from Symbiodiniaceae-DNA, and to obtain larger numbers of loci (Shafer et al., 2017), the filtered raw reads were mapped onto the reference genome of Heliopora coerulea obtained from Symbiodiniaceae-free larvae (10.6084/m9.figshare.14356418) using Stacks version 2.2.0 (Catchen et al., 2011; Rochette et al., 2019) with the reference-aligned pipeline for single nucleotide polymorphism (SNP) discovery and genotyping. As a first step, the unpaired reads were removed using the repair.sh software tool within BBtools software package (Bushnell, 2017). Next, the retained reads were aligned to the reference genome using Burrows-Wheeler Aligner (BWA) version 0.0.17-r118 (bwa index and mem with default parameters) (Li and Durbin, 2009). Alignments were then compressed, sorted, and indexed with SAMtools version 1.10 (Li et al., 2009). Finally, SNP genotyping for each individual was carried out using the gstacks program in Stacks with default setting. SNP calling for each genetic analysis was carried out using the populations program in Stacks with the option -r to obtain different minimum proportions of genotyped individuals per locus (0.1 for phylogenetic analyses, 0.9 for assigning clones and population genetic analyses), and the option –ordered_export to ensure that only one of the overlapping SNPs is output from reference aligned data. All the output files from Stacks were converted to the appropriate formats for each genetic analysis using PGDSpider version 2.1.1.5 (Lischer and Excoffier, 2012).

We used the software GenoDive version 3.04 (Meirmans, 2020) to assign possible clones within each population based on the infinite allele model (IAM) omitting all missing data. Less than 12 differences of multi-locus genotype within the same population were identified as possible clone mates (estimated by inspection of the pairwise distance histogram in GenoDive). All clone mates except for one individual with the least missing data were removed from this study. After removing possible clones, we selected up to three least missing data samples per sampling site. In case multiple different gross morphologies and/or genetically different lineages were found within the same site, up to three samples per morphology were selected. Finally, we rerun Stacks using this selected dataset excluding possible clones to obtain new SNPs for phylogenetic and population genetic analyses.

We estimated phylogenetic trees based on the Maximum likelihood (ML) method. The program IQ-TREE2 version 2.0.6 (Minh et al., 2020), using the TVM + F + R5 model selected based on the Bayesian Information Criterion (BIC), was used to reconstruct a phylogenetic hypothesis for Heliopora spp. The bootstrap analysis was performed with 1000 replicates using UFBoot, Ultrafast Bootstrap Approximation (Minh et al., 2013). The final tree was drawn using FigTree version 1.4.4 (Rambaut, 2012).

We conducted Principal Coordinates Analysis (PCoA) to examine genetic relationships among Heliopora individuals using GenAlEx version 6.502 (Peakall and Smouse, 2012). In addition, Discriminant Analysis of Principal Components (DAPC) (adegenet package version 2.1.3 in R version 4.0.2) (Jombart, 2008; Jombart et al., 2010; Jombart and Ahmed, 2011) was performed to visualize the genetic structure. The lowest BIC was used to detect the optimal number of K clusters.

In total 339,818,358 raw reads with an average of 264,039 reads per sample were obtained for 1,287 individuals. After filtering low-quality reads, 335,164,915 reads with an average of 260,423 reads per sample were obtained. After excluding index, adapter, and short reads less than 80 bp, 139,837,592 reads with an average of 108,654 reads were obtained. In total 67,119,478 reads were mapped onto the reference-genome with an average of 52,152 reads per sample. We used 795 SNPs across 1,287 individuals (r = 0.9) to find possible clones, and then excluded 576 individuals from the subsequently genetic analyses. We finally selected 245 individuals from 73 sampling sites from 7 regions (DRA Accession No. DRA012077) for phylogenetic (ML tree; 24,741 SNPs, r = 0.1) and population genetic analyses (PCoA and DAPC; 1,092 SNPs, r = 0.9). For the latter, we selected 238 individuals from the final dataset after excluding seven individuals with a missing rate of more than 50% per 1,092 SNPs.

ML phylogenetic analysis revealed three major groups of Heliopora spp. in the Indo-West Pacific with 100% bootstrap values (Figure 1). The phylogenetic hypothesis presented here indicated that there is one new undescribed group (hereafter Heliopora sp. group) genetically different from the other two groups, each of which includes typical morphologies of the previously described species (slender-branching growth form H. hiberniana whose type locality is in the northwestern Australia; lobate growth form H. coerulea whose type locality is in the Indian Ocean). The “Heliopora sp.” group (shown in blue in Figure 1) comprises the HC-A lineage found in Kuroshio Current regions including Japan and Taiwan (Yasuda et al., 2010, 2014, 2015; Taninaka et al., 2019), and Singapore. The “H. hiberniana” group (shown in green in Figure 1) includes the holotype of H. hiberniana found in northwestern Australia (Richards et al., 2018) and the eastern Indian Ocean samples from Thailand and the Maldives. The “H. coerulea” group (shown in red in Figure 1) comprises mainly lobate growth morphs HC-B lineage found in Kuroshio Current regions including Japan and Taiwan (Yasuda et al., 2010, 2014, 2015; Taninaka et al., 2019), Guam, and H. coerulea found in northwestern Australia (Richards et al., 2018). In addition, there were genetically isolated subclades within each group that clustered together by different growth forms or geographic regions (Figure 1; sp.1–3, hib.1–2, and coe.1–3). The Heliopora sp. group was divided into three subclades: Singapore subclade (sp.1), HC-A with only columnar growth form subclade of Japan (sp.2), and HC-A with both encrusting and columnar growth forms subclade found in both Japan and Taiwan (sp.3). Heliopora sp.3 also included a sample collected from 50m depth in Green Island, Taiwan. A sample collected from Yaku Island in Japan, the northernmost distribution site, looks like an outgroup of sp.2 and sp.3 in the phylogeny, possibly because it had a hybrid genotype between Heliopora sp.2 and sp.3 based on STRUCTURE analysis (data shown in Supplementary Figure 1). The group “H. hiberniana” was divided into two subclades: northwestern Australia H. hiberniana subclade (hib.1) and a subclade consisting of Thailand and the Maldives samples (hib.2). The group “H. coerulea” was divided into three subclades; a Guam subclade (coe.1), H. coerulea in northwestern Australia subclade (coe.2), and a subclade consisting of HC-B found in Japan and Taiwan samples (coe.3).

Consistent with the phylogenetic analysis result, population genetic analyses revealed three genetic isolated groups of Heliopora with eight subclades in the Indo-West Pacific (Figures 2A,B). The PCoA showed three clusters based on 11.73% (PC1) and 7.3% (PC2) of the total variability (Figure 2A). The DAPC also showed clear genetic separation of the three clusters and eight subclades (Figure 2B). K = 3 was selected supporting the ML phylogeny result, while K = 8 was determined as the optimal number of clusters by calculating the first 200 PCs for the maximum of 50 clusters based on BIC. The posterior DAPC assignments of both K = 3 and 8 were consistent with the prior clusters (all blue crosses were on red rectangles).

The three groups of Heliopora indicated uneven geographical distribution (Figure 3A). The distribution of Heliopora sp. group was biased to western Pacific Ocean side including Kuroshio Current region (Japan and Taiwan) and the South China Sea (Singapore). On the contrary, the distribution of H. hiberniana group was mainly in the Indian Ocean side such as Kimberley region and Christmas Island in northwestern Australia, the Andaman Sea (Thailand), and the Maldives. The H. coerulea group was distributed both in the Indian Ocean and western Pacific ranging from the Guam and the Kuroshio Current region (Japan and Taiwan) to Kimberley and Christmas Island in northwestern Australia. The distribution of H. coerulea group was partly overlapping with the other two species groups in Japan, Taiwan, and northwestern Australia (Figure 3A, coe.2 and hib.1; coe.3 and sp.2-3).

The gross morphologies of the three groups of Heliopora were highly variable (Figure 3B). H. coerulea group included two typical morphs; lobate growth form (Figure 3B, a–d) and thick-branching growth form (Figure 3B, e). H. hiberniana group included two typical morphs; slender-branching growth form (Figure 3B, f,g) and lobate growth form (Figure 3B, h–j). Heliopora sp. group included two typical morphs; columnar growth form (Figure 3B, k–m) and encrusting growth form (Figure 3B, n,o). The typical growth forms of Heliopora sp. group were relatively distinct from other two species, while the morphologies of H. hiberniana group and H. coerulea group partly overlapped and were difficult to identify by gross morphology.

Phylogenetic and population genetic analyses using genome-wide SNPs by MIG-seq showed that the Indo-West Pacific Heliopora was divided into three groups: H. coerulea group, distributed in the Indo-West Pacific and including HC-B in the Kuroshio Current region; H. hiberniana group distributed mainly in the Indian Ocean including the type locality of H. hiberniana in northwestern Australia; and a new undescribed Heliopora sp. group distributed mainly in the Western Pacific Ocean, including HC-A in the Kuroshio Current region (Figures 1, 2 and 3A; Supplementary Table 1).

The new undescribed Heliopora sp. group was distributed mainly in the Western Pacific (Japan and Taiwan) and Singapore in this study, while H. hiberniana group was distributed mainly in the Indian Ocean (northwestern Australia, the Andaman Sea in Thailand, and the Maldives in this study; distributed also in Indonesia, Richards et al., 2020). Such uneven geographical distributions imply the origin of Heliopora sp. group could be in the Pacific Ocean and that of H. hiberniana could be in the Indian Ocean, but this remains to be tested with a dated phylogeny. Cenozoic tectonic events (Hall and Holloway, 1998) and sea-level fluctuations associated with climate change (Pillans et al., 1998) in the area between the Indian and Pacific Oceans are thought to have promoted allopatric speciation in many coral reef organisms (Carpenter et al., 2011) including starfish (Benzie, 1999; Vogler et al., 2008) and reef fishes (Gaither et al., 2011; Bowen et al., 2013). It is, therefore, possible that the genus Heliopora has also allopatrically speciated into the two groups, Heliopora sp. and H. hiberniana, due to the physical division between the Indian and Pacific Oceans caused by past tectonic and climatic changes. Morphological and molecular analyses on other Indo-Pacific reef-building coral species also found sibling species between the Indian and Pacific Oceans (Veron, 1995; van Oppen et al., 2001; Pinzón et al., 2013; Huang et al., 2014; Arrigoni et al., 2020; Wepfer et al., 2020).

The H. coerulea group is distributed in both the Indian and Western Pacific Oceans. The distribution of H. coerulea group partially overlaps with the Heliopora sp. group in the western Pacific and with the H. hiberniana group in northwestern Australia. Previous studies revealed sympatrically distributed Heliopora spp. have different reproductive timing: the reproductive timing of Heliopora sp. group (HC-A) and H. coerulea group (HC-B) in Japan (Saito et al., 2015: Taninaka et al., 2018) and the Philippines (Villanueva, 2016) differ by almost one month. Heliopora hiberniana and H. coerulea colonies in northwestern Australia also have been found to have different reproductive timings even when found in sympatry (Richards et al., 2018). These patterns suggest that endogenetic regulation may contribute to differences in reproductive timing. A previous study revealed that there are fixed species-specific substitutions in the dopamine receptor 2-like gene (Isomura et al., 2013) and cryptochrome-1 (Oldach et al., 2017), both of which were identified as the genes that regulate reproductive timing in Acropora species (Iguchi et al., 2019). It is possible that allochronic speciation have occurred between these sympatrically distributed Heliopora spp. groups, or at least allochronism in reproductive timing plays a central role to keep their species boundaries after speciation. The difference in reproductive timing is an important mechanism as a prezygotic isolation for marine species that release gametes into the water column (Knowlton, 1993; Palumbi, 1994). Indeed, many of the sympatrically distributed broadcast spawning coral species have different reproductive timing in the same genus (e.g., Orbicella, Knowlton et al., 1997; Szmant et al., 1997; Acropora, Fukami et al., 2003; Nakajima et al., 2012; Ohki et al., 2015). Generally, the reproductive timing of corals is also strongly associated with environmental factors (e.g., moon phases, moonlight, solar radiation, water temperature) in addition to genetic factors (e.g., circadian clock genes, photoreceptor proteins) that correlate with the environmental factors (e.g., Levy et al., 2007; Brady et al., 2009; Crowder et al., 2017; Oldach et al., 2017). Thus, future studies on reproductive timing of Heliopora spp. that include both environmental and genetic factors may provide further insights into speciation of sympatrically distributed Heliopora spp.

The Heliopora sp. group was found both in Yaku Island in Japan (N 30°16′21.45″, Nakabayashi et al., 2017), the northernmost site and at a depth of 50 m on Green Island in Taiwan, the deepest habitat records of the genus Heliopora. Therefore, contrary to previous knowledge on the ecological traits of H. coerulea that it prefers warm shallow water habitats (Zann and Bolton, 1985), Heliopora sp. group could be found in colder and deeper habitat; the northernmost known habitat of H. coerulea group (HC-B) is Miyako Island (N 24°51′52.10″, Yasuda et al., 2014) and that of Heliopora sp. group (HC-A) is Yaku Island (N 30°16′21.45″, Nakabayashi et al., 2017). In the field, contrasting ecological differences have been previously observed for Heliopora sp. group (HC-A) and H. coerulea group (HC-B). Colonies of Heliopora sp. group (HC-A) are relatively small and more abundant outside the well-developed fringing reef, where water temperature is lower, while colonies of H. coerulea group (HC-B) can be larger, forming micro-atolls, and dominating inner reef habitats where water temperature is higher (Yasuda et al., 2010; Taninaka et al., 2018), providing evidence of niche differentiation which further support the hypothesis of ecological speciation between these two groups.

There are some specific growth forms within each Heliopora group. Columnar and encrusting growth forms are relatively distinct from other morphologies and specific to Heliopora sp. group (Figure 3B, k–o). Slender-branching growth form with whitish skeleton (Figure 3B, f–g) is also specific to H. hiberniana group. However, lobate growth form (Figure 3B, a–d and h–j) can be found in both H. coerulea and H. hiberniana groups. Like other reef-building coral species, young and/or juvenile stages are hard to distinguish among closely related species (Babcock, 1992; Babcock et al., 2003). Previous studies have also reported intermediate growth forms sometimes rendering it difficult to distinguish even between HC-A (Heliopora sp. group) and HC-B (H. coerulea group) possibly due to infrequent hybridization and phenotypic plasticity in the field (Yasuda et al., 2014; Taninaka et al., 2018). Therefore, not all colonies can be clearly classified into the Heliopora groups based on field observation alone.

Phenotypic plasticity and continuous morphological variation among closely related coral species are a common phenomenon and have hindered accurate morphological species identification (e.g., Forsman et al., 2009; Marti-Puig et al., 2014; McFadden et al., 2017). Light intensity and water movement are the most influential environmental factors responsible for coral phenotypic plasticity (Todd, 2008), while other physical environmental factors (e.g., dredging, sediment disturbances, turbidity, competition) are also correlated with morphological diversity, leading to continuous morphological variation of different species under similar environmental conditions (e.g., Darling et al., 2012; Erftemeijer et al., 2012; Swierts and Vermeij, 2016). However, it is still unclear what environmental factors influence gross morphological diversity in Heliopora spp.

Richards et al. (2018) previously reported the key morphological characteristics to distinguish H. coerulea from H. hiberniana, including gross morphology, the intensity of bluish skeletal color, and skeletal microstructure (e.g., highly elaborated coenchymal echinulations and smaller and more numerous autopores), indicating skeletal color can be a key morphological characteristic. In the western Pacific, a few colonies with whitish skeletons have been found in the Heliopora sp. group in Japan (H. Taninaka and N. Yasuda pers. obs.). Contrary to the case in the northwestern Australian region, there was no genetic difference between the whitish and other bluish skeleton colonies, indicating skeletal color cannot be a universal key morphological character to distinguish Heliopora sp. from the other species.

Based on these findings, it is necessary to clarify the relationship between morphological diversity including fine skeletal structure and environmental and/or genetic factors among Heliopora spp. in more detail in the future.

Our phylogenetic tree showed substantial phylogenetic structure and multiple genetically isolated subclades exist within each Heliopora species group. These subclades partially correspond to typical growth forms and geographical distributions (Figures 1, 2 and 3A). In the Heliopora sp. group for example, the geographically distant Singapore subclade (sp.1) first diverged from Japan-Taiwan clades (sp.2 and 3). Then, the Japan-Taiwan clades were split into only columnar growth form subclade (sp.2) found only in Japan and encrusting and columnar growth forms subclade (sp.3) found in both Japan and Taiwan. These two subclades were genetically distinct despite being sympatric. In the field, some ecological differences can be observed between the two subclades: Heliopora sp.2 tends to be found in shallow inner reef areas with strong sunlight and relatively weak waves. Heliopora sp.3 tends to be found in relatively deep darker areas such as outer reef slope with fast water currents and turbidity (H. Taninaka pers. obs.). Further ecological and physiological differences between Heliopora sp.2 and Heliopora sp.3 would clarify the species status of the two subclades. The Heliopora sp.2 and Heliopora sp.3 subclades could not be distinguished in the previous studies using traditional genetic markers (e.g., mitochondrial DNA, ITS2, microsatellite and microsatellite flaking region sequences) (Yasuda et al., 2014, 2015; Taninaka et al., 2019), while this study corroborated the effectiveness of using genome-wide SNPs to elucidate detailed genetic boundaries between closely related lineages/species as in the case of other octocorals (Herrera and Shank, 2016; Quattrini et al., 2019; Takata et al., 2019).

In the H. hiberniana group, two subclades (hib.1-2) showed deep genetic divergence corresponding to geographical distance, namely the northwestern Australian subclade (hib.1) and the Maldives-Thailand subclade (hib.2). This genetic isolation could be the result of long-term restriction of gene flow between the two subclades. The current circulation system in the Indian Ocean mainly consists of currents that move from east to west that does not directly connect northern (Maldives and Thailand) and southern (northwestern Australia) sides of the mid-east Indian Ocean (Schott and McCreary, 2001; Schott et al., 2009). Due to such current patterns, previous studies on other marine organisms with pelagic larval phase also showed similar genetic isolation in the Indian Ocean (Benzie et al., 2002; Bay et al., 2004; Hui et al., 2016). On the other hand, geological research (Schettino and Turco, 2011; Keith et al., 2013; Obura, 2016), fossils, and phylogenetic data (Obura, 2012, 2016; Veron et al., 2015) propose a hypothesis that tectonic changes occurred during the Paleogene and the Neogene and may have promoted genetic isolation and speciation events in the Indian Ocean. Such historical geographic events might have promoted genetic divergence of Heliopora spp. in the Indian Ocean, but this requires further testing within a calibrated phylogenetic framework.

In the H. coerulea group, three subclades (coe.1–3) corresponding to geographically distinct regions were found, namely Guam subclade (coe.1), the northwestern Australia subclade (coe.2), and Japan-Taiwan subclade (coe.3). The most genetically isolated subclade that first split from the others is the Guam subclade (coe.1). It is notable that Japan-Taiwan subclade (coe.3) is genetically more closely related to the northwestern Australia subclade (coe.2) in the Indian Ocean than to the geographically closer Guam subclade (coe.1). It is likely that genetic connectivity has been restricted between Guam and Japan-Taiwan since the Guam subclade (coe.1) diverged from the other two subclades early in the history of H. coerulea. This is unusual because genetic connectivity between Mariana Islands region including Guam and the Kuroshio Current region is often strong in most of the marine invertebrate species with larval dispersal period (e.g., Palumbi et al., 1997; Williams and Benzie, 1998; Lessios et al., 2001; Pinzón et al., 2013; Arrigoni et al., 2020; Wepfer et al., 2020; but see exception in Wörheide et al., 2008). Although the reasons for the genetic break between Guam and Kuroshio Current region in Heliopora spp. are unclear, it could be caused by limited larval dispersal capability of Heliopora spp. (Taninaka et al., 2019) and/or local adaptation of each subclade. On the other hand, northwestern Australian populations of other marine species are genetically more similar to west Pacific populations than other Indian Ocean populations (e.g., Williams and Benzie, 1998; Lessios et al., 2001; Pinzón et al., 2013; Arrigoni et al., 2020; Wepfer et al., 2020) but contrary examples also exist where Indian and Pacific Ocean coral populations are divergent (e.g., Richards et al., 2016). Williams and Benzie (1998) discussed that this might be due to Western Australia being recently colonized by recruits from the Western Pacific lineage after the end of the last ice age. Recent genome-wide phylogeographic analysis of reef-building corals have also found such large discrepancies between geographic and genetic distances (Arrigoni et al., 2020; Wepfer et al., 2020), which need to be examined in closer detail in the future.