*Article* **Artificial Intelligence as a Tool to Study the 3D Skeletal Architecture in Newly Settled Coral Recruits: Insights into the Effects of Ocean Acidification on Coral Biomineralization**

**Federica Scucchia 1,2,\*, Katrein Sauer 3, Paul Zaslansky 3,\*,† and Tali Mass 1,4,†**


**Abstract:** Understanding the formation of the coral skeleton has been a common subject uniting various marine and materials study fields. Two main regions dominate coral skeleton growth: Rapid Accretion Deposits (RADs) and Thickening Deposits (TDs). These have been extensively characterized at the 2D level, but their 3D characteristics are still poorly described. Here, we present an innovative approach to combine synchrotron phase contrast-enhanced microCT (PCE-CT) with artificial intelligence (AI) to explore the 3D architecture of RADs and TDs within the coral skeleton. As a reference study system, we used recruits of the stony coral *Stylophora pistillata* from the Red Sea, grown under both natural and simulated ocean acidification conditions. We thus studied the recruit's skeleton under both regular and morphologically-altered acidic conditions. By imaging the corals with PCE-CT, we revealed the interwoven morphologies of RADs and TDs. Deep-learning neural networks were invoked to explore AI segmentation of these regions, to overcome limitations of common segmentation techniques. This analysis yielded highly-detailed 3D information about the RAD's and TD's architecture. Our results demonstrate how AI can be used as a powerful tool to obtain 3D data essential for studying coral biomineralization and for exploring the effects of environmental change on coral growth.

**Keywords:** coral reefs; coral recruits; biomineralization; skeletal structure; synchrotron phase contrast-enhanced microCT; PCE-CT; artificial intelligence; ocean acidification

### **1. Introduction**

In tropical and subtropical oceans, symbiotic stony corals provide an ecological framework that retains nutrients, support high rates of autotrophic production, and harbor extensive biological diversity. The skeletons, created by tiny coral polyps in colonies, are known to be among the largest bioconstructions in the world. Understanding the details of structure generation and principles of coral skeletal formation have been a common topic among diverse fields of study including biology, geochemistry, geology, paleontology, and materials science. The skeletons of stony corals consist of two main skeletal regions known as Rapid Accretion Deposits (RADs; also known as "calcification centers") and Thickening Deposits (TDs; also called "fibers"). RADs correspond to areas of the skeleton with rapid calcium carbonate deposition where new skeletal growth zones start to form [1–3]. They contain an organic matrix and consist of granular aggregates that comprise randomly oriented calcium carbonate crystals containing both aragonite and stable amorphous calcium carbonate [4–6]. A second calcium carbonate form, Mg-calcite, has been shown to also be present in the early stages of RADs formation, typical of initial mineral development

**Citation:** Scucchia, F.; Sauer, K.; Zaslansky, P.; Mass, T. Artificial Intelligence as a Tool to Study the 3D Skeletal Architecture in Newly Settled Coral Recruits: Insights into the Effects of Ocean Acidification on Coral Biomineralization. *J. Mar. Sci. Eng.* **2022**, *10*, 391. https://doi.org/ 10.3390/jmse10030391

Academic Editors: Hildegard Westphal, Justin Ries and Steve Doo

Received: 13 January 2022 Accepted: 5 March 2022 Published: 9 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

in early coral life stages [7]. Outward from RADs, the mineral grows into TDs, forming elongated aragonite crystals with a fibrous morphology appearing compact and dense with significantly reduced amounts of organics as compared to RADs [2,3,8–11].

The coral skeleton has been well characterized at the micro-structural level [2–5,8–13], but less is known about how changes in the environment may affect its development, in particular regarding the formation of RADs and TDs. Recent predictions suggest that future ocean acidification (OA) conditions, resulting in a chemical shift toward higher seawater pCO2 and lower pH, will significantly hinder the coral mineral formation and overall skeletal growth [14–18]. In the early stages of coral life, we have recently shown that exposure to future OA conditions causes a reduction in the abundance of RADs in the skeletal spines [19]. However, previous results relied solely on scanning electron microscope (SEM) images; hence, they were limited only to two-dimensional (2D) surface observations. Indeed, most studies on coral skeletal structures make use of 2D-based microscopy analyses [1–3,8,9,11–13,20], which lack volumetric quantifications. Therefore, the lack of an appropriate methodology to precisely visualize the volumetric configuration of RADs and TDs has limited investigation of the 3D structure of these two skeletal regions [3,11,21].

Following the widespread development of high-resolution X-ray μCT in medical and industrial applications, non-destructive imaging of 3D structures has become popular. It is thus possible to monitor mineral internal structural changes [22] and characterization [23] at resolutions down to a few tens of nanometers. Conventional X-ray μCT imaging methods mainly rely on X-ray absorption as a mechanism to generate images based on the attenuation of the X-rays by the material, yielding spatial information about density distributions within the sample. This technique, however, is prone to normalization artifacts and does not distinguish between materials with very similar X-ray attenuation [24,25]. Other methods such as phase contrast-enhanced microCT (PCE-CT) exploit differences in X-ray interaction to create high-contrast images where even small differences in density along the X-ray path can lead to strong contrast enhancement at interfaces between material phases [24]. PCE-CT imaging can thus highlight edges and internal boundaries in intact samples, even in the case of quite similar material compositions [24]. These characteristics of PCE-CT make it a good candidate for enhancing the visibility of RADs and TDs and ease their identification within the coral skeleton.

For a detailed analysis of microstructures, volumetric quantification techniques must be applied with the capacity to accurately differentiate and represent individual material phases. In recent years, several works have shown that artificial intelligence (AI) algorithms can be applied for visual classification tasks such as in the case of materials with complex microstructures [26–28]. These applications make increasing use of machine learning methods based on Deep Learning (DL) neural networks. DL networks comprise a large number of interconnected processors (neurons) that work in parallel [29]. For visual classification, the network is trained to identify structures of interest using a training "ground truth" dataset, which is used by the computer as a reference to classify similar images. The process is known as supervised classification, and has been employed in μCT images of minerals to visualize and characterize 3D structures with high accuracy [30–33].

In the present study, we combine AI and PCE-CT to obtain quantitative information about the density and 3D distribution of RADs and TDs within coral skeletons imaged with phase contrast enhancement. Such information is currently unavailable using conventional 2D or other imaging methods. We used coral primary polyps as the study system for the development of a robust analytic 3D quantification methodology, since rapid calcification during initial life stages of stony corals provides a unique opportunity to study the formation of calcium carbonate that involves extensive skeleton morphological changes. In addition, the early stages of coral development represent a good candidate to test the performance of our approach under experimental conditions that we know induce alterations in the development of RADs, such as exposure to OA conditions [19].

To develop our new methodology and optimize its performance across varied conditions, we analyzed high resolution scans of the skeletal structures of recruits of an abundant stony coral species *Stylophora pistillata* from the Red Sea, following settlement and growth under both natural and OA conditions. In particular, we simulated the global mean surfaceocean decline in pH predicted to occur under the RCP8.5 scenario [34] by the end of this century (pH drop ca. 0.1–0.4 units, from ~8.0 to 7.6; [35]). Employing PCE-CT, we identified the μm-sized edges of RADs and TDs within the skeleton. The high-contrast images generated with PCE-CT allowed us to use AI to segment RADs and TDs, under both natural conditions and under what we expect to be morphologically-altered OA conditions, and to reproduce these structures in 3D with quantitative detail.

#### **2. Materials and Methods**

#### *2.1. Sample Collection and OA Experiment*

Coral larvae were collected from 15 randomly selected adult colonies of the stony coral *S. pistillata* on the reef adjacent to the Interuniversity Institute of Marine Sciences (IUI, 29◦30 06.0" N 34◦54 58.3" E) in the Gulf of Eilat (Israel), under a permit from the Israeli Natural Parks Authority. Larvae were collected using dedicated traps (160 μm plankton net with a plastic collection container) that were placed on adult colonies from the shallow reef (depth 6–7 m) for several nights of spawning during April 2020.

All actively swimming larvae were pooled together and were transported to a controlled environment aquarium system at the Leon H. Charney School of Marine Science (University of Haifa, Haifa, Israel). Larvae were placed in custom-made polypropylene plastic chambers (~20 larvae per chamber) consisting of a central cylinder sealed at the ends by plankton netting. Settlement chambers were used in this experiment to avoid losing larvae through the water recycle system of the aquariums before settlement, as previously reported [19,36,37]. Prior to the insertion of the larvae, the chambers were washed to remove any potential chemical released by the plastic, and they were left to soak in seawater for one week.

The settlement chambers (1 chamber per aquarium) were placed in a system of 6 flowthrough aquariums with artificial seawater (Red Sea Salt, Red Sea Ltd., Houston, TX, USA) replicating the spring northern Red Sea water conditions as previously described [19]: pH of 8.19 ± 0.01 (seasonal mean ± s.d.), salinity of 40.63 ± 0.03 g L−<sup>1</sup> (seasonal mean ± s.d.), and temperature of 22.18 ± 0.17 ◦C (seasonal mean ± s.d.). An irradiance of 250 μmol photons m−<sup>2</sup> s−<sup>1</sup> on a 12-h light-dark cycle was provided by a Mitras LX 7206 LED aquarium light system (GHL, Kaiserslautern, Germany).

Prior to the insertion of the settlement chambers, the carbonate chemistry of seawater was manipulated in 3 of the experimental aquariums by injecting CO2 to reduce the ambient pH 8.2 (pCO2~487 μatm) and obtain the target value of pH 7.6 (pCO2~1938 μatm). Temperature and salinity were monitored continuously throughout the experiment by electrodes linked to a monitoring system (GHL, Kaiserslautern, Germany) that controlled CO2 injection, and they were also measured three times per day using a portable Orion Star™ A222 conductivity meter (Thermo Fisher Scientific, Waltham, MA, USA). Measurements of pH (NBS scale) were performed three times per day using a pH glass electrode (Metrohm, Herisau, Switzerland) where buffer solutions (Rocker Scientific, Taiwan) were employed to perform calibration (pH 4, pH 7, and pH 10). Measurements of total alkalinity (TA) were made once per day (triplicates) via titration with 0.1 N HCl containing 40.7 g NaCl L<sup>−</sup>1, using an automatic alkalinity titrator (855 Robotic Titrosampler, Metrohm, Herisau, Switzerland) controlled by Tiamo (Software version 2.0, Metrohm, Herisau, Switzerland). Automated titrations of 50 mL samples were performed. Parameters of seawater carbonate system were calculated from pH, TA, temperature, and salinity using the CO2SYS package [38] with constants from [39] as refit by [40] (Supplemental Table S1). Throughout the experiment, corals were fed once per day with 2 mL of concentrated planktonic suspension (Microvore, Brightwell R aquatics, Fort Payne, AL, USA). The experiment lasted for a total

of 9 days, following our previous experimental procedure [19], after which primary polyps were gently removed from the chambers and stored with 90% ethanol in 2 mL tubes.

#### *2.2. X-ray μCT: Image Acquisition and Tomographic Reconstruction*

Primary polyps were fixed with EpoFix resin (Agar Scientific, Stansted, UK) on top of polypropylene micropipette tips. Tomographic scanning was conducted at BAMline [41,42], the imaging beamline of BESSY II (the synchrotron storage ring of HZB—Helmholtz-Zentrum Berlin, Germany). Each sample (N = 3 per pH treatment) was attached to a metal stub and scanned with incremental rotation (multiple projections spanning 180◦) (Figure 1a) using the high-resolution imaging setup (Supplemental Figure S1) [43] with exposure times set to 1 s. Projection images were acquired with a final pixel size of 2.2 μm using PCE-CT imaging mode at an energy of 24.5 keV. Additional absorption scans were performed for comparison, using an energy of 16.26 keV, to ensure correct identification of the skeletal structures avoiding contrast enhancement.

**Figure 1.** Workflow for acquisition and processing of tomographic datasets of coral recruit skeletons. (**a**) X-ray projections are acquired from different angles of the skeleton with incremental rotation (0◦ to 180◦); (**b**) following normalization (gray star = additional details are shown in Figure 2), tomographic cross-sectional slices are reconstructed; (**c**) cross-sectional slices are stacked and processed into 3D views of the recruit skeletons; (**d**) AI-based segmentation and classification are performed so that the features of interest can be analyzed and visualized in 3D. An example distribution of RAD thicknesses is shown.

Prior to reconstruction, data were normalized to account for beam inhomogeneities using an in-house Octave-based reconstruction pipeline in the labs of the Charité, Universitätsmedizin (Berlin, Germany). Specifically, for each scan, the radiograms were background-corrected by normalization with best-fitting (minimum variance) empty beam (flat-field) images, obtained both before and after each scan, and corrected by subtraction of dark-current images (Figure 2). Reconstruction was performed by the filtered back projection method using nRecon (v 1.7.4.2, Brucker micro-CT, Kontich, Belgium) to generate 3D datasets from cross-sectional 2D images (Figure 1b). Tomographic datasets were visualized and further processed in 3D using Dragonfly (v 2021.3, Object Research Systems—ORS, Montreal, Quebec, Canada) (Figure 1c,d) [44], with additional stack analyses performed in FIJI [45].

**Figure 2.** Example X-ray μCT radiographs of a recruit skeleton as obtained on the beamline (**left**) and after normalization (**right**). Radiographs were background-corrected by normalization with empty beam (flat-field) images and by the subtraction of dark-current images. The dark silhouette of the coral is seen to be suspended within a faint halo of epoxy resin, visible in the normalized (**right**) image.

#### *2.3. Deep-Learning-Based Image Segmentation*

To objectively identify structures in the PCE tomographic data, each reconstruction was segmented using the Deep Learning functions built into the 3D commercially-available software Dragonfly (v2021.3, Object Research Systems—ORS, Montreal, Quebec, Canada) [46]. Segmentation is necessary for the processes of analysis and quantification of the volume of RADs and TDs with respect to the total volume of the skeleton (mineral excluding air cavities) of each primary coral polyp scanned. This process is required to partition digital images into several sub-volumes, each representing a distinct 3D feature, obtained by identification and grouping of similar volume pixels (voxels) [47]. Image segmentation is thus typically used to locate objects of interest in large volumes of data needing quantification. In the present study, we applied semantic image segmentation by labeling and grouping pixels in each slice-image according to the feature that they were part of (RADs, TDs, or air spaces). For this classification, we used the U-Net network architecture algorithm, which has been shown to be very efficient in segmenting biomedical images [48]. When applying supervised classification, the network is manually trained by providing examples of how different features within the data should be classified. This is achieved by manually selecting groups of voxels and assigning them to a specific feature within a set of representative images, considered to be a training dataset. The computer then uses the manually-segmented training datasets as reference for automatic classification of all other regions in all images. We therefore manually delineated RADs, TDs, and air spaces in several (<10) 2D slices within each reconstructed tomographic dataset. These were used to train the network to recognize different features within all other images (Figure 3a–d).

**Figure 3.** Overview of the DL (deep-learning-based) image segmentation process. (**a**) Example tomographic cross-sectional slices across the skeleton of a primary polyp before segmentation, revealing RADs (dark gray areas, yellow arrows) enclosed within TDs (light gray areas, pink arrowheads). (**b**,**c**) Manually delineated RADs (**b**; marked in yellow) and TDs (**c**; marked in pink) defined in a confined region within the training dataset. (**d**) Complete region manually delineated as training dataset, showing all three features (RADs in yellow, TDs in pink, and air spaces in green). (**e**) Result of the complete dataset processed by automatic DL-based segmentation. The AI network was applied to all slices within the stack, and it has automatically successfully distinguished RADs (yellow) and TDs (pink). Scale bars: 300 μm.

Training requires optimizing a series of parameters that define the performance of the U-Net network workflow. Within the Dragonfly platform, several parameters can be tuned: (1) patch size (images are split into smaller rectangular patches that contain the feature of interest [49]), (2) stride-to-input ratio (a parameter defining the typical position relationships between neighboring patches), (3) batch size (the number of patches evaluated concomitantly), (4) number of epochs (cycles of iteration that involve computation of all batches of the training dataset), (5) loss function (which measures the difference between groups of voxels in the training set images, considered to be "the ground truth images", and the output value predicted by the network [50]), (6) optimization algorithm (different approaches to minimize the loss function, and thus to reduce the training error), and (7) data augmentation (how many times each data patch is augmented during a single training epoch). In this work, optimization of the DL network performance was achieved by making trial runs with different combinations of number of training slices and of the above-mentioned values of network parameters (starting with default settings provided at first by the software), until reaching reproducible accuracy of the segmentation (a single user ensured by visual examination that the trained network accurately segmented skeletal regions, by comparison to the ground truth images, i.e., the training set of images). The final optimized workflow thus required training datasets of 8 to 11 slices and a series of parameters (1–6 as described above) set to: (1) 32, (2) 1, (3) 32, (4) 100, (5) OrsDiceLoss, (6) Adadelta, and (7) 2. The trained network was then applied to segment all slices within each tomographic stack (Figure 3e). The segmented skeletal structures were visualized in 3D (Figure 1d) and their thickness was computed.

#### *2.4. Performance of the DL-Based Image Segmentation and Evaluation of the RADs/TDs Shape Variability within and between Tomographic Datasets*

Each tomographic dataset (corresponding to a single coral recruit) was used to train a DL network independently. Representative slices from each dataset were trained to identify the 3 image features, as shown in Figure 3. Thereafter, networks optimized for each dataset were checked for segmentation of other datasets (e.g., DL network optimized on dataset A applied to segment dataset B). The degree to which image features were correctly classified as TDs or RADs was assessed by visual inspection and by computing the degree of mismatch (underestimation or overestimation in % of the RADs area with respect to the total area of the skeleton - measured as the combined cross-sectional RAD + TD areas in each slice). Performance was also compared with a DL network trained using all tomographic datasets (3 from the control pH and 3 from the low pH, 8–11 training slices were selected within each dataset) and using it to segment datasets from both the control and low pH conditions.

Because the performance of AI-based image segmentation is strongly influenced by the degree of variability of the target-objects shapes [51–54], we estimated the variation in the RADs-TDs relative shape configuration by using basic shape descriptors (perimeter and area per slice) and comparing them across datasets. In particular, we computed the RADs-TDs relative shape configuration according to Equation (1).

$$C = \frac{\frac{P(r)}{A(r)}}{P(t)/A(t)}\tag{1}$$

where *P*(*r*) and *A*(*r*) are the perimeter and area of RADs, respectively, and *P*(*t*) and *A*(*t*) are the perimeter and area of TDs, respectively. Such values were computed per each slice within each single dataset and ultimately expressed as a percentage. We then calculated the skewness of C distributions and used it as a metric to evaluate the degree of the RADs-TDs shape variation across datasets, as well as to compare it with the results of the segmentation performance.

#### *2.5. Statistical Analysis*

The per-slice percentage of RADs area over the cross-sectional area of the skeleton computed for different segmentation outcomes, as well as the skewness of the distribution of the per-slice RADs-TDs configuration, were tested for normality (Shapiro–Wilk test) and homogeneity of variance (Brown–Forsythe test). For the skewness measurements, the unpaired *t*-test was used (to test the differences between pH treatments), in which significant groups have a value of *p* ≤ 0.05 (*n* = 3 primary polyps per pH condition). Nonparametric equivalents of tests were used in cases where the assumptions of normality and homogeneity of variance were violated. For such cases (i.e., mean per-slice percentage of RADs area over the total area of the skeleton), a Wilcoxon matched-pairs signed rank test was used, in which significant groups have a value of *p* ≤ 0.05. The GraphPad Prism software version 9.0.0 (GraphPad Inc., San Diego, CA, USA) was used to perform the statistical tests. All results are presented as mean ± SE.

#### **3. Results**

#### *3.1. Detection and Distinction of Coral Skeletal Structures*

We reconstructed high-resolution tomography datasets of up to 1000 cross-sectional slices per scanned coral primary polyp, each amounting to a total of 54 GB of data. Each reconstructed image presents the complete form of the original coral skeleton including details of the internal skeletal architecture. In the PCE-CT images, strong differences in contrasts within the skeleton revealed clear internal variations in mineral distribution (Figure 4). Such variations coincide with previous observations of RADs possessing a less compact arrangement of the aragonite aggregates and a higher amount of organics compared to TDs [2,3,8–13]. Absorption images further confirm variations in mineral density between skeletal regions, though with blurry edges due to the poor signal. These scans show that RADs are indeed less dense than the surrounding TDs (~20% lower density; Supplemental Figure S2). This clear difference in density allowed us to ascertain that the contrasts observed by PCE-CT correspond to RADs and TDs in all phase contrast enhanced datasets. To allow for volumetric quantification and 3D visualization, image segmentation has to be applied, which requires clear boundaries to be present between different features [55]. The PCE-CT imaging mode appears, therefore, to be perfectly suited for the purpose, as the interfaces between RADs and TDs are accentuated and easy to identify, especially when compared to absorption images (Figure 4).

**Figure 4.** *Cont.*

**Figure 4.** Typical coral skeleton details observed in cross-sectional tomography slices and comparison between absorption and PCE imaging modes. Differences in X-ray absorption within the skeleton reveal variations in both mineral distribution and density between RADs (darker regions correspond to lower density) and TDs (lighter regions, higher density); the red marked insets are shown on the right. Strongly defined boundaries can be observed between RADs and TDs in slices from PCE-CT datasets, whereas absorption images provide a blurry, poorly distinguishable separation between the two skeletal regions. Scale bars: 300 μm (**left** images), 100 μm (**right** images).

#### *3.2. Segmentation of RADs and TDs*

Conventional segmentation methods typically rely on choosing cut-off gray-values and are often based on the voxel intensity histograms. By applying a thresholding, the image can be segmented into the features that are below or above defined threshold values [56]. A gray-value histogram of an example coral skeleton shows that no clear distinction exists in gray-values corresponding TDs and RADs, as there is an overlap between the intensity ranges (for example, see gray-values of the margin of TDs facing the exterior of the skeleton and facing the RADs in Figure 5a,c). Phase-contrast images in fact are edge enhanced, meaning that both the brightest and lowest intensities of the image are found in the margins delineating the edges and internal boundaries of different densities and interfaces in the sample, and meaning that different features of interest may have the same "color".

Although PCE-CT imaging creates a sharp visual discrimination of RADs and TDs, a simple "binary" threshold applied to this type of image does not yield precise segmentation and reliable distinction between the two skeletal regions, due to the overlap between the intensity range of RADs and the TDs margins. This is shown by the mismatches between the distribution of RADs and TDs as observed in the reconstructed image and the distribution outlined by simple thresholding (Figure 5a,b). Thus, segmentation and ultimately quantification and 3D visualization of RADs and TDs require the use of a more sophisticated technique. Figure 6a depicts results obtained by the automatic segmentation performed using DL networks, which reveals a far better distinction of RADs and TDs within the recruit skeletons. The impressive separation between features becomes evident within the complex shapes of RADs across the skeleton (Video S1), and provides a means for overcoming the limitations imposed by the overlap in gray-values intensity ranges. Such detail-sensitive segmentation of RADs and TDs across an entire stack of slices ultimately allows for accurate volumetric quantification and 3D reproduction of the entire recruit's skeleton, including fine details of local thickness within RADs (Figure 6b, Video S1).

**Figure 5.** Slices and typical histogram of gray-values in PCE image and quality of non-DL-based image segmentation. (**a**) Close-up of the skeleton, showing RADs (dark gray areas), TDs (light gray areas), and delineating the outer and inner edges of TDs (green and pink lines, respectively). (**b**) Result of a simple threshold applied to segment RADs (RADs segmented by the thresholding method are colored in red). Note that, compared to their actual shape in the ground truth image (**a**), RADs are not precisely segmented (see arrows for example mismatch, where TDs regions are classified as RADs) and they are not distinguished from the outer edge of TDs, which are also thresholded as RADs by this method. (**c**) Gray-value histogram for an example tomographic stack of an imaged coral skeleton. The gray-value ranges for all regions depicted in (**a**) are highlighted in different colors. Note the overlap between the gray-value ranges of RADs and of the inner/outer edges of TDs, which makes it impossible to reliably segment these regions by simple threshold, as demonstrated in (**b**).

**Figure 6.** Deep-Learning-based image segmentation and quantitative 3D reproduction of segmented RADs and TDs. (**a**) Two example 3D renderings of a primary polyp skeleton virtually sliced in the transverse plane (upper row) and in the longitudinal plane (bottom row). Red dotted lines show enlargements of the sliced skeleton that have been segmented using DL (segmented TDs and RADs are colored in pink and yellow, respectively). DL, Deep Learning. (**b**) Segmentation of RADs and TDs across the entire stack of slices allows 3D separation of features and quantification, including details of local thickness (μm).

#### *3.3. Evaluation of the DL-Based Segmentation Performance*

Figures 6 and 7c,d show the intricate details made available due to the DL segmentation based on training an independent network (e.g., network A trained and applied on

dataset A). Visual inspections of the segmentation outcomes reveal that applying a network optimized on a single tomographic dataset does not always provide a precise segmentation of other datasets (cross segmentation), as demonstrated by small mismatches between the real arrangement of RADs and TDs (as observed in "ground truth" images, Figure 7a,b) and the distribution outlined by the DL segmentation (Figure 7e,f).

**Figure 7.** Evaluation of the DL-based segmentation performance. (**a**,**b**) Details of the skeletal structures in representative coral recruits from the control and low pH conditions before applying the DL-based segmentation (ground truth images). RADs and TDs in example areas have been manually delineated in yellow and pink, respectively. (**c**,**d**) Segmentation of RADs and TDs based on training a network independently per each single dataset. Specifically, (**c**) shows the outcome of the segmentation of network A trained and applied on dataset A (a representative control pH sample), whereas (**d**) shows the outcome of the segmentation of network X trained and applied on dataset X (a representative low pH sample). (**e**) Segmentation obtained by employing a DL network that was trained using dataset B (another representative control pH sample) and applied to dataset A. The white arrows point to mismatches in the segmentation as compared to the more accurate segmentation in (**c**). Note how RAD regions have been segmented as TDs. (**f**) Segmentation obtained by training the DL network on dataset A and using it to segment dataset X. The white arrows point to mismatches in the segmentation as compared to the correct segmentation in (**d**). (**g**,**h**) Segmentation obtained by employing a DL network that was trained using all tomographic datasets (all 6 datasets from both control and pH conditions). White arrows identify example points where different segmentation is seen in (**g**) as compared to (**e**), and the better segmentation in (**h**) as compared to (**f**). Scale bars: 100 μm.

Such visual segmentation mismatches for the cross segmentation cases lead to significant differences in the mean per-slice percentage of RADs area over the total area of the skeleton, for both control pH dataset (mean per-slice difference of −1.60% ± 0.07 between the % RADs area/total area obtained with independent segmentation and the one obtained with cross segmentation; Wilcoxon test, W = −20060, *p* < 0.0001) (red dots compared to green dots in Figure 8a) and for low pH dataset (mean per-slice difference of +7.90% ± 0.38 between the % RADs area/total area obtained with independent segmentation and the one obtained with cross segmentation; Wilcoxon test, W = 13861, *p* < 0.0001) (purple dots compared to blue dots in Figure 8b). Importantly, however, for the normal pH, the differences are less than 2%, which is impressive given the huge number of slices processed. The generation of such discrepancies may be due to the sensitivity of AI in segmenting images where there are large variations in the shape of the objects of interest between different

datasets [51–53]. We indeed found substantial variations in the shape of skeletal features among datasets, as indicated by the different skewness values of the distributions of the RADs-TDs shape configurations (Figure 8c,d and Supplemental Figure S3). Curiously though, we did not find a significant correlation between these skewness values and the degree of overestimation or underestimation of the RADs/total skeleton area. This suggests that, on average, DL segmentation does not skew the data and is certainly no less accurate than other segmentation methods.

**Figure 8.** Evaluation of the DL-based segmentation performance based on the per-slice estimation of the area of the coral skeletal features and skewness of the distributions of the RADs-TDs shape configuration across samples. (**a**,**b**) Measured per-slice area of RADs over the total area of the skeleton (RADs + TDs) estimated using independent, cross, and all-inclusive segmentations (as shown in Figure 7) in the case of a control pH sample (**a**) and a low pH sample (**b**). The slice number indicates the progression through the recruit's skeleton going from the base to the top, as illustrated by the coral skeletons on the right side of the graphs. Note how, for both control and low pH datasets, the per-slice area of the RADs obtained with all-inclusive segmentation gets closer to the value obtained with independent segmentation (red dots versus yellow dots and purple dots versus light blue dots, respectively) compared to the per-slice area obtained with cross segmentation (red dots versus green dots and purple dots versus blue dots, respectively). (**c**,**d**) Skewness of the distributions of the per-slice RADs-TDs shape configuration in control pH (**c**) and low pH (**d**) samples (see also Supplemental Figure S3). The ratio of the RADs perimeter/area over the TDs perimeter/area ("C") was computed as an estimate of the relative shape of the skeletal features across the recruit skeletons. The skewness values are used as a metric to evaluate the degree of the RADs-TDs shape variation across datasets.

Notably, differences in the skewness of distributions are particularly evident when comparing datasets across pH conditions. Low pH datasets in fact have a significantly lower mean skewness (5.47 ± 0.48) compared to control pH dataset (mean of 11.13 ± 1.36; Unpaired *t* test, *t* = 3.922, *p* = 0.017) (Figure 8c,d), suggesting an underlying shaping effect of acidic seawater on the RADs-TDs distribution across the recruits skeleton, requiring further in-depth investigation.

Further visual evaluation of the different segmentation outcomes reveals that applying a network trained using all tomographic datasets (named all-inclusive segmentation) provides a much improved segmentation performance compared to the cross segmentation (Figure 7e–h). These observations are supported by the measurements of the mean perslice percentage of RADs area over the total area of the skeleton (Figure 8a,b). For both control and low pH datasets, the all-inclusive segmentation yielded a significantly better estimation of such percent area than the one produced by the cross segmentation (for the control pH dataset: mean per-slice difference of −0.95% ± 0.10 between the % RADs area/total area obtained with independent segmentation and the one obtained with allinclusive segmentation, Wilcoxon test, W = −13081, *p* < 0.0001; for the low pH dataset: mean per-slice difference of +2.65% ± 0.15 between the % RADs area/total area obtained with independent segmentation and the one obtained with all-inclusive segmentation, Wilcoxon test, W = −13651, *p* < 0.0001)(see red dots compared to yellow dots in Figure 8a and purple dots compared to light blue dots in Figure 8b). These results indicate that variations in the coral geometries themselves increase the variability of shapes encountered by the model during training, which in turn substantially increases the robustness of the DL-based segmentation performance.

#### **4. Discussion**

Combining PCE-CT and AI-based image segmentation, we have developed a highly sensitive methodology to obtain quantitative volumetric information about RADs and TDs, the two main skeleton domains of reef-building corals (Figure 6, Video S1). As is the case for many tomographic analyses, image segmentation of the coral is a necessary step for quantitative analytical investigations of the datasets. Our work provides a pipeline that allows to extract the regions corresponding to volumes and densities of interest and to obtain their detailed 3D geometry quantification. Image contrast differences strongly influence the segmentation performance, as they determine how well objects of interests are differentiated by both human and machine vision, thus improving the object perceptual quality [57,58]. In the present study, we show that PCE-CT yields the needed high contrast difference that is essential to visually identify the fine intricate interwoven structures of TDs and RADs within the intact coral skeleton. Indeed, PCE has proven to be a powerful tomographic method in both biomedical research and geological sciences to increase image contrast and enhance image segmentation performance [59].

Although PCE-CT imaging greatly facilitates the visual distinction between RADs and TDs by accentuating their μm-sized boundaries, it introduces difficulties to precisely segment the two skeletal regions using simple segmentation techniques. Specifically, edge enhancement and strong artefacts totally preclude any phase retrieval such as the paganin method [60]. For many datasets, thresholding is the simplest and more common method of segmenting images, which assumes that it is possible to find a single threshold value coinciding with the grayscale intensity of the pixels that correspond to the same structure in the volumetric image [56]. However, accurate image thresholding is only possible in cases where the different features in an image have clearly different peaks in the gray-value histogram. The significant overlap observed between the gray-value ranges of RADs and of the TDs edges makes it impossible to achieve a complete segmentation of the skeletal structures using a simple threshold (Figure 5). It is for this reason that DL algorithms are of great interest, as they may help overcome this difficulty by iteratively learning from data, allowing computers to digitally dissect images through automatic detection and classification of pixels belonging to the structures of interest [61]. In spite of the high complexity of the RADs structure across the recruit skeletons, the automatic segmentation performed by training DL networks delineated and distinguished RADs and TDs with great detail, an essential step on the path for volumetric quantification and visualization (Figure 6 and Video S1). Though computationally intensive, automated image segmentation has been widely reported in recent years for biomedical image segmentation (reviewed in [62]) as well as for mineral characterization [30–33]. However, up to now, the use of AI for retrieval of biologically interpretable results in the context of coral skeleton morphology has never been reported.

The impressive, robust performance of the DL-based segmentation in identifying and distinguish RADs and TDs across the coral skeleton particularly applies in the case of DL networks optimized for single datasets independently. In contrast, a network that was trained on one specific tomographic dataset yields severe mismatches when applied to segmentation of a different dataset, presumably because it is presented with previously unencountered shape configurations of the target skeletal features (Figures 7 and 8, Supplemental Figure S3). This suggests that the variability in coral geometries must be considered and possibly quantified, when generalizing findings, e.g., TDs or RADs morphologies. Furthermore, this suggests that the performance of AI in segmenting coral skeletal structures is quite sensitive to the large variability of the shape of RADs and TDs that exist among different recruits. Such shape variability is especially evident when comparing control pH recruits to low pH recruits, suggesting a morphologically-altering effect of OA on coral skeletal features that will require further in-depth examination to be deciphered. We speculate that the significant morphological differences tip the ratio between TDs and RADs in the skeletons in such a way that AI needs much more information to fully characterize the effects of OA. Note that the number of regions and slices presented to the DL for training in the normal pH and OA affected recruits was identical. This clearly indicates that normal pH corals are more similar to one another than to the pH affected corals.

Indeed, varying appearance of target features poses a great challenge in image segmentation (reviewed in [54]). For example in medical research, there can be wide differences between different medical images that affect the adaptability of DL network models during segmentation [62,63]. This is mainly due to the heterogeneous appearance of target organs for segmentation, which may vary hugely in size, shape and location from one patient to the other [64]. Another aspect that may restrict the performance of DL-based image segmentation is the limited diversity of data available to train a network [54,62,63]. Here, we show that widening the spectrum of possible appearances of the target skeletal features to be recognized using AI across multiple datasets considerably improves the DL-based segmentation performance (Figures 7 and 8). This showcases the power of AI when fed with diverse training data, and the potential of application across multiple tomographic datasets.

A commonly adopted method to generate new samples and increase the size of the training datasets is data augmentation, which is the application of a set of image transformations (e.g., rotation or flipping) to the dataset [65]. In this way, previously unencountered target features can potentially be approximated by the network. However, this method also adds sources of error and variability to the general methodological framework, as so far there is no strict consensus on which transformation to introduce and generally the transformation parameters are inconstantly chosen [54]. In this regard, we aim at improving the performance and generalizability of DL models across diverse coral datasets by collecting a larger diverse array of coral tomography data. If such datasets are made publicly available through, for example, open-access online databases, they will form a source of data that can be used to train robust, efficient, and generalizable DL models for use in coral research. It should be noted that image processing and segmentation of entire coral polyps with resolutions of ~2 μm and, in particular, using DL networks require very large computation memory and power, for which high-performance workstations may not be ubiquitously available. Yet, providing tools to help researchers analyze and utilize such data is essential, thus forming a data-collection-analysis feedback loop. Our work presented here is a first step on this ambitious path.

#### **5. Conclusions**

To the best of our knowledge, this is the first study that applied AI approaches to investigate the internal 3D geometry of RADs and TDs in coral skeletons based on highsensitivity, non-destructive PCE-CT data. The approach presented here opens the possibility of using AI to reconstruct and quantitatively analyze the internal skeletal network of reefbuilding corals. In particular, such technology has the potential to improve ecological monitoring programs, by assessing coral skeletal anatomy and providing morphological descriptors in order to taxonomize hard corals to genus and species levels. In addition, it can expand our knowledge on coral biomineralization by providing a means of characterizing the microscale interfaces within the different coral skeletal regions. Lastly, such a tool allows us to study the dynamics of coral skeletal formation under a variety of environmental conditions, including predicted OA scenarios.

In the last decades we witnessed a significant decline in coral reefs all over the world, and as global climate changes are becoming more severe, the need to protect coral reefs is becoming urgent. Hence, applying materials science contrast-enhanced methods to structural research in corals would be of great benefit for the development of highly targeted restoration strategies, as it would critically expand our knowledge on the effects of environmental change on coral skeleton development.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/jmse10030391/s1, Figure S1: Coral sample fixed in the sample holder in front of the X-ray source. Figure S2: Density variation among RADs and TDs. Figure S3: Skewness of the distributions of the per-slice RADs-TDs shape configuration across tomographic datasets. Table S1: Parameters of the carbonate chemistry across experimental pH conditions. Video S1: 3D view of coral skeletal structures before and after DL-based-segmentation.

**Author Contributions:** Conceptualization, F.S., T.M. and P.Z.; Formal analysis, F.S.; Funding acquisition, T.M. and P.Z.; Investigation, F.S., K.S. and P.Z.; Methodology, F.S. and P.Z.; Project administration, T.M. and P.Z.; Supervision, T.M. and P.Z.; Visualization, F.S.; Writing—original draft, F.S. and P.Z.; Writing—review & editing, F.S., K.S., T.M. and P.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from GIF, the German-Israeli Foundation for Scientific Research and Development (I-1496-302.8), from the Israeli Binational Science Foundation (BSF 2016321 to T.M.) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 755876 to T.M.). The experiment was performed in a controlled aquarium system which was funded by Institutional ISF grants 2288/16.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplemental Materials. Raw X-ray μCT data are available upon request from the corresponding authors.

**Acknowledgments:** We thank the HZB for granting us the use of the BAMline, specifically Henning Markötter and Michael Sintschuk for beamtime setup. We also thank Maayan Neder for her help in collecting coral larvae from the wild, Tal Zaquin for his assistance with setting up the OA experiment and Ole Lenz for his help in establishing the reconstruction pipeline and guidance in dataset processing.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

