We assessed the power of 3D quantitative morphology to discriminate morphogroups using skeleton specimens of three closely related tabular Acropora species previously delineated using different lines of evidence (i.e., morphology, breeding trials, and molecular analyses): A. cf. bifurcata (n = 28), A. cf. cytherea (n = 21) and A. aff. hyacinthus (n = 25), hereafter species A, B and C respectively (for further information and comparison to type material see Table 1 at Ramírez-Portilla et al., 2022). Briefly, morphospecies were identified in the field following Veron (2000), particularly using the branch taper (either gradually narrowing or cylindrical) and the radial corallites shape (all labellate either with round, straight or flaring lips, see zoom in branches in Figure 1 in this paper and Figure 33 in Wallace, 1999). Subsequently, multivariate morphological analyses of qualitative and quantitative variables, cross-fertilization experiments, and molecular analyses using target capture and Sanger sequencing were used to identify species boundaries in the data set. For the 3D morphology assessment in this study, we documented a total of 74 skeleton fragments deposited as vouchers at the Sesoko Station, Tropical Biosphere Research Center (TBRC); collected in 2015, 2018, and 2019 from the outer reef south of Sesoko Island (26.6288 North, 127.8622 East, Okinawa, Japan). Documented specimens corresponded to medium-size fragments (min. area 8x8 cm) collected from adult colonies with similar sizes ( Supplementary Table S1 , photos available in Morphobank Project 4065, http://morphobank.org/permalink/?P4065).

The 3D scanning of the coral fragments was performed using a handheld Artec 3D Space Spider Scanner coupled with the software Artec Studio v10 (Artec 3D, Luxembourg). For a preliminary assessment of the 3D model quality, scanning was completed using the real-time fusion mode for all fragments ( Supplementary Figure S1 ). Following Reichert et al. (2016; 2017), the Artec Studio software was used to render and clean up the 3D models for which fusion was performed with 0.2 mm resolution ( Supplementary Materials and Methods ). Meshes were then exported as triangulated mesh files (either.stl or.obj; available in Morphobank http://morphobank.org/permalink/?P4205) for downstream analyses derived either from triangulated polygon meshes, medial axis skeleton graphs, or a combination of both ( Figure 1 and Supplementary Figure S2 ).

The resulting triangulated polygon meshes were further analysed and visualized using the Visualization Toolkit v9.1.0 (VTK; Schroeder et al., 2006) in Python v3.8 (Van Rossum and Drake, 2009) and the Insight Toolkit v5.1.2 (ITK; Ibáñez et al., 2003) in c++20 (ISO/IEC, 2020). The surface area (SA) and the volume (V) of the specimens were obtained with vtkMassProperties from which the surface-area-to-volume (S/V) ratio and the sphericity (φ) were estimated directly ( Supplementary Materials and Methods ).

Four common characteristics of surface curvature were estimated for each vertex of the polygon mesh following Meyer et al. (2003). First, a discrete approximation of the Gauss-Bonnet theorem and the Laplace-Beltrami operator were implemented to obtain the Gaussian (K) and mean (H) curvature for all vertices respectively ( Supplementary Materials and Methods ). The sign of H was then determined based on the direction of the normal vectors, which were obtained using vtkTriangleMeshPointNormals. Finally, the two principal curvatures, maximum curvature (k1) and minimum curvature (k2), were estimated considering that the Gaussian curvature (K) is defined as the product of the two principal curvatures at that location, and the mean curvature (H) corresponds to average of the two principal curvatures ( Supplementary Figure S3 ).

To capture the topological branching structure of corals and facilitate the estimation of measures related to this type of morphology, we extracted the medial axis skeleton from the previously rendered 3D models using a voxel thinning algorithm. For this purpose, the polygon mesh was first smoothed using the vtkWindowedSincPolyDataFilter module (iterations: 100, pass-band frequency = 0.005), thereby reducing the details of the surface and potential noise while still maintaining the general shape of the coral specimens. Next, the smoothed mesh was transformed into a binary voxel image (resolution = 0.5mm×0.5mm×0.5mm) using vtkPolyDataToImageStencil (tolerance = 0). Finally, voxel thinning was performed with itkBinaryImageThinningFilter3D (Homann, 2007), an implementation of the algorithm of Lee et al. (1994) that results in single-voxel thin skeletons. The voxel skeletons were then transformed into graphs (G) by translating each voxel to a vertex (v) with coordinates corresponding to the represented location and connected by edges using the method described by Reinders et al. (2000). In this graph, a branch (b) was considered to be the set of neighbouring vertices and edges between two successive junction vertices (with a vertex degree higher than two), or between a junction and a terminal vertex (with a degree of one, Supplementary Materials and Methods ).

The branches were then identified from the graph and three different morphological characters were estimated ( Supplementary Materials and Methods ). Branch length (br_length ) was calculated as the sum of all the edge lengths. Branching rate (br_rate ), or how often a coral branches, was defined as the distance between the first (v₀) and last vertex (v_N ) of the branch. Then, two definitions were used to estimate branch spacing. First, following Kruszyński et al. (2007), branch spacing (br_spacing_v1) was defined as the shortest distance between the tip (v_T ) of the terminal branches and any vertex in the skeleton graph not belonging to the current branch. Finally, following Wallace et al. (1991), a second proxy of branch spacing (br_spacing_v2) was defined as the shortest distance between the tip (v_T ) of a terminal branch and any other v_T .

To obtain information about the branch width (br_width ), once the medial skeleton axis and the smoothed polygon mesh were obtained for each specimen, each vertex of the skeleton graph was associated with a medial thickness parameter (d(v)), which represents the diameter of the branch at v. There, the medial thickness was estimated as twice the distance to the closest point on the smoothed polygon mesh. Following Kruszyński et al. (2007), the medial thickness of the individual vertices was translated to three metrics related to the branch width (i.e., the diameter of a sphere at a certain point of the branch): the width at the base of a branch or the junction vertices (da, a−sphere), the width adjacent to the junction a, towards the midsection of the branch (db, b−sphere), and the terminal width (dc, c−sphere) or medial thickness at the tip of the branches (v_T ) ( Figure 2 , upper left). As an additional parameter, the average thickness of each branch (d_avg ) was obtained by averaging the medial thickness of all vertices in the branch. The location of the a−sphere and b−sphere were also used to obtain the angle of the branches (b_angle ). The angles were obtained for all terminal branches (b_T ) and were defined as the smallest angle at its associated a−sphere between its b−sphere and the b−sphere of a neighbouring branch.

Curvature features were also estimated at the branch tips where it was defined from the subset of the vertices on the polygon surface that were located on the tips of the branches. To identify the branch tip vertices the skeleton graph was used ( Supplementary Figure S2 ). For each b_T , a cylinder that had a diameter of da and an axis that followed the direction of the vector between v_T and v_N/2 was placed on v_T together with a plane orthogonal to the cylinder axis. The vertices and faces that were located within the cylinder and exceeded the plane were selected using with the vtkExtractPolyDataGeometry module. From this selection, the set of connected vertices closest to v_T (obtained with vtkPolyDataConnectivityFilter) were considered to be the branch tip vertices (v_tip ) of the polygon mesh. For these vertices, curvature values were calculated for the specimens as previously described for the polygon meshes (see section 2.3).

Three variable types were estimated from the 3D data set; complexity, curvature, and branching ( Figure 2 and Table 1 ). For complexity variables, global values for each one of the coral fragments were obtained (i.e., a single value per specimen). For curvature and branching variables, estimation methods yielded results per branch in the skeleton or per vertex of the polygon mesh. Therefore, to transform these distributions into univariate measures and obtain the average values, certain features were assessed. For the branch-related measures, outliers (|Z| > 3) of each coral were removed and the mean (_mean) and variance (_var) of the distributions were obtained. For the general curvature measures, values within the 2.5–97.5^th percentiles were analysed to obtain the weighted mean (_mean), variance (_var), skewness (_skew) and kurtosis (_kurt) of each distribution. For curvature measures at the branch tips, first the distribution of the number of v_tip per b_T was analysed to remove branches with too many v_tip (Z > 3), as it indicates that reliable estimation of v_tip failed. The remaining branch tip vertices were assembled and analysed in a similar fashion as the general curvature distributions.

To provide a quantitative comparison of the estimated morphological variables, both univariate values and distributions were analysed using R v4.1.0 (R Core Team, 2018) through the RStudio console v1.4.1103 (RStudio Team, 2017). The three species previously delineated in this data set (Ramírez-Portilla et al., 2022) were used as a three-level factor for the subsequent analyses.

To inspect clustering using the complete set of variables and the preliminary selected subset, the most likely number of groups was estimated according to 30 different indices (NbClust v3.0; Charrad et al., 2014), followed by a hierarchical clustering analysis (HCA; cluster v2.1.2; Maechler et al., 2021) using Euclidean distance and three different clustering methods (i.e., Ward, complete, and average) in which p-values were calculated via multiscale bootstrap resampling (pvclust v2.2-0; Suzuki et al., 2019).

A principal component analysis (PCA) was also performed to evaluate the ordination of the subset (stats v4.1.0; R Core Team, 2018). For this purpose, unbiased feature selection was performed using Gaussian model-based clustering (clustvarsel v2.3.4; Scrucca and Raftery, 2018) according to the Bayesian information criterion (BIC). Briefly, a set of variables that best discriminated groups using normal mixture models (NMMs) without a priori information was defined using the greedy algorithm both in forward and backward directions (Raftery and Dean, 2006; Scrucca, 2010). This set of variables was then used to reduce the dimensionality of the data using a PCA ( Supplementary Figure S4 ).

Congruence between morphogroups discriminated using these multivariate approaches were contrasted to the three species previously delineated in this data set by mapping each coral specimen to its corresponding taxonomic assignment in each of the analyses (Ramírez-Portilla et al., 2022).

The discriminative potential of the 3D-based variables was gauged by removing highly correlated features from the complete variable subset according to their variance inflation factor (VIF < 10; usdm v1.1-18; Naimi et al., 2014) to perform a multivariate analysis of variance (MANOVA; stats v4.1.0; R Core Team, 2018), and a linear discriminant analysis (LDA) with the maximum likelihood (ML) estimator method (MASS v7.3-54; Venables and Ripley, 2002). The accuracy of the discriminant approach was assessed by randomly partitioning the data set in a training (n = 50, 67.6% of specimens) and testing (n = 24, 32.4% of specimens) subsets and calculating the corresponding prediction accuracy tables or confusion matrices.

Overall, 53 univariate variables from one of three types (complexity, curvature, and branching, Figure 2 ) were estimated from the data rendered by the 3D models of the 74 coral fragments ( Table 1 ). For each specimen, a single measure of its surface to volume ratio (S/V), fractal dimension (FD), and sphericity (φ) captured the geometric complexity of its colony shape, the irregularity of its surface, and the compactness of its volume (complexity variables). Contrastingly, average values per branch or branch tip either estimated traits such as spacing, length, width, and angle (branching) or characterized the topological concavity/convexity of coral surfaces (curvature). A total of 19 variables were derived from the polygon meshes, 8 from the medial axis skeletons, and 26 using both the polygon meshes and the medial axis skeletons. In addition, probability density functions using kernel estimates were obtained for 15 of these variables: 8 curvature variables and 7 branching variables ( Table 1 ). Forty-one of the univariate variables did not conform one or both of the assumptions of normality and were subsequently transformed ( Supplementary Table S2 ). Four variables were removed from downstream analyses as they did not conform either to the normality or the homogeneity of variance assumption, even after transformation (i.e., K_tip_mean, K_tip_skew, k2_tip_mean, and k2_tip_skew).

To test the potential of the variables for species-level differentiation, both univariate values (i.e., global and average) and Kernel density (KD) distributions of the 3D-estimated data were examined. Using an analysis of variance (ANOVA) for univariate values, significant interspecific differences were found in the means of 42 variables when the three species previously delineated in this data set were used as a three-level factor (p-value < 0.05, df = 2; Supplementary Table S3 ). Further exploration using post-hoc Tukey tests (p-value < 0.05; Supplementary Table S4 ) indicated differences between the means of all species in 10 of the characters and between at least two pairs of species in 21 cases. In summary, more than half of the variables derived from the univariate values exhibited significant differences in the ANOVA in at least two-species comparisons in one of the two-sample tests (n = 29; preliminary selected subset).

Although no complexity or branching variable could statistically differentiate the means of all three species, at least 50% of the branching variables could differentiate two species (9/18 according to the post-hoc test). Likewise, significant differences in two to three interspecific comparisons were detected in 23 of the 33 curvature variables. In addition, although the interspecific morphospaces tended to overlap, the probability density profiles of each species were visibly distinct for most univariate variables with significant differences ( Figure 3 and Supplementary Figure S5 ). Overall, the pair of species for which more significant differences were found in the central tendencies was A vs. C, with 64% of the comparisons with p-values < α in contrast to 47% between A vs. B and 42% between B vs. C. A high degree of correlation was found between most of these variables ( Supplementary Table S5 and Supplementary Figure S8 ), which mainly corresponded to curvature features (20 curvature and 9 branching variables).

The kernel densities (KD) analysis did not reveal significant interspecific differentiation in the distribution profiles, contrasting with the univariate variables results obtained. Instead, a high degree of overlap was apparent in the distributions of most of the 3D-based morphological variables ( Figure 4A and Supplementary Figure S6 ). The cumulative analysis of the distribution comparison scores (DCS) based on the Mann-Whitney U tests of all variables and curvature variables was able to discriminate two morphogroups ( Supplementary Figure S7A and Supplementary Figure S7C , respectively). However, among the individual heatmaps plotted for each variable, there was a considerable degree of clustering between samples of the same species when assessing interspecific differences using the minimum curvature (k2) distributions ( Figure 4B ). Here, three clusters were observed, each comprising mainly individuals of one of the three species in the data set (from top to bottom: cluster 1 = 61.90% of species B, cluster 2 = 56% of species C, cluster 3 = 82.14% of species A).

Despite the overall significant differentiation in central tendency of the set of 3D-based morphological features (see section 3.2 and Supplementary Figure S9 ), none of the clusters identified by the exploratory hierarchical clustering analysis (HCA) were entirely congruent with the species delineation achieved using other lines of evidence ( Supplementary Figure S10 , see Ramírez-Portilla et al., 2022). Likewise, when feature selection for Gaussian model-based clustering allowed finding the optimal subset of features containing information, only two morphogroups were identified ( Supplementary Table S6 ). Consequently, the ordination recovered by the principal component analysis (PCA), based on such feature selection methodology, showed a considerable degree of overlap between the morphospaces as defined by the reduced dimensions considered in this analysis ( Supplementary Figure S11 ).

From the 21 variables that did not present collinearity (Variance Inflation Factor (VIF) < 10; Supplementary Table S7 ), 10 were also present in the preliminary selected subset and showed significant interspecific differentiation (MANOVA; Supplementary Table S8 ). The linear discriminant analysis (LDA) using these variables was able to distinguish three groups that matched the previously supported species delimitation with 97.30% of accuracy ( Figure 5 ; scatter plot) and was able to predict correctly at least 75% of the observations once the model was trained and later tested ( Supplementary Table S9 ). Variables such as the skewness of the Gaussian curvature (K_skew), the mean curvature (H_mean), and the branch angle (b_angle_mean) were the ones that mainly contributed to the discrimination according to the ordination coefficients of each linear component ( Figure 5 , bar plots).

While evaluating the performance of individual 3D variables, both univariate and multivariate analyses were able to identify 3D-based morphological features with significant differences between a priori delimited species. This trend was particularly evident when using curvature variables, which showed differentiation in their central tendencies in most pairwise species comparisons ( Supplementary Tables S3, S4 ), in discriminant analysis ( Figure 5 ), and partly in comparisons between Kernel density distributions ( Supplementary Figure S6 and Supplementary Figure S7C ). Indeed, significant differences in attributes such as the skewness, kurtosis, and mean values of these variables suggest that they attain sufficient resolution to capture the morphological differences between the specimens of the three complex-shaped coral taxa used here in for validation. By enabling the characterization of surface profiles and owing to the relationship between curvature variables and functional traits (Hyde et al., 1997; Ankhelyi et al., 2018), these results suggest that the estimation of curvatures holds promise for improving our understanding of the relationship between morphology and potentially specific ecological traits.

Contrastingly, analyses of branch-related variables did not provide enough resolution at the species level, likely due to the high similarity between the branching patterns of these closely related taxa (Wallace, 1999). Although only branches with a minimum of four vertices were taken into account to reduce the likelihood of including spurious ones in the estimation of branching variables, these features seemed to be highly variable within species and individuals ( Supplementary Figure S7 ). These results suggest that branch-related variables are taxonomically uninformative in this particular case, as it has been observed that species-specific patterns can emerge when comparing such variables between more distantly related taxa (Kaandorp, 1999). Besides, estimating branching variables can be more relevant to understand the function, evolution, and plasticity of an organism’s shape, particularly when studying marine taxa and their response to environmental fluctuations (Kaandorp et al., 2003; Kaandorp et al., 2005; Chindapol et al., 2013; Paz-García et al., 2015a).

Overall, analyses based on global and average univariate variables exhibited morphological differentiation consistent with a priori delineated species ( Supplementary Figures S5, S9 ). In contrast, most non-averaged density distribution analyses did not display congruent discrimination patterns ( Figure 4 and Supplementary Figure S7A ), seemingly due to the high intraspecific variation and consequent overlap of the probability density functions between a priori delineated species ( Supplementary Figure S6 ). These trends can be related to the different estimations performed in each case. The univariate values rely on average values obtained per specimen, while the probability density distributions were estimated per branch in the skeleton or per-vertex of the polygon mesh using kernel densities. Thus, since kernel densities exhibit high correspondence to data (Pradlwarter and Schuëller, 2008), they can both contain a wealth of information and display a wide range of variance that can potentially conceal the main species-specific trends in quantitative phenotypic data. Moreover, in the case of complex-shaped organisms such as the tabular Acropora corals used in this study, the variation of probability density functions can be related to the influence that shape complexity exerts on scanning reproducibility (Bythell et al., 2001). Despite the efforts to avoid the effect of self-shading ( Supplementary Materials and Methods and Supplementary Figure S1 ), it has been observed in previous studies that the coefficient of variation between iterative scans can increase in branching corals due to a higher rate of potentially overlapping structures (Reichert et al., 2016). As a result, the shape complexity of the coral specimens could have affected the variability of the observed data.

Broadly, species delimitation approaches can be differentiated into validation or discovery tools according to whether or not the samples are partitioned into taxonomic categories before performing the analysis (Carstens et al., 2013). In this study, results showed that a well-justified combination of novel 3D-based variables can aid discrimination of morphogroups of irregularly shaped organisms when based on a priori assignment of samples to categories (Ence and Carstens, 2011). In comparison to the quantitative morphological characters previously assessed from the same set of specimens by Ramírez-Portilla et al. (2022), the 3D-based variables evaluated in the current study were able to discriminate three phena congruent with other species delimitation approaches with higher overall accuracy ( Figure 5 and Supplementary Table S10 ; 97.30% in this study vs. 94.94% in the previous). However, when morphogroups were delineated without a priori information in this study ( Supplementary Figures S9–S11 and Supplementary Table S6 ), they were not congruent with species boundaries assessed using other sources of information. The high phenotypic heterogeneity detected within the a priori delineated species (particularly B and C) using density distributions ( Supplementary Figure S5 and Supplementary Table S10 ) could have hampered the unambiguous and unbiased delimitation of three phena congruent with the species boundaries previously delineated using other lines of evidence (see section 4.4).

Although these results may seem paradoxical given the significant interspecific differences found using 3D-derived variables (see section 4.1), phenotypic differentiation in central tendencies between species defined a priori does not necessarily count as evidence of species boundaries when assessed in light of evolutionary theory (Luckow, 1995; Zapata and Jiménez, 2012; Cadena et al., 2018; Cadena and Zapata, 2021). Instead, distinct distributions of phenotypic characters (e.g., those derived from fitting quantitative data to NMMs) can constitute support for species hypotheses as long as they do not result from intraspecific polymorphisms (e.g., González-Espinosa et al., 2018) or morphological plasticity (e.g., Paz-García et al., 2015a). Here, this trend only became evident once the features that yielded information congruent with a priori species boundaries were selected and used to perform discriminant analyses ( Figure 5 and Supplementary Tables S9, S10 ).

The quandary of feature selection in multidimensional data sets is often the limiting factor for extending the applicability of approaches such as 3D-derived variables to a wider variety of organisms (Poon et al., 2013). Certainly, many issues in detecting species boundaries from morphological and phenotypic analyses derive from the potential exclusion of important characters during dimensionality reduction (Cadena et al., 2018). Here, given the large number of variables derived from the 3D analyses and the fact that not all of them provided discriminative and non-redundant information ( Supplementary Tables S3–S7 ), the process of feature screening proved to be key to selecting features used to discriminate between a priori delimited species (Ramírez-Portilla et al., 2022). Feature screening and selection, however, would substantially rely on the overall morphology of the studied organisms. Assessing branching features, for example, would be inadequate for describing the 3D morphology of massive or encrusting growth forms.

Although the morphospaces overlapped when performing bivariate comparisons and other multivariate graphical representations ( Figure 3 and Supplementary Figure S5 ), linear combinations of features after variable screening were useful to identify potentially informative characters and a combination of them that enabled discrimination of three morphogroups congruent with the species delineated a priori in the data set ( Figure 5 ). These results support the notion that not only technological advances in 3D data acquisition and model rendering, but also feature screening and selection actually provide prospective variables to quantify morphology and discriminate groups (Valcárcel and Vargas, 2010), particularly of complex-shaped organisms lacking traditional landmarks.

Regardless, the methodology used here to estimate morphological variables from the 3D models, can be applied to understand a wider variety of phenomena such as morphological plasticity, development, and environmental effects on shape and biodiversity. Therefore, the approaches implemented in this study do not only intend to inform taxonomy, but also to provide tools that can support evolutionary, ecological, and biomonitoring aims to characterize and understand form in complex-shaped taxa in forthcoming studies.

Independent lines of evidence that previously delineated taxonomic units in the data set (i.e., morphology, breeding trials, and molecular analyses) robustly supported the identification of three species (Ramírez-Portilla et al., 2022). Although the 3D-based variables assessed here provided enough power to discriminate morphogroups congruent with such species boundaries delineation, it was not able to delimitate the same three groups when unbiased clustering and feature selection were performed. Only two clusters or components were identified and supported by the multivariate analyses ( Supplementary Figures S9, S10 ), even after a set of variables that best discriminated groups using normal mixture models (NMMs) without a priori information was employed ( Supplementary Table S6 ). In this regard, the results obtained from the density distributions suggest that the heterogeneity of the phenotypic variation detected by 3D-morphological analyses within some of the species could have confounded the identification of components in the mixture ( Supplementary Figure S5 and Supplementary Table S10 ). This would be consistent with the close ties between intraspecific variability and interspecific similarity that have hampered the widespread use of morphological features to delineate taxonomically intricate taxa (Sites and Marshall, 2004), particularly in speciose groups such as the coral genus Acropora where both high morphological intraspecific variability and interspecific similarity, particularly between closely related species, has been reported (Wallace and Willis, 1994; Wallace, 1999).

Alternatively, the scanning quality achieved in the present study could have hindered the potential of species-level delimitation given that important differences in the microstructure could be masked by the technical resolution of the underlying 3D mesh (Gutierrez-Heredia et al., 2016; Reichert et al., 2017). Indeed, features at micromorphological level, such as corallite shape and dimension, have been deemed crucial skeletal characters for discriminating between complex-shaped coral species like those of the genus Acropora (Wallace, 1999; Wolstenholme et al., 2003; Ramírez-Portilla et al., 2022). The potential masking of these features and the effect of the sample sizes used for validation, which in several cases were lower than n = 10 (Finch and Schneider, 2006), could explain the relatively low prediction accuracy achieved in this study when randomly partitioning the data set in training and testing subsets ( Supplementary Tables S9, S10 ). These results suggest that further refinements of 3D-morphological analyses, such as increased precision and reduced errors, might even enable a priori delimitation of phena rather than a posteriori confirmation. For instance, the integration of 3D models from structured light scanning or photogrammetric approaches with high-resolution 3D methodologies such as CT scanning could solve some of the present accuracy issues (Laforsch et al., 2008; Naumann et al., 2009; Veal et al., 2010; Gutiérrez-Heredia et al., 2015; Aston et al., 2022). In the meantime, our results support the discriminative value of implementing 3D-based variables either as hypotheses testing or validation approaches rather than discovery ones.

Progress in the development of methods to delineate species using morphological data is urgently needed, particularly to improve modelling of phenotypic variation in agreement with evolutionary theory (Cadena and Zapata, 2021). Due to the potential for morphological plasticity of marine taxa such as corals (Todd, 2008), assessing the discriminative power of 3D- morphological variables to distinguish species throughout the environmental ranges they occupy ought to be explored further, particularly for species with relatively wide distributions where geographical differences can mislead morphology-based species delimitation (Fukami et al., 2004; Forsman et al., 2015). Moreover, the link between 3D-based phenotypic variables and their biological, ecological, and functional significance need to be addressed in future studies that will not only intend to assess interspecific variations in morphology and their taxonomic relevance, but also their potential role in speciation and adaptation, particularly for complex-shaped organisms such as corals (Zawada et al., 2019a; Zawada et al., 2019b; Torres-Pulliza et al., 2020; Aston et al., 2022; Siqueira et al., 2022), which inhabit one of the ecosystems most threatened by climate and anthropogenic disturbances (Hughes et al., 2017; Hughes et al., 2018).