The presence–absence data regarding the coralligenous formations in the Greek Seas were mined by the archive of the Laboratory of Marine Geology and Physical Oceanography (https://oceanus-lab.upatras.gr/), Department of Geology, University of Patras, Greece. The data were derived from 33 marine geophysical expeditions and cover areas all around the Greek Seas, spanning three decades of marine geo-acoustic surveying for mostly geology-oriented, commercial, or research projects, as listed in the Supplementary Document ( S1 ). Those data can be segmented into two major categories: analogue (printed records in thermal papers) that were collected between 1990 and 2007 and digital ones collected from 2007 and onwards. The leading marine geo-acoustic system used to acquire the analogue data was a 100-kHz EG&G 272 TD SSS, while a dual-frequency (100 kHz and 400 kHz) EdgeTech 4200 SSS was used in the more recent works. Thus, SSS was the main data source for seafloor characterization, although it was very often supplemented by 3.5 kHz or Chirp SBP data, ROV, or towed underwater camera (TUC) and sediment samplings for visual census and ground truthing. Although SSS and MBES backscatter data analysis are the most prominent techniques of marine remote sensing used for the detection of coralligenous formations (Roberts et al., 2005; Ballesteros, 2006; Collier and Humber, 2007; Cogan et al., 2009; Le Bas and Huvenne, 2009; Martin et al., 2014; Fakiris et al., 2019; Sañé et al., 2021), SSS and SBP multi-platform hydroacoustic investigation has proven to be an efficient method, as described thoroughly before (Georgiadis et al., 2009; Fakiris et al., 2018; Fakiris et al., 2019; Dimas et al., 2022).

In cases where historical seafloor mappings already existed in the archived survey reports, they have thoroughly been re-examined to apply corrections to any misinterpreted areas. Seafloor biogenic formations were once unexploited seabed features, often resembling geological components (i.e., rocky outcrops), and as such, they might have in some cases been misclassified. Un-processed archived analogue or digital data were suitably processed and interpreted to produce thematic maps in common geospatial vector data format for geographic information system (GIS) (*.shp polygons). While curation of the digital data was a straightforward process that was limited to collecting the data into single GIS data groups and applying, when necessary, geographic projection conversions into World Geodetic System (WGS) 84, processing of analogue seafloor maps involved digitization through paper scanning, map stitching, and georeferencing in an ArcMap environment. The final vector file of the archived data covered an area of 3,197.68 km², including adequate parts of the Ionian and South Aegean Seas and key areas of the North Aegean Sea, as presented in Figure 1 .

In Dimas et al. (2022), the guidelines for detecting coralligenous formations in SSS and SBP are provided in detail, including characteristic echo-signatures for various coralligenous formations, such as reefs, banks, maerl, or rhodolith beds, and the same criteria have been applied in the present study. Even since, Pérès and Picard (1951) observed that coralligenous banks and reefs are often surrounded by rhodolith or maerl sedimentary beds and expressed the opinion that they were developed from the coalescence of the latter ones. Early research (Got and Laubier, 1968; Laborel, 1987) further suggested that considerable in-height buildups grow almost always upon rocky outcrops. It is also supported by our knowledge (Georgiadis et al., 2009; Dimas et al., 2022) that transformation from superficial coralligenous formations, such as maerl or rhodoliths, to buildups, such as banks or reefs, is strongly connected to local geological settings, taking place in sub- or out-cropping bedrock areas, creating a more or less continuous transition from one type to the other. Figure 2 shows SSS record examples of coralligenous formations in the Cyclades Plateau, central Aegean Sea, with various transitions from the reef or “minute reef” (term by Georgiadis et al., 2009) to maerl or rhodolith bed types. It has been noted that, in many cases, high buildups are surrounded by a peripheral transition zone with maerl or rhodolithic composition and that sometimes these zones connect to form an extended maerl/rhodolith bed with scattered buildups within, forming various reef aggregation densities and transition zone intensities. Fakiris (2012) quantified those aggregations in selected areas in the Cyclades Plateau, depending on the ratio between reef density, transition zone intensity, and connectivity using region of interest (ROI) analysis (Fakiris and Papatheodorou, 2012), and suggested four main groups of aggregation types (see Figure 2 ). The first group refers to dense “minute” reef aggregations without well-developed maerl/rhodolithic transition zones and is found on the coastal slopes of Cyclades islands in relatively shallow depths (38–55 m). The second group refers to isolated reefs with high-intensity buffer zone each, located in the deepest limits of the coralligenous plateau (70–100 m) laying on a flat seabed. The third group is also present in deep and flat beds, but it is composed of dense and smaller “minute” reefs with buffer zones that connect to form a uniform maerl/rhodolithic bed. The fourth group is found in intermediate depth, slope, and proximity to the coasts and is characterized by high “minute” reef densities and un-connected buffer zones. Although most of our data resemble the echo-signatures of one of the above aggregation type groups, there were cases where coralligenous banks were the main detected formation, especially in the Ionian Sea, while in the North Aegean Sea, there were limited cases of maerl beds (high backscatter areas validated through visual census) without the presence of any buildup formation. Given the heterogeneity in the resolution and quality of the data used in this study (composed of both digital and analogue sonar records), the ambiguous transitions from one formation type to another, and the relative sparsity of samplings in the Greek Seas, it was decided that only two seabed types will be used in the analysis, namely, the “coralligenous” (presence) and the “non-coralligenous” (absence) ones ( Figure 1 ). In the case of a coralligenous aggregation, its outer boundary, including any reefs or superficial maerl or rhodolith areas that seem to be spatially delimited, was used as a uniform classification thematic vector layer in the GIS digitization process.

After evaluating and curating the archived mapping data, 6.67% (211.67 km²) of the mapped seafloor area was found to be covered by coralligenous formations, independently of which type (reefs, maerl, and rhodolith beds) they have. To provide the means for ecoregion-scale analysis, the Greek Seas were partitioned into three distinct regions, namely, the Ionian, S. Aegean, and N. Aegean Seas, following Panayotidis et al. (2022), as illustrated in Figure 1 .

Environmental variables were retrieved from open data repositories and model products. Bathymetry and light attenuation, the latter including both photosynthetically active radiation (PAR) on the seafloor and average diffuse attenuation PAR (KDPAR), were retrieved from the European Marine Observation and Data Network (EMODnet) bathymetry (https://www.emodnet-bathymetry.eu/) and seabed habitats (https://emodnet.ec.europa.eu/en/seabed-habitats) EU portals, respectively. Physicochemical attributes of the seawater were retrieved from the Copernicus Marine Services EU data portal, using the model products of the Sea Physics and the Biogeochemistry Reanalysis systems. More data regarding light attenuation and specifically Lee's euphotic depth (Lee et al., 2007), as estimated through satellite-based remote-sensing, were retrieved from NASA's OceanColor Web (https://oceancolor.gsfc.nasa.gov/) (NASA, 2022), part of NASA's Ocean Biology Processing Group (OBPG), supported within the framework and facilities of the NASA Ocean Data Processing System (ODPS).

Apart from EMODnet products that have been provided in ASCII (*.csv) raster format of 100-m pixel sizing, all other datasets were in NetCDF (*.nc) format with 4-km resolution. Copernicus Sea Physics and Biogeochemistry Analysis and Forecast Systems’ products regarded depth layers from the surface to 1,400-m water column depth and monthly averages of three consecutive years (April 2017–April 2020), finally accounting for 13 variables with 36 monthly time layers and 84 depth layers (3,024 total layers) each. Those were reprocessed through customized Matlab scripts to keep only depth cells closest to the seafloor and to estimate the root mean squared (RMS) value of each variable over the 3 years' monthly data. All Copernicus and OceanColor Web data were resampled to 100-m cell size (i.e., the resolution of EMODnet bathymetry and light indices) using the nearest neighbour interpolator in QGIS. An up-sampling of such an order was decided to honour, as much as possible, the high-resolution mapping archive and achieve more data points within each coralligenous zone. Examples, indicative of the mapping resolution used, are shown in Figures 1B, C , where it is shown that a 4-km sampling resolution would greatly decrease the spatial fidelity of our mapping dataset, allowing an insignificant number of samples with coralligenous formations to be included in the training dataset.

The retrieved bathymetry raster was further processed to extract a range of topographic indices. The Benthic Terrain Modeler 3.0 toolbox for ArcGIS (Walbridge et al., 2018) was used to derive slope, broad, and fine-scale bathymetric position index (BPI). BPI compares the elevation of a cell in a digital elevation model to the mean elevation of adjacent cells of a defined area. BPI datasets are produced through a neighbourhood analysis function. Positive BPI values illustrate the higher elevation of a specified cell compared to neighbouring cells, while negative values represent lower elevation. Values close to zero describe flat areas (Weiss, 2001). Thus, BPI is used to highlight depression and ridge areas. BPI scale factor is calculated by multiplying the outer circle radius (in cells) by the cell size resolution (Lundblad et al., 2006). For the present study, two BPI indices were extracted: a fine one and a broad one. Broad BPI (bBPI) concentric circles' radii were 25 and 250 with a scale factor was 25 km. For the fine BPI (fBPI), the radii were 5 and 25 with a scale factor of 2.5 km. Those scales have been chosen upon a trial-and-error process to capture both fine and large-scale bathymetric features owing to geology (out/sub-crops, basins, and faults) that potentially favours the development of coralligenous formations either directly, by offering suitable substrate, or indirectly, by controlling local ocean dynamics and water mixing processes. Caution should be taken for depth inaccuracies concerning the EMODnet bathymetry products, as deviations of up to ±20 m were identified in areas where the research team held bathymetric data.

The final list of the 21 environmental variables initially considered in this work is presented in Table 1 , along with brief descriptions, data sources, and source resolutions, while the predictors’ distribution maps are provided in Supplementary Data ( S2 ).

The thematic map containing the binary classes, i.e., “coralligenous” and “non-coralligenous” formations, areal information was converted to a 100 × 100 m cell size raster file, concentred with the environmental variable rasters, which was chosen as the “working” resolution. This raster was further converted to a point vector file (point *.shp) containing a total of 296,236 entries. For each point, the raster values of the 21 environmental variables were joined using the point sampling tool in QGIS along with their x–y coordinates in the metric UTM 34N geographic system. Separate fields were created regarding the binary classification scheme (“coralligenous” and “non-coralligenous” classes) and the relevant ecoregion (Ionian, S. Aegean, and N. Aegean Seas) for each data point. The environmental variable vectors, especially the modelled ones derived from the Copernicus portal, had important gaps in shallow or semi-enclosed areas of the Greek Seas due to their coarse horizontal and vertical resolutions. These gaps coincided with mapped areas in the archive, namely, the Amvrakikos and Evoikos Gulfs. Those gaps in spatial coverage were retained, so the above-mentioned Gulfs, although of high importance due to presenting coralligenous formations in shallow, turbid waters, were unable to be included in the predictive modelling and any other statistical analyses.

Modelling techniques are often sensitive to the multicollinearity of regarded predictor variables used. The 21 environmental variables were first tested for value distribution normality and then for linear correlation using Pearson’s correlation coefficient in PAST 4.03 statistical software (Hamer et al., 2001). Only PARbed was found to have a significantly lognormal distribution, so the logarithmic transformation was applied to it. The correlation matrix for the environmental variables (correlogram in Figure 3 ), clearly suggested 11 mutually un-correlated groups of variables. Each group was composed of one to six correlated variables (corr. Coeff > 0.7), and only one variable was selected as a representative for each group, finally forming a vector of 11 variables ( Table 2 ), from now on called the “predictors”. The data points, including the 11 predictors along with the assigned binary bottom classes, constituted the training dataset for the development of the spatial distribution predictive model.

The random forest (RF) classification algorithm was used in this work for predictive model development. RF has recently been used for habitat suitability modelling works (Pearman et al., 2020; Jackson-Bué et al., 2022) exhibiting significant results. It is preferred against other modelling techniques because it makes no underlying assumptions of the variables’ distributions, it is robust to overfitting, it allows for non-linear interactions between the response and environmental variables (Cutler et al., 2007; Zhang et al., 2019), and it can handle existence–absence data.

The Waikato Environment for Knowledge Analysis (WEKA V3.9.6) software, developed by the University of Waikato, Hamilton, New Zealand, was used for performing RF classification on the absence–presence data to build the spatial distribution probability model of the coralligenous formations in the Greek Seas. RF is an ensemble learning method for classification, which is based on multiple decision trees, thus the forest, for the final classification of a sample based on the majority vote. Through this procedure, the major drawback of overfitting the results based on training data by a single tree is diminished using multiple trees of random variables. The major component of its decision tree is its root, which is the start of the tree and contains the number of variables used (N) out of the total summary of them (K). The formula N = log2 K + 1 is mostly used for the selection of the number of variances N that each tree will adopt. Each step added to the decision tree leads through two branches into a leaf or to a new decision section called the internal branch. The summary of the steps needed until the final decision is the depth of the tree. The number of trees used for the creation of a forest is the iterations. One of the most important features of the RF is the Gini index, which ranks variable importance based on average impurity decrease. In other words, the Gini index describes the ability of each variable to divide the data present in the root or internal branch into two separate sets in order to advance into the next branch or leaf. All the available data have been used for training the RF model, while, through trial and error, the number of trees for each forest was set to 200, and each tree included any variable with a maximum depth of 10. The performance of the model was assessed employing both n-fold and spatial block cross-validation (Valavi et al., 2019) principles, as described in Section 2.6.

After the classifier was built and its validation metrics were estimated in the Weka explorer environment, the exported serialized classifier file was applied to the predictors’ vector file of the entire region of the Greek Seas. This vector file was created by converting the resampled 100 × 100 m cell size rasters of the predictors into point vector files and merging them under a common data table in QGIS. The final classification output contained the probability distribution for both “coralligenous” and “non-coralligenous” bottom classes, assigned to the predictors’ point vector file for the Greek Seas along with each point’s coordinates. Using a suitably selected distribution threshold, as discussed in Section 2.6, a final classification map, classifying points with probabilities of occurrence over the threshold as “coralligenous”, was created, which was suitable for further geostatistical compensation and geographic representation.

The output of the spatial predictive model is a raster of probability distribution values for coralligenous formations’ occurrence. This means that areas with higher distribution values have higher probabilities of being covered by coralligenous formations. As mentioned in Section 2.5, the RF model was trained on the full mapping dataset to take advantage of the wide geographic span of our data and to produce predictions that account for the heterogeneity among different geographic regions. Model performance was initially assessed through a 10-fold cross validation, where 10 different RF models were trained, leaving each time a random 10% of the data out for evaluation. Given that our model predicts the existence of coralligenous formations in vast areas, often far away from the training data areas (as shown in Figure 1 ), it was really important to take into account spatial autocorrelation. The “random” split of data in the 10-fold cross validation infers high spatial autocorrelation in the process, which might inflate the model performance. This is because of the high resolution of the mapping data, leading the random data points to be similar to each other. To test the spatial consistency of the generated decision rules, the model performance was also evaluated via spatial block cross validation (see Valavi et al., 2019). According to this, instead of generating random cross-validation folds, spatially separated folds (blocks) are created, and the RF model is trained each time onefold and tested on the others. In our case, the three ecoregions were used as “blocks”, and the ability of the model to correctly predict coralligenous occurrences when trained on one (Ionian Sea, N. Aegean, or S. Aegean Seas) or another was examined. The accuracy of the above models was evaluated and compared to each other using the area under the receiver operating characteristic (ROC) curve (AUC). The ROC curve is a graph showing the performance of a classification model, expressed as the ratio between the true-positive rate and false-positive rate, at all classification (probability distribution) thresholds. The AUC value is eventually a threshold invariant model evaluation metric that indicates how well the model fits the test set data, from 0 (model worse than random predictions) to 1 (ideal model). An AUC = 0.5 indicates random predictions, and an AUC < 0.5 indicates a model that is worse than when random classes are assigned.

When the probabilistic model needs to be reduced to a binary classification model, where there are areas “with” or “without” coralligenous formations, then a probability distribution threshold (classification threshold) needs to be set and assume that positive predictions occur when the probability distribution values are higher than this threshold. An optimized threshold would be able to respect the initial observations with a reasonable trade-off between true predictions and misclassifications. Specification of this threshold is not straightforward, but it can be decided upon observing the model accuracy change when using different threshold values. For each threshold value, an error matrix and a selection of metrics should be considered, which document the deviation between predicted and observed classes of the test set, giving an estimate of the model performance. A selection of standard performance metrics was created from the error matrix. No single measure can fully describe the performance of a classification model, and Precision, Recall (or sensitivity), and F-Measure are estimated. Precision answers the question: “What proportion of positive identifications was actually correct?” Recall: “What proportion of actual positives was identified correctly?” F-Measure is the harmonic mean of Precision and Recall, so it is a more robust measure of the model’s performance. When these measures are plotted against the classification threshold values, one can decide on the optimal threshold that keeps Precision and Recall balanced without sacrificing the overall model’s accuracy, as quantified through the F-Measure (Esposito et al., 2021). Selecting a high threshold would make Precision rise, as more and more predictions would be correct, but Recall would decrease as some true observations would be neglected. A good threshold is located where the F-Measure vs. threshold curve has been stabilized close to its maximum values until Recall starts becoming lower than Precision. The lowest possible threshold within the above range has been chosen so as to retain high positive prediction populations. Predictions with probability densities over this threshold were considered positive “coralligenous” classification, forming the final test set, for which evaluation metrics have been reported. Kappa, the true skill statistic (TSS), overall accuracy (OA), and specificity evaluation metrics have also been estimated, as suggested by Zhang et al. (2019).

Finally, it is essential to know how similar the conditions in the predictions are to those found in the training data, while in areas where they differ much, predictions should not be considered trustworthy. An out-of-range check, comparing the training to the full dataset variable numeric ranges, was made to indicate areas where at least one predictor is not within the ranges of the equivalent of the training set.

Statistical analyses in the context of the present study were performed for both the observations and the model predictions. Principal component analysis (PCA) with biplot interpretation was performed with the PAST statistical software to reveal any underlying mechanisms that link specific environmental/oceanographic variables to higher probabilities of “coralligenous” encounters. Violin plots were also created in the same software to compare the value distribution of the selected variables between the areas with or without coralligenous formations and to investigate any direct classification boundaries that each one offers. The latter was performed on the predicted “coralligenous” occurrences using the classification threshold decided through the model validation process.

We prepared one ROC curve for the 10-fold cross-validation process and another three for when training the model on each spatial block (i.e., ecoregion) and evaluating the others. Figure 4A presents the comparison between these ROC curves. All models exhibited high AUCs (>0.95), except for the one trained on the N. Aegean ecoregion block (AUC = 0.79), whose performance was apparently affected by the relatively low quantity of training data located therein. Nonetheless, the spatial autocorrelation of available data has been proven high, allowing the model to perform well when trained on either the Ionian or S. Aegean block, let alone when trained on the full dataset (exhibiting AUC = 0.99).

The reported model evaluation metrics in Table 3 correspond to the classification scheme that occurs when training the model on the full training dataset and validating it on the full dataset or either ecoregion after having selected the optimized classification threshold. As mentioned in Section 2.6, threshold selection was performed based on Precision, Recall, and F-Measure, plotted against the threshold, as illustrated in Figure 4B . This plot suggests an optimized threshold choice between 0.23 and 0.41. The minimum suggested threshold of 0.23 has eventually been selected to achieve the maximum population of positive predictions, retaining a good balance between Precision and Recall and a high F-Measure. For this threshold, all metrics indicate very good performance of the model, reaching overall accuracies between 95.89% and 98.81%, kappa values from 71.2% to 85.26%, and low false-positive rates (0.97%–3.34%), implying that the predictions meet the great majority of the observed classes, producing insignificant numbers of false predictions.

The spatial distribution model, showing the occurrence probabilities of coralligenous formations in the Greek Seas, is illustrated in Figure 5 , in which predictions accounting for the selected classification threshold (>0.23) are indicated with red colour, while any areas where at least one predictor is not within the ranges of the equivalent in the training set are shaded out as of low reliability. A citable GeoTiff raster image of the model output has been made available in the figshare (https://figshare.com/) open data repository (https://doi.org/10.6084/m9.figshare.22276456.v2).

The Gini index coefficient estimated within the RF algorithm ( Figure 6 ) indicated that Lee’s euphotic depth was the variable offering the highest information gain to the learning process, followed by CurrBott and PARbed with values ranging from 0.30 to 0.28. The topographic indices, depth, bBPI, fBPI, and slope exhibited intermediate index values between 0.26 and 0.24. Chlorophyll (0.21), pH (0.19), salinity (0.18), and alkalinity (0.17) offered the least information gain. The above implies that light radiation, bottom currents, and seabed morphometry seem to be the main factors controlling the spatial distribution of coralligenous formations.

In Figure 1 , the water depth probability density function for observations with coralligenous formations is presented per ecoregion (Ionian, S. Aegean, and N. Aegean Seas), clearly showing that they have not been detected deeper than 140 m. The same has been shown through spatial predictive modelling, as shown in the depth probability densities in Figure 5 . To closely investigate which environmental variables control the spatial distribution of coralligenous formations, the model predictions have thus been narrowed to the 0–140-m depth range, and PCA has been applied to the updated predictors’ vector. A PCA biplot ( Figure 7 ) has been created using the first two PCs, mixing PC scores with variables’ loadings in the same PC axes. Minimum volume ellipsoids enclosing 60% of the PC scores of each bottom class were also drawn to visualize the separation between “coralligenous” and “non-coralligenous” classes and to investigate which variables they are controlled by.

The environmental variables that have been found, through the above analyses, to be controlling the most the spatial occurrences of coralligenous formations were further analysed to detect any distinct data thresholds separating them from areas without such formations. For this, the probability density of each “important” predictor’s values is visually compared between areas with or without coralligenous formations through violin plots, as illustrated in Figure 8 .

In total, a mapped area of 3,197.68 km² has been archived, 6.62% (211.67 km²) of which was covered by coralligenous formations. In the 20–140 m depth range (1,051 km²), where the coralligenous formations are exclusively narrowed to, their coverage reaches 16.2% (169.94 km²). In the Ionian and S. Aegean Seas, the mapped depth distributions are similar, illustrating a mean depth of 90 m, while in N. Aegean, the mean depth is approximately 40 m (see Figure 1 ). However, the depth distribution of coralligenous formations, as estimated through the spatial probability modelling, indicates that there is a common distribution of depths in all ecoregions, with N. Aegean having a bimodal statistical distribution ( Figure 5 ). The shallower mode of this corresponds to the northern parts of N. Aegean, on the coasts of Chalkidiki and Kavala-Thassos, where the Black Sea water outflow causes significantly different physicochemical properties in the water column.

The area of Greece with the most extensively mapped areas of coralligenous formations was the S. Aegean Sea. In a total of 1,455.21 km² mapped seafloor, 7.82% (113.8 km²) had coralligenous formations. The Cyclades Plateau is where the most extensive and well-formed reefs of all Greek sites were found. Significant formations were also mapped in the Kythira Elafonisos Strait (southwestern edge of the south Aegean Sea) and around Leros Island (Dodecanese, S. Aegean Sea). In the Saronic Gulf, formations were mapped in its central part and its southern border. The N. Aegean Sea was the least mapped area of the Greek Seas by LMGPO (524.36 km² where coralligenous formations were found in 4.18% (21.93 km²) of the surveyed areas. Coralligenous formations were also mapped in the Evoikos, Pagasitikos, and Toroneos Gulfs and Kavala Bay. In the N. Evoikos Gulf, the most well-formed buildups were mapped. In the Ionian Sea, a total area of 906.85 km² was mapped, and coralligenous formations were found in 75.92 km² (8.37%). In the northern part of the Ionian Sea, most of the formations were found north of Corfu Island and towards the mainland. Major formations were also mapped in the S. Ionian Sea and outer Patraikos Gulf. Formations were also mapped in the Amvrakikos Gulf, although it is a semi-enclosed, fjord-like marine environment. In the Evoikos and Amvrakikos Gulfs, there was no data coverage regarding the environmental variables, and thus, no predictive mapping has been performed there.

Model predictions honoured in great precision the observations and indicated areas where coralligenous formations are likely to have significant extents ( Figure 5 ). In the Ionian Sea, the model projected three areas with high probabilities of coralligenous formations' existence, namely, N. Ionian (Corfu and Paxos Islands), central Ionian (outer Patraikos Gulf and Zakynthos–Kefalonia Straits), and SW Peloponnese (Kyparissiakos Gulf to Pylos). In the S. Aegean Sea, four areas exhibit high coralligenous formations’ occurrence probabilities, namely, the Cyclades Plateau, N. Crete, the Dodecanese Islands, and the Kythira–Crete strait. In the N. Aegean Sea, Kavala Bay, N. Aegean Islands (Limnos, Lesvos, and Chios), outer Thermaikos, N. Pagasitikos, and Ierissos Gulfs are likely to have extensive coralligenous formations.

The spatial patterns of coralligenous formations’ occurrence probabilities in Greece as estimated in the present work are generally similar to those in Martin et al. (2013). Some major deviations though can be detected in Lesvos and Chios (NE Aegean), between Ikaria and Kos (SE Aegean) and in S. Corfu (N. Ionian) islands, where more abundant formations are predicted in the present work, while in E Limnos (N. Aegean), no formations are predicted herein, as justified also by our field data. The validity of any of the above models could only be assessed with the acquisition of more mapping data via targeted samplings.

The learning process of the RF predictive model involved, through the Gini index, weighting the importance of each environmental predictor regarding the degree it aided in producing correct decisions. As shown in Figure 6 , light irradiance, expressed via the Lee euphotic depth and the seabed PAR, and hydrodynamics (bottom currents) exhibited top gains in the classifier’s learning process, followed by all the topographic indices. The PCA biplot ( Figure 7 ) explaining 57% of the data variance projected almost all predictors with high loadings in the PC1 vs. PC2 axes. Lee and salinity vectors were positively correlated to the PC scores of the coralligenous formations; alkalinity, pH, and chlorophyll vectors negatively correlated to them, while all topographic indices were orthogonal (perpendicular) to the aforementioned ones. The predictive power of the selected variables is thus proven significant, given that the predictions’ scores exhibited significant separation between “coralligenous” and “non-coralligenous” classes in the PC1 vs. PC2 space, especially when combining at least two variables orthogonal to each other PC loading, forming well-separated prediction classes in the Euclidean space.

Regarding light attenuation, the violin plots of Lee and PARbed in the 0–140 m depth range ( Figure 8 ) indicate that coralligenous formations are expected where the euphotic depth is high (greater absorption of light irradiance in the water column) and the PAR on the seafloor is moderate. This is attributed to the sciaphilic nature of the coralligenous algae. In this study, coralligenous formations have been found and predicted in areas where the irradiance reaching the seafloor (PARbed) was between 0.005% and 3% of the surface one, which is also supported by Ballesteros (2006). Moreover, the violin plots indicate that the majority of coralligenous formations are found where Lee is in a narrow range between 65 and 68 m and the depth is below −70 m. Lee euphotic depth is derived as the depth in the water column that PAR is 1% of its surface value; thus, the above ranges imply that formations are mostly found where PAR is even lower.

Bottom currents (CurrBott), although not exhibiting high loadings in the PC1 vs. PC2 biplot (it had greater loadings in PC4), nor any distinctive range thresholds can be drawn in the violin plots ( Figure 8 ), have been attributed to one of the highest importance scores due to the Gini index. The importance of hydrodynamics in coralligenous formation development has previously been reported in the literature. Unattached corals, like maerl and rhodoliths, are subject to mechanical abrasion under strong currents (Basso et al., 2009). Water movement, to some level, however, is crucial, as it enables oxygen circulation and keeps their surface clean (Larkum et al., 2003) while simultaneously protecting them from poisoning through water stagnation (Basso, 1998; Nelson et al., 2012). Moreover, as bottom currents are the means for sediment transport, they are controlling sedimentation rates on the seafloor. Since coralligenous formations have very slow growth rates (Ballesteros, 2006; Caragnano et al., 2016), areas where the sedimentation rate is higher than the growth rate are not suitable for their development (Basso, 1998; Dethier and Steneck, 2001; Foster, 2001; Riul et al., 2009; Villas-Bôas et al., 2014). Thus, coralligenous formations build in a specific range of bottom current speeds, as they require currents strong enough to mix the seawater but not that strong as to detach the thalli of the formations from the seafloor. More closely examining bottom current speed distribution in the violin plots, no formation has been in areas with less than 0.5 cm/s water speeds, while few develop in areas with over 3 cm/s of water speed.

Topographic variables have been given relatively high importance in the model development, as expressed through the Gini index, and they also exhibited high PC loadings in the PCA biplot. Especially slope, bBPI, and fBPI exhibit considerably different value distributions in areas with coralligenous formations than in areas without ones ( Figure 8 ), showcasing the importance of topographic indices for the suitability modelling of coralligenous formations in the Mediterranean. Buildups seem to develop in areas with exclusively positive bBPI values, i.e., on large-scale ridges and mounds, avoiding flat basins and valleys. The same is supported by fBPI, which exhibits 25th–75th percentiles in the −1 to 61 range for areas classified as “coralligenous” and −3 to 7 for “non-coralligenous”, implying that buildups developed on finer-scale bathymetric ridges and mounds and their lower slopes. Likewise, coralligenous formations seem to avoid slopes close to 0 and develop on areas with slopes in the range of 1° to 5° and mainly on mild slopes (mean of 1.5°). The above is in agreement that coralligenous formations (especially reefs) are known to develop in sub-cropping or out-cropping areas (Got and Laubier, 1968; Laborel, 1987; Georgiadis et al., 2009; Dimas et al., 2022), which inherently form ridges and mounds on the seafloor, including the coastal slopes where the bedrock is also sub- or out-cropping. This finding is an addition to the previous knowledge. For instance, in Georgiadis et al. (2009) and Martin et al. (2013), the Cyclades Plateau was considered to have high occurrence probabilities of coralligenous formations throughout its 70–100 m bathymetric range. This was due to the insufficient presence–absence data used for drawing assumptions or developing predictive models, making it difficult to draw the dependency between coralligenous formations and local seabed morphology. In the present study, the spatial extents of coralligenous formations in the Cyclades Plateau are predicted on peripheral to the Plateau large-scale ridges and upper slopes, exactly validating the plethora of observations that have been archived in the area.

Both salinity and alkalinity values in areas with coralligenous formations exhibited very distinctive distributions in the violin plots of Figure 8 . Their loading vectors were also significant in the PCA biplot ( Figure 7 ), implying a preference for coralligenous formations to develop in bottoms with high salinity and low alkalinity. More specifically, the “coralligenous” class shows a clear preference for the upper range of the modelled bottom salinity values, with its 25th–75th percentiles being 39.1–39.2 PSU out of a total variation in the 0–140 m depth range of the Greek Seas of 37.8–39.3 PSU. Alkalinity exhibited a clear threshold at 2.72 mmol/m³, with coralligenous formations appearing to occur in the narrow 2.71–2.72 mmol/m³ range, and without ones almost exclusively over it.

pH distribution, although obtaining middle-to-low values (more acidic), is in the narrow range of 8.106–8.116 for seafloors classified as “coralligenous”, in comparison to the 8.11–8.146 range for “non-coralligenous. The same occurs for chlorophyll, with coralligenous formations preferring moderate values, between 0.086 and 0.153 mg/m³ out of the 0.11–0.19 mg/m³ total range in the 0–140 m depth range. Even though chlorophyll and pH value distribution did not show any clear separatory thresholds between areas with and without coralligenous formations, they seem to have been more important than salinity and alkalinity in the classification process according to the Gini index ( Figure 6 ), while they also take high PC loading in the PCA biplot ( Figure 7 ). This is because the formations mapped on the Thracian coastal areas on the N end of N. Aegean are directly influenced by the exchanges with the high pH and chlorophyll Black Sea water masses, so the RF model ranked them high in an effort to correctly classify them in space. Finally, although the bottom temperature has not been selected as a model predictor (being linearly correlated to depth), it is an important environmental variable, and as such, its distribution is examined. It exhibits similar median values for areas with or without coralligenous formations (16.8°C), not forming any sharp threshold between them. Nevertheless, coralligenous formations show preference in the narrow range between 16°C and 17°C (distribution peak), while its main distribution in the Greek shallow bottoms (depth < 140 m) is between 15°C and 18°C.

There is a discussion about whether Mediterranean coralligenous formations in deep waters (>50 m) are still active or not. As discussed in Ballesteros (2006), the growth rate in formations found in the NW Mediterranean, at depths between 10 and 60 m, presented a great variety between 0.006 and 0.83 mm/year, depending on the depth they were found and the geologic period they were formed. Two growth rates were distinguished: a) higher growth rate (0.11 to 0.42 mm/year), placed in relatively shallow areas of 10 to 35 m, and b) low growth rates, reaching 0 in the 50–60 m depth. The latter corresponds to formations developed between 8k and 5k years BP, a geologic period when they have reached their maximum growth rate, estimated at approximately 0.20 to 0.83 mm/year.

Our field data show a maximum coverage of coralligenous formations between 63 and 110 m depth (mean depth of 90 m). If we assume a water depth of 35 m for maximum coralligenous algae growth (Ballesteros, 2006) and based on the global eustatic sea level curves (Lambeck et al., 2014), then those coralligenous formations should have been developed approximately 13 and 9 ka BP. This interval corresponds to a transitional period between the Last Glacial low sea-level stand and the Holocene high sea-level stand. Sedimentological studies have shown that in this period, the seafloor in the studied area started to experience more active bottom currents, together with the sea level rise (Tripsanas et al., 2015) and more eutrophicated conditions linked to increases in precipitation, high levels of marine productivity, and strengthening of water-column stratification (Rohling et al., 2015).

However, in the Thracian Sea (N end of N. Aegean), in Amvrakikos Gulf and N Evoikos Gulf, the shallowest coralligenous formations have been detected in depths of between 20 and 50 m. These areas are characterized by very turbid waters and shallow euphotic depths, as they are affected by river outflows and the Black Sea water mass exchanges, or they are semi-enclosed gulfs with high chlorophyll productions, likely forming suitable conditions for coralligenous algae (known to be sciaphilic). Whether these formations are active or not today and, if they are, which is the trigger mechanism for their growth, especially in the deeper waters, are high-priority research topics that have not yet been studied.