Phylogeography of Blue Corals (Genus Heliopora) Across the Indo-West Pacific

Introduction

Coral reef ecosystems have the highest biodiversity of all marine ecosystems (Fisher et al., 2015) and are important in terms of their scientific, economic, and social value (Moberg and Folke, 1999; Spalding et al., 2001; Cesar et al., 2003). However, reef-building corals, the foundation organisms in coral reefs, have been rapidly deteriorating (Hughes et al., 2003; Altieri et al., 2017; Eyre et al., 2018) due to coral mass bleaching and other threats associated with the on-going climate change (e.g., global warming and ocean acidification) (Hoegh-Guldberg et al., 2007; Hoegh-Guldberg and Bruno, 2010; Hoegh-Guldberg, 2011; Hughes et al., 2018); local anthropogenic impacts (e.g., red-soil runoff, nutrient enrichment, and marine plastic pollution) (Halpern et al., 2008; Omori, 2011; Huang et al., 2020; de Oliveira Soares et al., 2020); outbreaks of corallivorous animals (Birkeland and Lucas, 1990; Yasuda, 2018; Montalbetti et al., 2019) and coral disease (Bruno et al., 2003, 2007; Burge et al., 2014). Approximately one-third of the world’s reef-building corals are facing an elevated risk of extinction (Carpenter et al., 2008; Wilkinson, 2008), but many more cryptic species may be undescribed (Bickford et al., 2007) and if the natural range of species is misunderstood, their threatened status may be under-estimated (Richards et al., 2016).

The classification of reef-building coral species has traditionally been based mainly on morphological characteristics (Veron, 2000; Fabricius and Alderslade, 2001). In some cases, however, morphological classification has been difficult due to a high level of phenotypic plasticity across geographic and depth gradients, and in some genera, there is a lack of discrete morphological features that can be used to underpin species identifications (e.g., Forsman et al., 2009; Marti-Puig et al., 2014; McFadden et al., 2017). More recently, genetic information based on several molecular markers has shown morphology can conceal cryptic evolutionary relationships (Richards et al., 2013) and it is increasingly common for genetically distinct cryptic species or lineages to be found within widespread coral species (e.g., Nakajima et al., 2012; Pinzón et al., 2013; Warner et al., 2015). Such discordance between morphology and genetics has led to confusion regarding the species boundaries in many groups (e.g., Acropora, van Oppen et al., 2001; Marquez et al., 2002, Pocillopora, Flot et al., 2008; Souter, 2010; Schmidt-Roach et al., 2013). Such ambiguity threatens to undermine effective biodiversity conservation efforts and prevents the true complexity of coral reef ecosystems from being properly understood.

Molecular data is not, however, always able to provide a definitive answer about species boundaries within reef-building corals. Molecular systematic studies have been partly hindered by the lack of appropriate genetic markers and the traditionally used mitochondrial markers have suffered from inadequate analytical resolution (Shearer et al., 2002; Ridgway and Gates, 2006; Bilewitch and Degnan, 2011). The ability to interpret nuclear datasets is also hampered by the multi-copy nature of genes (Vollmer and Palumbi, 2004) even those that were originally thought to be single copy (e.g., Pax-C, Rosser et al., 2017), introgressive hybridization events (Veron, 1995; Richards et al., 2013) and/or incomplete lineage sorting (e.g., van Oppen et al., 2001; Souter, 2010; Yasuda et al., 2015). Recent advances in High-Throughput Sequencing (HTS) have enabled the use of genome-wide polymorphisms such as restriction site associated DNA sequencing (RAD-seq, Miller et al., 2007; Baird et al., 2008; Rowe et al., 2011; and Dart-seq, Rosser et al., 2017) and Multiplexed inter simple sequence repeat (ISSR) genotyping by sequencing (MIG-seq, Suyama and Matsuki, 2015) to infer phylogenetic relationships. Such techniques often successfully delimit species of non-model organisms including corals with higher resolution than traditional genetic markers (RAD-seq, e.g., Pante et al., 2015; Herrera and Shank, 2016; Quattrini et al., 2019 and MIG-seq, e.g., Tamaki et al., 2017; Park et al., 2019; Hirai, 2019).

MIG-seq is an easy, cost-effective, novel method to obtain a moderate number of single nucleotide polymorphisms (SNPs) of non-model organisms with polymerase chain reaction (PCR) and HTS technology. The number of SNPs from MIG-seq analysis is generally less than those obtained using other techniques such as RAD-seq. However, MIG-seq has several advantages, since the method can be performed with small amounts and/or low-quality DNA and is relatively easy and cheap. Indeed, MIG-seq successfully revealed species boundaries of octocoral species that could be undetectable by traditional genetic markers (Richards et al., 2018; Takata et al., 2019). While some studies examined genetic relationships of closely related reef-building corals using HTS techniques in geographically restricted regions (e.g., Forsman et al., 2017; Johnston et al., 2017), only a few studies have analyzed speciation and/or genetic divergence patterns covering the Indo-Pacific scale (but see Warner et al., 2015; Nakajima et al., 2017; Arrigoni et al., 2020; Wepfer et al., 2020). For a more comprehensive picture of the phylogeographic structure of coral species in the Indo-Pacific, and to understand the species diversity of corals, it is necessary to obtain more extensive representation of species across their geographic range.

Increasingly recognized is the importance of integrating information about reproductive timing and physiology with genetic and morphological data to obtain more robust estimates of coral species (see Ohki et al., 2015; Rosser, 2015; Villanueva, 2016; Luzon et al., 2017; Richards et al., 2018). Timing of reproduction in broadcast spawning species is particularly important because gametes are viable for only a few hours, so individuals that spawn more than a few hours apart are unlikely to cross-fertilize (Levitan et al., 2011) and in this regard, reproductive timing heavily influences evolutionary patterns. Whilst integrated datasets that include reproductive data and cover a wide geographic area are emerging as the gold standard in phylogenetic studies, practical and logistical constraints often mean it is not possible to obtain complete sampling coverage in every region (Keyse et al., 2014; Wepfer et al., 2021). Nevertheless, examining the species boundaries, ecological preferences, and distribution of reef-building corals is fundamental for promoting effective conservation and management strategies.

The blue corals (genus Heliopora) are members of octocorals, though it has a well-developed aragonite skeleton like those of scleractinians (Quattrini et al., 2020). Owning to their hard skeleton, Heliopora spp. play an important role in coral reef accretion in various Indo-West Pacific reefs (e.g., Zann and Bolton, 1985; Planck et al., 1988; Abe et al., 2008; Takino et al., 2010), and therefore is an important reef-building coral taxon. Heliopora coerulea (Pallas, 1766) has mainly three growth forms, namely lobate, columnar, and encrusting (Dana, 1846), but with high morphological plasticity and continuous variations among growth forms (Eguchi, 1948; Colgan, 1984; H. Taninaka pers. obs.). Heliopora coerulea has long been regarded as a living fossil, being the sole extant member of the family Helioporidae (Colgan, 1984) until the recent description of Heliopora hiberniana Richards et al., 2018, a morphologically, genetically and reproductively distinct species from northwestern Australia in the Indian Ocean (Richards et al., 2018). To date, H. coerulea is known to be distributed throughout the Indo-Pacific realm (Wells, 1954; Bouillon and Houvenaghel-Crévecoeur, 1970; Zann and Bolton, 1985), whereas H. hiberniana is known from the north-west shelf of Western Australia, the Maldives, and the Wakatobi and Gili Islands in Indonesia (Richards et al., 2018, 2020).

Heliopora spp. are gonochoric, brooding corals (Babcock, 1990) with short larval dispersal durations (Harii et al., 2002; Harii and Kayanne, 2003). Due to such a low larval dispersal ability, genetic connectivity among geographically separate populations may be limited, a pattern that was evident within the scale of several hundred square kilometers (Taninaka et al., 2019). Therefore, it is possible that Heliopora includes further genetically distinct lineages in the Indo-Pacific. Previous genetic studies have indicated that additional cryptic species may indeed be present within H. coerulea along the Kuroshio Current region of the Philippines, Taiwan and Japan (Yasuda et al., 2014). Study from Yasuda et al., 2014 hypothesized that there are at least two genetically distinct lineages (HC-A and HC-B) of H. coerulea in the Kuroshio Current region that can also be distinguished by different gross morphologies (small branch for HC-A and flat shapes for HC-B, Yasuda et al., 2014, 2015; Iguchi et al., 2019). Subsequent observations (Saito et al., 2015; Villanueva, 2016) and histological examination (Taninaka et al., 2018) revealed that the reproductive timing of the two lineages is different by almost one month even under similar environmental conditions. However, the phylogeographical relationships among HC-A and B in the Kuroshio Current regions, and with H. hiberniana and H. coerulea in northwestern Australia are still unknown. The aim of this study was to clarify the phylogenetic relationships and distributions of the genus Heliopora across the Indo-West Pacific to gain insights on the speciation patterns of Indo-Pacific coral species. To achieve this goal, we collected Heliopora samples widely from the Indo-West Pacific region and applied genome-wide phylogenetic, phylogeographic, and population genetic analyses using the MIG-seq method.

Materials and Methods

Sample Collection and DNA Extraction

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.

MIG-seq Library Preparation and Sequencing

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).

Sequence Processing and SNP Genotyping

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).

Assign Clones

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.

Phylogenetic Analysis

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).

Population Genetic Analysis

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.

Results

Number of Reads and SNPs From MIG-seq Analysis

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.

Phylogenetic Relationships

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).

Individual-Based Genetic Differentiation

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).

Geographical Distribution

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).

Gross Morphology

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.

Discussion

In the present study, we examined the phylogenetic relationships, gross morphological variations, and the geographical distributions of Heliopora spp. collected from the Indo-West Pacific region. Our phylogenetic reconstruction indicates that there are three genetically isolated species groups in the genus Heliopora and we provide evidence that there is substantial subcladal genetic structuring within each group.

We do not describe the newly found Heliopora sp. group or its subclades as new species of the genus Heliopora in this study. Cladistic analysis and formal species description are being undertaken separately.

Phylogeography and Speciation Patterns of Heliopora

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.

Gross Morphological Diversity of Heliopora

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.

Subclades Within Each Heliopora Group

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.

Conclusion

The present study provides a comprehensive picture of phylogenetic relationships, distribution, and growth form patterns among closely related Heliopora species that gives new insights into speciation patterns of Heliopora spp. in the Indo-West Pacific. Our genome-wide genetic data of Heliopora spp. revealed three major groups with eight subclades. Heliopora sp. group was genetically distinct from the other two groups (H. coerulea group and H. hiberniana group) containing the samples from type localities of the previously described species. Genetic diversification in the genus Heliopora could be attributed to allopatry, allochrony, and local and/or ecological adaptation. Further studies incorporating more genome-wide markers for phylogeographic analysis, geohistorical information, physiological and ecological data would more precisely delineate the boundaries of Heliopora species.

Our study also demonstrates the effectiveness of the MIG-seq method for clarifying the species boundaries of octocoral species that were indistinguishable using traditional genetic markers such as mtDNA and ITS2. In addition, we highlight the importance of wide geographic sampling to generate a more complete picture of the phylogeographic relationships and speciation patterns among closely related octocoral species and how they diversified across the Indo-Pacific. It is expected that the MIG-seq method can be applied to other coral species and non-model organisms to resolve longstanding questions in marine evolution.

Field Study Permission

Domestic sampling except for Amami Oshima was conducted under each prefectural government permits; Yaku Island (Kagoshima Prefectural Government permit, No. 2–56), Okinawa including Okinawa Mainland, Kume Island, Ishigaki Island, Sekisei Lagoon, and Iriomote Island (Okinawa Prefectural Government permits, Nos. 18–34, 19–39, 19–60, 21–18, 22–15, 25–44, 26–10, 27–29, 27–81, 28–78, 30–42, and 31–18), and Amami Oshima was conducted under personal permit by Mr. Katsuki Oki. Foreign sampling was conducted under each local government permits; Taiwan (samples collected from Dabaisha and Gongguan, the Taitung County Government, No. 1074150258; a 50 m specimen collected in Chaiko, Taitung County Government, No. 1040000285, and that is registered at the Zoological collection of the Biodiversity Research Museum in Academia Sinica AZIC0001313), Guam (SC-19-006 (occ), University of Guam Marine Laboratory), Singapore (National Parks Board Permit NP/RP16-156), Australia (samples in Western Australia were collected under permit number AU-COM2013-203 and Australian Government approval number 006-RRRW-130723-01), Thailand (the permission for collecting and sending the specimens is granted by the Department of Fisheries), Maldives (Permit No. (OTHR)30-D/INDIV/2017/123 released by the Maldivian Ministry of Fisheries and Agriculture). All the 1st PCR except for those from Thailand were conducted on site and only 1st PCR products were brought or sent to Japan. Ethanol-preserved coral fragment samples from Thailand were imported to Japan in 2011 (CITES number: AC.0510.2/407), and we used them in this study under an additional permit from the Department of Marine and Coastal Resources in Thailand.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: DDBJ (accession: DRA012077), FigShare (doi: 10.6084/m9.figshare.14356418 and doi: 10.6084/m9.figshare.14814231).

Author Contributions

HT and NY conceived the study and drafted the manuscript. All authors except for TK and HY collected samples and conducted molecular genetic experiment. HT, TK, and HY conducted genetic analysis. All authors edited the draft and approved it for publication.

Funding

Funding of this study was provided by the Environmental Research and Technology Development Fund (4RF–1501 and 4–1304) of the Ministry of the Environment, Japan; a Grant-in-Aid for Young Scientists (A) (17H04996); a Grant-in-Aid for Research Fellows (19J21342) by The Japan Society for the Promotion of Science (JSPS); and the Marine Science R&D Programme (MSRDP-P03) supported by the National Research Foundation, Prime Minister’s Office, Singapore. ZR was supported by a JSPS short-term visiting Fellowship to the University of Miyazaki and samples from Western Australia were collected with the support of Woodside Energy on the Woodside Kimberley Collection Project.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to Eriko Nagahiro for her assistance for genetic analysis. DM would like to thank Marco Casartelli and Julian Sitemba for their help in laboratory activities.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.714662/full#supplementary-material

Supplementary Figure 1 | A bar plot of Heliopora sp. estimated by STRUCTURE ver. 2.3.4 assuming the number of cluster K = 3 using 225 SNPs obtained from the Stacks (r = 0.9 with the option –write_single_snp). Ten independent runs based on 200,000 burn-in followed by 200,000 Markov chain Monte Carlo (MCMC) replications were conducted using admixture models assuming correlated allele frequencies among 84 samples with uniform prior. The x-axis indicates each individual for Heliopora sp. that were classified into three subclades (sp.1–sp.3) in the phylogenetic tree (Figure 1). Different colors of the y-axis indicate the probability of assignment to different clusters. Red allows indicate possible hybrid individuals.

Supplementary Table 1 | Summary of sample information in this study: name of sampling country “Country,” name of sampling region “Region,” name of sampling site “Location,” longitude and latitude “Coordinates,” number of collected samples “Ns” and analyzed samples “Na,” ratio of assigned possible clones within each population “Rc,” name of Heliopora species group “Spp. group” and subclade “Subclade,” observed growth forms within each subclade “Major form (minor form),” and sampling years with published references “Year (References).”

References