Speciation among sympatric lineages in the genus Palythoa (Cnidaria: Anthozoa: Zoantharia) revealed by morphological comparison, phylogenetic analyses and investigation of spawning period

Zoantharians are sessile marine invertebrates and colonial organisms possessing sexual and asexual reproductive ability. The zooxanthellate zoantharian genus Palythoa is widely distributed in coral reef ecosystems. In the Ryukyu Archipelago, Japan, sympatric Palythoa tuberculosa and P. mutuki are the dominant species of this genus in the intertidal zone. Previous phylogenetic analyses have shown that these two species are closely related, and additionally revealed a putative sympatric hybrid species (designated as Palythoa sp. yoron). In this study, we attempted to delineate Palythoa species boundaries and to clarify the relationships among these three groups plus another additional putative sympatric species (P. aff. mutuki) by multiple independent criteria. The morphology of these four lineages was clearly different; for example the number of tentacles was significantly different for each species group in all pairwise comparisons. From observations of gonadal development conducted in 2010 and 2011, P. sp. yoron and P. aff. mutuki appear to be reproductively isolated from P. tuberculosa. In the phylogenetic tree resulting from maximum likelihood analyses of the ITS-rDNA sequence alignment, P. tuberculosa and P. sp. yoron formed a very well supported monophyletic clade (NJ = 100%, ML = 95%, Bayes = 0.99). This study demonstrates that despite clear morphological and/or reproductive differences, P. tuberculosa and P. sp. yoron are phylogenetically entangled and closely related to each other, as are P. mutuki and P. aff. mutuki. Additionally, no single molecular marker was able to divide these four lineages into monophyletic clades by themselves, and a marker that has enough resolution to solve this molecular phylogenetic species complex is required. In summary, the morphological and reproductive results suggest these lineages are four separate species, and that incomplete genetic lineage sorting may prevent the accurate phylogenetic detection of distinct species with the DNA markers utilized in this study, demonstrating the value of morphological and reproductive data when examining closely related lineages.


INTRODUCTION
Zoantharians are sessile marine invertebrates and colonial organisms possessing sexual and asexual reproductive ability (Ryland, 1997). Zoantharians belong to subclass Hexacorallia (Cnidaria, Anthozoa) and they have the significant feature of embedding small particles (sand, detritus) into their body column. Zooxanthellate zoantharian species are found worldwide in tropical and subtropical shallow water areas (Trench, 1974;Reimer, Takishita & Maruyama, 2006).
Traditionally, zoantharian classification has been based on morphological characters such as the relative degree of coenenchyme development, number of tentacles per polyp, oral disk diameter, and position and features of the sphincter muscle (Ryland & Lancaster, 2003). However, sand encrustation (Reimer et al., 2010) and large intraspecific variation have often made histological classification difficult (Muirhead & Ryland, 1985;Mueller & Haywick, 1995;Reimer et al., 2010). Phylogenetic work using mitochondrial 16S ribosomal DNA and cytochrome oxidase subunit I (mtCOI) and the nuclear internal transcribed spacer region of ribosomal DNA (ITS-rDNA) as molecular markers have begun to reveal evolutionary relationships in this group (e.g., Reimer et al., 2004;Sinniger et al., 2005;Reimer et al., 2007b).
The zooxanthellate zoantharian genus Palythoa Lamouroux, 1816 is widely distributed in coral reef ecosystems as a common group of organisms. In the Ryukyu Archipelago of southern Japan (Fig. 1), Palythoa tuberculosa (Esper, 1805) and P. mutuki Haddon & Shackleton, 1891 are the dominant species of this genus in the intertidal zone (Irei, Nozawa & Reimer, 2011). Reimer et al. (2007a) showed that these two species are closely related with phylogenetic analyses based on ITS-rDNA and mtCOI. Furthermore, they revealed a putative hybrid species (designated as Palythoa sp. yoron), which was presumed to have originated via interspecies hybridization between P. tuberculosa and P. mutuki, based on shared additive patterns of nucleotide polymorphisms of ITS-rDNA sequences, and indicated a potential reticulate evolutionary history in these three species groups. A subsequent investigation conducted by Shiroma & Reimer (2010) revealed that P. sp. yoron was sympatric in the intertidal zone with these two other species in Okinawa, but also was present in a different microenvironment than P. tuberculosa and P. mutuki. As well, P. sp. yoron is intermediate in morphological form between P. tuberculosa and P. mutuki. (Fig. 2, Table 1), with all three species readily distinguishable from one another (Shiroma & Reimer, 2010).
In this study, we attempted to determine the delimitation of Palythoa species boundaries and to clarify the relationships among species groups using multiple independent criteria. We first made primary hypotheses of species delimitation based on morphology and habitat preference. We then re-examined these hypotheses via genetic data and investigated ovary

Specimen collection
Specimens of Palythoa species were collected in the intertidal zone from several sites in the Ryukyu Archipelago, including Okinawa-jima Island, Yoron-to Island, Okinoerabu-jima Island, and Tokunoshima Island (Fig. 1, Table 2) between March 2010 to October 2012. All specimens were stored in 99.5% ethanol for DNA analyses or 5% formalin-SW solution for morphological and anatomical analyses. Each specimen was identified according to morphological classification methodology (Pax, 1910), supplemented with a key to field identification (Reimer, 2010), and ecological and morphological aspects of P. sp. yoron (Shiroma & Reimer, 2010). Characters employed for identification of Palythoa species were environment (habitat), coenenchyme development, polyp structure, number of polyps per colony, and numbers of tentacles per polyp. All specimens were identified preliminarily as Palythoa tuberculosa ( Fig. 2A), P. mutuki (Fig. 2B) and P. sp. yoron (Fig. 2C). During collection, it was noticed that certain specimens had a similar external appearance with P. mutuki but with less well developed marginal ridges and larger polyp sizes. Such specimens were found sympatrically with other specimens, and these were designated as P. aff. mutuki (Fig. 2D). In addition, spawning timing investigations for all species groups were carried out between June 2010 to December 2010, and from June 2011 to February 2012 at Kaminomine, Tokunoshima, Kagoshima (27 • 46 09 N, 129 • 02 16 E) by monthly sampling. In particular, for collecting P. tuberculosa, investigation was conducted in a wide area from lagoon tide pools to the outer reef in 2010. However, in 2011-2012 investigations were conducted only in tide pools due to rough sea conditions. At least five different colonies of approximately ten polyps for each species were collected in whole or partially.

External anatomy
Fixed specimens were cut horizontally at the oral disk height by surgical knife and tweezers under stereomicroscope (S8APO, Leica, Tokyo) and the number of tentacles, which is one of the characters for Palythoa species (e.g., Ryland & Lancaster, 2003), were counted (Table 3). To eliminate pseudo-replication in comparison among species, a single polyp was chosen with the table of random number from each colony. The mean numbers of tentacles per polyp for each species pair were compared using Mann-Whitney U test with Bonferroni correction.

Spawning period investigation
Ovary development of all preserved colonies was observed via cross sections made by cutting polyps vertically through the mouth located at the center of oral disk under a Table 3 The mean number of tentacles ± standard deviation and results of Mann-Whitney U test with Bonferroni correction between each Palythoa species pairs. N = total number of examined polyps for each species (one per colony). stereomicroscope. During anthozoans' oogenesis, oocytes form a single-layered germinal ribbon down the mesoglea of the central third of the septa. Subsequently, the germinal ribbon develops a sequence of swollen nodes where the septum folds locally in an S and the layers fuse (Ryland, 1997;Ryland, 2000). When we observed a germinal ribbon in a polyp, we counted the polyp as ''possessing developing ovaries'', and the number of polyps possessing developing ovaries were totaled. To evaluate the spawning period of each species, the ratio of the number of polyps possessing developing and/or developed ovaries to the total number of polyps examined was calculated over time. When the calculated proportion of developed/developing ovaries dropped dramatically, we designated this as the start of the estimated spawning period. The end of the estimated spawning period was defined as the point where the number of developed/developing ovaries reached 0%.

DNA extraction, PCR amplification and direct sequencing
DNA from each specimen was extracted using a DNeasy Blood and Tissue Kit (QIAGEN, Tokyo, Japan) according to the manufacturer's instructions. A small amount of tissue from each specimen was removed using a surgical knife sterilized by open flame. Extracted DNA was subsequently stored at −20 • C, and then we amplified target sequences via polymerase chain reaction (PCR). Three molecular markers that have previously been used for differentiation of Palythoa were chosen; (1) the mitochondrial 16S of ribosomal DNA (mt 16S-rDNA), (2) the mitochondrial cytochrome c oxidase subunit I (mtCOI), and (3) the internal transcribed spacer region of nuclear ribosomal DNA (ITS-rDNA) (Reimer et al., 2004;Sinniger et al., 2005;Reimer et al., 2007a, etc.). Furthermore, a nuclear housekeeping gene, (4) asparaginelinked glycosylation 11 protein (ALG11) region, was also examined for the first time in zoantharians. This marker has been found to be more informative than mtCOI in examining sponge relationships and succeeded in solving previously debated nodes (Hill et al., 2013) and has also been considered to be useful for resolving cnidarian relationships (Belinky et al., 2012). Thermal cycler programs were set to the following conditions: (1) mt 16S-rDNA; an initial denaturing step at 94 • C for 2 min, followed by 40 cycles of 30 s 94 • C, 1 min annealing at 52 • C and 2 min extension at 72 • C, followed by 5 min final elongation at 72 • C with Zoantharia-specific primer set 16Sant1a (5 -GCC ATG AGT ATA GAC GCA CA-3 ) and 16SbmoH (5 -CGA ACA GCC AAC CCT TGG-3 ) (Sinniger et al., 2005); (2) mtCOI; 1 min at 95 • C, then 35 cycles: 1 min at 95 • C, 1 min at 40 • C and 90 s at 72 • C, followed by 7 min at 72 • C with the universal primers HCO2198 (5 -TAA ACT TCA GGG TGA CCA AAA AAT CA-3 ) and LCO1490 (5 -TAA ACT TCA GGG TGA CCA AAA AAT CA-3 ) (Folmer et al., 1994); and (3) ITS-rDNA; 1 min at 95 • C , then 35 cycles of 1 min at 94 • C, 1 min at 50 • C, and 2 min at 72 • C, followed by 10 min at 72 • C with Zoantharia-specific primers Zoan-f (5 -CTT GAT CAT TTA GAG GGA GT-3 ) and Zoan-r (5 -CGG AGA TTT CAA ATT TGA GCT-3 ) (Reimer et al., 2007a).
Amplification for the remaining coding region (ALG11) was performed by touch-down PCR and nested PCR because of low numbers of copies in the whole genome as this is a single-copy gene. For ALG11, although we basically followed the original protocols (Sperling, Pisani & Peterson, 2007;Belinky et al., 2012), some modifications were required to fit the thermal cycler we used, and the conditions were as follows: (4) ALG11 first touchdown, 2 min at 95 • C, then 13 cycles of 1 min at 95 • C, 1 min at 52-40 • C (dropping one degree for each cycle), 1.5 min at 72 • C; followed by 20 cycles of 1 min at 95 • C, 1 min at 52 • C, 1.5 min at 72 • C; lastly 5 min at 72 • C with primers ALG11-D1 (5 -TTY CAY CCN TAY TGY AAY GCN GGN GG-3 ) and ALG11-R1 (5 -ATN CCR AAR TGY TCR TTC CAC AT-3 ), and (5) MAT-f (5 -GGN GAR GGN CAY CCN GAY AA-3 ). In the second touchdown procedure an amplicon of the first touchdown was utilized as the template, followed by 2 min at 95 • C, then 35 cycles of 1 min at 95 • C, 1 min at 52 • C, and 1.5 min at 72 • C. In the end, nested PCR was performed with 2 min at 95 • C, and then 35 cycles of 1 min at 95 • C, 1 min at 52 • C, and 1.5 min at 72 • C with primers ALG11-D2 (5 -TGY AAY GCN GGN GGN GGN GGN GA-3 ) and ALG11-R2 (5 -CCR AAR TGY TCR TTC CAC ATN GTR TG-3 ).
Amplicons were outsourced for sequencing to a private sequencing company (Fasmac Co., Ltd., Kanagawa, Japan) on an Applied Biosystems 3730xl DNA sequencer, using BigDye Terminator V3.1 and the same primer sets as for PCR as described above. Sequence data were edited using BioEdit v.7.2.0 (Hall, 1999).
As numerous indels (inserts and deletions) were confirmed in ITS-rDNA sequences, alignment was performed using ClustalW (Thompson, Higgins & Gibson, 1994) with gap penalties of 10 for open and 1 for extended, followed by manual fixing for obviously misaligned areas such as gap position. Sequences of the 5.8S rDNA region located between internal transcribed spacer 1 (ITS1) and internal transcribed spacer 2 (ITS2) were removed   (2006) from analyses because the substitution rate is apparently lower than ITS1 and ITS2, and an admixture of different substitution rates could lead to a misleading choice of the appropriate substitution model. Additionally, in order to not overestimate for genetic distance in following phylogenetic analyses, sites were removed if they had a percentage of gaps and/or ambiguous sites higher than 95% (partial-deletion option). Fifty-six out of sixty-five specimens had one or more degenerate codes in sequences of the ALG11 region. All degenerate codes were divided into two standard bases using PHASE v2.1.1, which implements a Bayesian statistical method for reconstructing haplotypes from population genotype data (Stephens, Smith & Donnelly, 2001;Stephens & Scheet, 2005). Furthermore, first and second codon positions were removed from the dataset by checking amino acid sequences after translation.

Substitution model selection
Substitution models for each gene were estimated by jModelTest v2.1.3 (Darriba et al., 2012) through the following steps. Initially, likelihood calculations were carried out for all substitution models with configurations of seven substitution schemes, equal or unequal base frequencies (+F), rate variation among sites with a number of rate categories (+G, nCat 5) and base tree topology (ML optimized). Subsequently, the most appropriate model for each marker was selected under (i) the corrected Akaike information criterion (AICc) for Maximum-Likelihood and neighbor-joining phylogenetic estimation, or (ii) Bayesian information criterion (BIC) for Bayes estimation. Thus, the (i)TrN/(ii)TrNef for mt 16S-rDNA, (i)F81/(ii)JC for mtCOI, (i,ii)K80+ for ITS-rDNA, and (i)K80+ /(ii)TPM1uf+ models for ALG11 were employed, respectively.

Gene tree estimations
For four distinct datasets (mt 16S-rDNA, mtCOI, ITS-rDNA, ALG11), phylogenetic analyses were applied independently with the optimal substitution model under AICc estimated by jModelTest. Maximum-Likelihood (ML) analyses were performed using PhyML (Guindon & Gascuel, 2003) and neighbor-joining (NJ) methods were performed using MEGA5.2.2 (Tamura et al., 2011). All other parameters besides substitution model and the discrete gamma distribution were implemented with the default value. Bootstrap analyses (Felsenstein, 1985) of 1,000 replicates were tested to evaluate the support of every branch.
Bayesian inference for gene trees was performed using BEAST v.1.7.0  with the optimal substitution model under BIC. All parameters were used as default values except for the molecular clock, in which the rate was changed to the log-normal relaxed model, while only the substitution model for ALG11 was modified to TPM1uf after generating the initial setting file. Four Markov chain Monte Carlo (MCMC) simulations were run for 10 million generations with sampling intervals of 1,000. Convergence of analyses and adequacy of the sample sizes, with ESS values above 200 (ESS = the number of effectively independent draws from the posterior distribution that the Markov chain is equivalent to) were confirmed in Tracer v.1.5. (Rambaut et al., 2013). Analyses were combined using LogCombiner v.1.8.0, which is included within BEAST, after excluding the first 10% as burn-in. Obtained trees were summarized in a maximum clade credibility tree using TreeAnotator v.1.8.0 and visualized in FigTree v.1.4.0.

Species tree estimations
*BEAST estimates the species tree directly from the sequence data, nucleotide substitution model parameters and the coalescent process (Heled & Drummond, 2010). The species trees were built by grouping all 235 sequences by putative species groups and simultaneously estimating each of three individual gene trees (mt 16S-rDNA, ITS-rDNA and ALG11), and the summary species trees using BEAST were drawn for two different species model; (1) a six species model including P. tuberculosa, P. sp. yoron, P. mutuki, P. aff. mutuki, P. sp. sakurajimensis sensu Reimer et al. (2007a) and Reimer et al. (2007b) and P. heliodiscus, and (2) a four species model combining P. sp. yoron with P. tuberculosa, and P. aff. mutuki with P. mutuki, along with P. sp. sakurajimensis and P. heliodiscus.
All parameters were used as default except for; (1) the molecular clock rate, which was changed to the log-normal relaxed model (Drummond et al., 2006), (2) the substitution rate for mt 16S-rDNA, for which the range was calibrated to between 0.001-0.002/Mya based on the reported substitution rate for mtCOI (Shearer et al., 2002), and (3) the substitution model for ALG11 was modified to TPM1uf after generating the setting file. MCMC analyses were run for 100 million generations with sampling intervals of 10,000 and excluding the first 10% as burn-in. All the parameters in the output file were confirmed in Tracer v1.5. Obtained trees were summarized in a maximum clade credibility tree using TreeAnotator v.1.8.0.

Morphological analyses
The numbers of tentacles were measured for single randomly selected polyps from eleven colonies of P. tuberculosa, eight colonies of P. sp. yoron, seven colonies of P. mutuki, and eight colonies of P. aff. mutuki. The mean number of tentacles ± standard deviation per polyp was 31.6 ± 3.4 for P. tuberculosa, 40.5 ± 2.6 for P. sp. yoron, 54.4 ± 7.4 for P. mutuki, and 71.0 ± 4.1 for P. aff. mutuki. Each respective mean number of tentacles was significantly different (p < 0.01) from all others in all pair tests (Table 3).
For cnidae, many subtle differences in sizes of the various types of cnidae present in different tissues were present (Table 4; Fig. 3). However, the most obvious differences were in small holotrichs, which were rarely observed in the tentacles of column of both P. aff. mutuki and P. mutuki, and additionally observed in the tentacles and pharynx of P. aff. mutuki, but were never observed in tissues of P. sp. yoron or P. tuberculosa (Table 4). However, these small holotrichs were only observed in one out of three specimens each of P. aff. mutuki and P. mutuki, and thus no diagnostic differences were observed in the cnidae of all four species-groups examined (Table 4).
In summary, we could clearly distinguish all four Palythoa species groups based on tentacle numbers (Table 3), as well as gross external morphology (Fig. 2), but not via cnidae analyses (Table 4).

Estimated spawning period
During the initial investigation of June to December in 2010, developed ovaries were observed in P. tuberculosa from the middle of June to the middle of September with decreasing numbers of polyps possessing ova (Fig. 4A, Table 6). Additionally, matured eggs were also observed multiple times (on 28 July and 20 September). In contrast, developed ovaries and matured eggs were observed (Figs. 5A, 5B) only one time (on 26 October) in P. sp. yoron. As well, developing ovaries were observed in P. mutuki from the end of July to the middle of September, however, no matured eggs were observed during this investigation.
In 2011, developed ovaries were observed in P. aff. mutuki on 15 June (Figs. 5E, 5F, 4B), and subsequently developed ovaries were observed in P. sp. yoron in early October and early November (Figs. 5C, 5D), for the second consecutive year. On the other hand, no fully developed ovaries were observed in P. tuberculosa and P. mutuki despite developing ovaries being observed continuously during the summer season (on 23 July, 22 August and 5 October), similar as observed in 2010.

Phylogenetic analyses
Molecular phylogenetic trees mtCOI. The phylogenetic tree resulting from maximum likelihood analyses of the mtCOI sequence alignment is shown in Fig. 6A. Palythoa tuberculosa, P. sp. yoron, P. mutuki and P. aff. mutuki formed one mixed clade with low bootstrap support (Maximum-Likelihood P. mutuki used in previous research (Reimer et al., 2007a;Reimer et al., 2007b; formed one group with sequences from P. sp. sakurajimensis. mt 16S-rDNA. The phylogenetic tree resulting from maximum likelihood analyses of the mt 16S-rDNA sequence alignment is shown in Fig. 6B. Palythoa tuberculosa, P. sp. yoron, P. mutuki and P. aff. mutuki formed one mixed clade with low bootstrap support (ML = 65%, NJ = 64%, B < 0.50). Within this mixed clade, P. mutuki and P. aff. mutuki formed a mixed subclade with low bootstrap support in ML and NJ analyses, however, this monophyletic clade was strongly supported in Bayesian analyses (ML = 64%, NJ = 64%, B = 1.0). Additionally, two sequences of P. mutuki from GenBank that were distinguished from other sequences of P. mutuki in previous research (Reimer, Takishita    ALG11. The phylogenetic tree resulting from maximum likelihood analyses of the ALG11 sequence alignment is shown in Fig. 6D. Compared to the above phylogenetic trees, this tree was the most admixed, regardless of morphospecies. For example, sequences from P. sp. sakurajimensis (used as outgroup here) appeared throughout the tree. Topology comparison between trees. Examining the two outgroups used in this study, Palythoa sp. sakurajimensis was phylogenetically much closer to P. tuberculosa, P. sp. yoron, P. mutuki and P. aff. mutuki compared to P. heliodiscus in every gene tree. There were few differences in sequences from the other four species groups, with only one base pair difference in the mtCOI tree, resulting in P. sp. sakurajimensis' sequences forming one group with some P. mutuki specimens, and only one to two base pairs' difference in  the mt 16S-rDNA tree for all four species groups. In particular, in the ALG11 tree, P. sp. sakurajimensis' sequences were admixed with the other four species groups. Palythoa tuberculosa and P. sp. yoron (designated as the ''Palythoa tuberculosa group'' here), and P. mutuki and P. aff. mutuki (designated as ''Palythoa mutuki group'' here) did not separate into four species groups in each DNA marker's tree. The P. tuberculosa group formed a monophyletic clade in the ITS-rDNA tree and one grouping in the mt 16S-rDNA gene tree with one base difference from the P. mutuki group. On the other hand, the P. mutuki group did not show any common pattern, i.e., admixed with all other species groups except for P. heliodiscus in the ALG11 gene tree, most sequences forming one monophyletic clade with the P. tuberculosa group due to no differences in sequences with some sequences forming one group with P. sp. sakurajimensis due to a one base pair difference from other specimens in the mtCOI tree, forming a monophyletic clade with one subclade in the mt 16S-rDNA tree, and forming a paraphyletic clade with a monophyletic subclade of P. tuberculosa in the ITS-rDNA tree.
Species trees. All hypothetical species were fully supported with posterior probability under both the four and six species models (Figs. 7A, 7B). The divergence time from the most recent common ancestor of P. tuberculosa, P. sp. yoron, P. mutuki and P. aff. mutuki, (divergence of P. sp. sakurajimensis in both cases), was calculated as 147,000 years before present with 95% credible interval [lower 30,900-upper 292,000] under the six species model and as 113,000 years under the four species model with 95% credible interval [lower 25,500-upper 231,000].

DISCUSSION
The purpose of this study was to re-evaluate the systematics of some Palythoa species using an integrative approach. Primary hypotheses of species delimitation were based on external morphology (phenetic criterion) and habitat preferences (ecological criterion). These hypotheses were then examined in the light of additional characters, namely the number of tentacles, spawning periods and genetic data.

Morphology and plasticity
The mean numbers of tentacles were significantly different among specimens of the four putative species; P. tuberculosa, P. sp. yoron, P. mutuki and P. aff. mutuki (Table 3). However, in previous research, the tentacle number of P. tuberculosa has been reported as various ranges, i.e., 30 to 40 (Klunzinger, 1877), up to 50 (Walsh & Bowers, 1971), 38 to 52 (Reimer & Todd, 2009), 30 to 37 (Shiroma & Reimer, 2010, or 30 to 50 (Hibino et al., 2013). A wider range of variations has been reported in P. mutuki, with 88 to 144 (Ryland & Lancaster, 2003), 60 to 74, approximately 80 for P. mutuki-related (Reimer & Todd, 2009), or 42 to 66 (Shiroma & Reimer, 2010 reported. Thus, the ranges of tentacle numbers can be assumed to be 30 to 52 for P. tuberculosa and 42 to 144 for P. mutuki, and therefore tentacle numbers of P. sp. yoron and P. aff. mutuki observed in this study are within ranges of previously reported intraspecific variation. These differences between tentacle numbers reported in the literature and our data may be partly explained by the fact that previous authors did not consider P. sp. yoron and P. aff. mutuki as different species. However, Ong, Reimer & Todd (2013) also demonstrated phenotypic plasticity in P. tuberculosa with high ability to acclimate against changes in light-induced environments. From in situ observations, P. sp. yoron seems to prefer locations exposed to strong current such as extensive reef flats where the back reef moat is widely developed. Correspondingly, P. sp. yoron is also often found in back reef moats, as Shiroma & Reimer (2010) mentioned, covered with sand or other loose detritus. High numbers of tentacles enable them to acquire nutritious detritus and feed on planktonic organisms, but strong-current environments repeatedly cover colonies with sand. From the viewpoint of its small, tetrapod colony shape,  P. sp. yoron seems have adapted to such an environment. Therefore, to ensure whether differences in tentacle numbers and colony form between P. tuberculosa and P. sp. yoron are caused by species differentiation, the observation of reaction norms of each species with transplantation experiments is needed.
Although previous cnidae research with detailed statistical analyses revealed finer-scale differences among Palythoa species (Ryland & Lancaster, 2004), we did not observe any useful diagnostic differences with utility for rapid identification of species groups in this study.

Spawning periods and reproductive isolation
Over the two years analyzed, P. sp. yoron consistently developed ovaries later than the three other putative species. If we assume a sharp drop in the proportion of developed ovaries as the consequence of the release of eggs, the annual spawning period estimated for P. sp. yoron was early to mid-November and that of P. aff. mutuki mid-to late June. The spawning period of P. tuberculosa in Okinawa-jima I. has been reported in early August (Yamazato, Yoshimoto & Yoshihara, 1973), from the end of July to middle August (Shiroma & Reimer, 2010), and on 19 and 20 August in 2009 (Hirose et al., 2011). In our study, spawning was estimated to have occurred in August in 2010 and possibly from early July in 2011. The reproductive season of P. mutuki was presumed that be synchronized with P. tuberculosa in 2010, although developed eggs were not confirmed. Little is known about the sexual reproductive ability of this species, and according to Ryland & Lancaster (2003) the only previous records of P. mutuki possessing developed oocytes are from Fiji and Tuvalu. To overcome this lack of knowledge, closer examinations via staging of histological sections for gonadal development (such as done by Polak et al., 2011) are required.
Interpreting these results in terms of putative reproductive isolation is not straightforward. Even assuming that a sharp drop in the proportion of developed ovaries translates into a major spawning event, which seems to be a reasonable hypothesis, this does not exclude the possibility of eggs being released much later than the initial peak. For example, while we estimated the spawning period of P. tuberculosa to have occurred in August in 2010, nearly 20% of individuals still had developing or developed ovaries on September 20th, which may have been released as mature eggs at any time from then until October 26th (Fig. 4A), and enabled potential cross-fertilization with P. sp. yoron. On the other hand, data thus far indicate spawning on one or two nights per year for brachycneminic zoantharians (Ryland, 1997), and reabsorption of oocytes (Ono et al. 2005) that did not spawn. More work is needed to determine exact spawning patterns of Palythoa tuberculosa and closely Palythoa related species, but the asynchrony of both ovary development (P. sp. yoron) and spawning peaks for P. tuberculosa and P. aff. mutuki suggest that at least partial pre-zygotic reproductive isolation is possible among P. sp. yoron, P. tuberculosa and P. aff. mutuki at Tokunoshima I.

Species boundaries in phylogenetic trees
The four genetic markers analyzed in this study displayed contrasting patterns. The two mitochondrial genes were relatively conservative, as has been reported for other anthozoans (Shearer et al., 2002;Huang et al., 2008), but mt 16S-rDNA allowed the recovery of P. heliodiscus, P. sp. sakurajimensis, the P. mutuki group and the P. tuberculosa group as four genetically homogeneous groups (phenetic criterion), and all species or species groups were reciprocally monophyletic with the exception of P. tuberculosa. ITS-rDNA showed a similar pattern with the P. mutuki group and the P. tuberculosa group represented in distinct clades, although the P. mutuki group was paraphyletic. This consistency across mitochondrial and nuclear markers also suggests that there is no genetic exchange (biologic criterion) between these four groups, and thus provides a first level of species delimitation. In contrast, all Palythoa spp. besides P. heliodiscus were largely mixed in the tree recovered from the ALG11 marker, which strongly suggests incomplete lineage sorting for this gene.
Despite obvious differences in morphology and reproductive season between P. tuberculosa and P. sp. yoron, as well between P. mutuki and P. aff. mutuki, no molecular marker was successful in dividing these species pairs into their own monophyletic clades. Palythoa sp. yoron formed a subclade from two specimens in Reimer et al. (2007a), however, in this study reconstructing phylogenetic trees based on the same genomic region with more specimens of P. sp. yoron, one mixed monophyletic clade was supported well with all the other P. tuberculosa specimens. The same pattern was observed with P. mutuki and P. aff. mutuki. These results imply either gene flow between each pair of nominal species or incomplete lineage sorting. Although these two alternative hypotheses are not mutually exclusive, the absence of intermediate morphotypes and the presence of distinct spawning periods lead us to favor the latter over extensive gene flow. Sequences from other single-copy nuclear markers like ALG11 are required to more thoroughly resolve these two species pairs.

Sympatric speciation timing
Recently, sympatric speciation has come to be understood as a major generator of marine biodiversity (reviewed in Bowen et al., 2013). Under such situations, ecological (e.g., behavior or microhabitat) boundaries lead to isolation. However, the hierarchy of timing of sympatric speciation processes (e.g., the order that separation occurs via phylogenetic, reproductive, and morphological criteria) as lineages diverge remains not well understood, with no clear consensus (Norris & Hull, 2012;Pabijan et al., 2017). For example, in tropical bivalves, phylogenetic differences (=cryptic species) have been observed without any clear evidence of morphological differences (e.g., Lemer et al., 2014). On the other hand, in many marine taxa, it has been proposed that during sympatric speciation, reproductive isolation is one driving force behind lineage divergence (Palumbi, 1994).
In this study, morphology and reproductive data sets showed four Palythoa lineages, while DNA markers showed either two lineages (ITS-rDNA, mtCOI, mt 16S-rDNA) or one admixed lineage (ALG11). Combined molecular analyses suggested either two or four lineages were equally possible (Fig. 7). Such varied results along a speciation continuum between different datasets reflect the patterns to be expected during ongoing or incomplete speciation events (Nosil, Harmon & Seehausen, 2009). As all four Palythoa lineages can be found in sympatry at Tokushima I., our results suggest that reproductive isolation, perhaps caused by past hybridization and back-crossing events (Reimer et al., 2007a;MacLeod et al., 2015), led to the generation of these different lineages and morphological differentiation. Phylogenetic differentiation currently remains incomplete due to the evolutionary recentness of these events, estimated as less than 200,000 years before present. Such confounding data, with reproductive isolation but incomplete genetic lineage sorting, can be expected due to the extended duration of speciation events (Norris & Hull, 2012).

CONCLUSIONS
Overall, the data imply that Palythoa species have a much more complex evolutionary history at the species level than previously expected (e.g., in Reimer et al., 2007a). However, natural hybridization between P. tuberculosa, P. sp. yoron and P. aff. mutuki seems to not be currently occurring, at least for populations at Tokunoshima I. observed in this study. In spite of ambiguous phylogenetic differentiation between P. tuberculosa and P. sp. yoron, and between P. mutuki and P. aff. mutuki, we consider these four lineages are all distinct species based on their morphological differentiation and distinct spawning periods. In situ observation of spawning events combined with genomic level examinations will help further clarify the hierarchy of timing in speciation events, and these four sympatric Palythoa lineages present a potential model system for such studies.