Predictive Mapping of Mediterranean Seagrasses-Exploring the Inﬂuence of Seaﬂoor Light and Wave Energy on Their Fine-Scale Spatial Variability

: Seagrasses are ﬂowering plants, adapted to marine environments, that are highly diverse in the Mediterranean Sea and provide a variety of ecosystem services. It is commonly recognized that light availability sets the lower limit of seagrass bathymetric distribution, while the upper limit depends on the level of bottom disturbance by currents and waves. In this work, detailed distribution of seagrass, obtained through geoacoustic habitat mapping and optical ground truthing, is correlated to wave energy and light on the seaﬂoor of the Marine Protected Area of Laganas Bay, Zakynthos Island, Greece, where the seagrasses Posidonia oceanica and Cymodocea nodosa form extensive meadows. Mean wave energy on the seaﬂoor was estimated through wave propagation modeling, while the photosynthetically active radiation through open-access satellite-derived light parameters, reduced to the seaﬂoor using the detailed acquired bathymetry. A signiﬁcant correlation of seagrass distribution with wave energy and light was made clear, allowing for performing ﬁne-scale predictive seagrass mapping using a random forest classiﬁer. The predicted distributions exhibited >80% overall accuracy for P. oceanica and >90% for C. nodosa , indicating that ﬁne-scale seagrass predictive mapping in the Mediterranean can be performed robustly through bottom wave energy and light, especially when detailed bathymetric data exist to allow for accurate estimations.


Introduction
Seagrasses are marine flowering plants that grow in shallow coastal waters and cover vast areas. The dense underwater meadows they create are some of the most biologically diverse and productive elements of coastal systems [1]. Since these meadows alter their abiotic environment (e.g., through build-up and modification of local hydrodynamic conditions), they are considered ecosystem engineers [2]. Five different species of seagrass can be found in the Mediterranean Sea [3]: Posidonia oceanica (endemic), Cymodocea nodosa, Zostera noltei, Zostera marina (glacial relict) [4], and Halophila stipulacea, an alien species introduced from the Red Sea [5]. In terms of distribution and abundance, P. oceanica is the most significant seagrass species out of the five [6]. According to [7], P. oceanica grows in depths between 0 m and 45 m. Approximately 1-2% of the 0-50 m depth zone is covered by this biomass, and it covers an estimated area of 2.5 to 5 million ha [8]. P. oceanica grows into extensive, densely covered meadows with 1 m tall leaves on average [9]. These meadows support a diverse community with many fauna species of economic interest and offer Remote Sens. 2023, 15, 2943 3 of 18 was shown that the near-bottom orbital velocity and P. oceanica coverage were inversely related, thus stating that the upper limit of seagrass distribution depends on the level of bottom disturbance by currents and waves. In 2010, it was found by [59] that wave energy is a major limiting factor for P. oceanica growth and survival at depths shallower than the wave-breaking depth. The landscape and spatial patterns of the seagrass habitat meadow were found to be altered by exposure to wind-wave energy in [60]. It was observed that, compared to locations exposed to low wave energy, high-energy environments support patchier P. oceanica meadow landscapes with less seagrass cover, more intricate patch shapes, and less complex intra-patch architecture. The differences in the plant's tolerance to wave forcing depending on the substratum type were identified by [61]. The colonization of rocky outcrops is crucial for P. oceanica's survival in areas with moderate to intense hydrodynamics as these substrates provide stability and greater attachment strength, which are essential elements for the successful establishment of seedlings [62,63]. C. nodosa, on the other hand, with its lower total leaf surface area compared to P. oceanica [64], experiences lower drag and so it can be met in areas with higher hydrodynamics [65].
The present work aims to provide a fine-scale case study correlating the spatial distribution of two Mediterranean seagrass species, namely P. oceanica and C. nodosa, with wave and light energies dissipated on the seafloor, in a marine protected area, i.e., Laganas Bay, Zakinthos Island, Greece, where extensive P. oceanica and C. nodosa meadows exist. Precise and fine-scale acoustic habitat mapping in the study area has previously been conducted through the utilization of multibeam echosounder, sidescan sonar, and towed underwater camera for ground truthing [31]. The derived bathymetry was employed to simulate the dominant wave propagation patterns in the given region, establishing correlations between the respective orbital velocities and backwash current speeds with seagrass distribution. Photosynthetically active radiation (PAR) on the seafloor was estimated via open access satellite sea-surface derivatives, reduced to the seafloor by taking the light diffuse attenuation coefficient into account. Clear correlations between wave energy and light at the seabed with both P. oceanica and C. nodosa distributions were found, which further allowed the validation of a habitat predictive model using depth, wave, and light as the only predictors. A random forest predictive model was developed, realistically trained over areas with visual census ground truthing or over random data splits, exhibiting predicted habitat distributions with overall accuracies of more than 80% and up to 92%, implying worthy use in un-mapped coastal areas of the Mediterranean Sea.

Study Area
The study area includes the seafloor in the marine protected area of Laganas Bay, Zakynthos Island, Greece. The National Marine Park of Zakynthos (NMPZ) (Figure 1) is the first marine protected area (MPA) in the Mediterranean for the protection of sea turtles, including no-take zones [66]. It was established in 1999 with the goal to protect and manage one of the most important nesting rookeries of loggerhead sea turtles, Caretta caretta, in the region [67]. The NMPZ mainly contains four habitat types that are all included in the EU Habitats Directive 92/43/EEC, namely: Posidonia oceanica (EU habitat code-EUhc: 1120-"Posidonia beds"), rocky seafloor (EUhc: 1170-"Reefs"), sandbanks (EUhc 1110), and two types of soft substrates, i.e., "Unvegetated soft bottoms" (EUhc 119A) and "Vegetated soft bottoms" (119B), with the latter corresponding to Cymodocea nodosa beds. The first marine habitat mapping effort in the NMPZ was carried out by the Hellenic Center for Marine Research (HCMR) in 1999, in the context of the Annex I of the EU Habitats Directive 92/43/EEC (http://gis.ices.dk/geonetwork/srv/eng/catalog.search#/metadata/ce013 464-e509-40fc-b2dd-859e1e6daa0f, accessed on 15 May 2023). This was derived through aerial photos of scale 1:5000, data from the RoxAnn echosounder, phytobenthos samples and in situ observations by scuba diving, and is still used as the official marine habitat map for the NMPZ. In 2004, Pasqualini et al. [68], constructed a 2.5 m pixel-size marine habitat map of the NMPZ through analyzing SPOT 5 satellite multispectral imagery, including sandy and rocky seafloors and two types of P. oceanica, namely "patchy" and "continuous" ones. In 2019, the Oceanus-Lab (https://oceanus-lab.upatras.gr/, accessed on 15 May 2023) constructed the most detailed habitat map existing of the NMPZ, accompanied by a detailed bathymetry of the area, both derived via marine geophysical techniques described in [31] and in Paragraph 2.2. The latter constitutes the baseline source of information used in the present work, representing the most up-to-date and complete mapping effort in the area. echosounder, phytobenthos samples and in situ observations by scuba diving, and is still used as the official marine habitat map for the NMPZ. In 2004, Pasqualini et al. [68], constructed a 2.5 m pixel-size marine habitat map of the NMPZ through analyzing SPOT 5 satellite multispectral imagery, including sandy and rocky seafloors and two types of P. oceanica, namely "patchy" and "continuous" ones. In 2019, the Oceanus-Lab (https://oceanus-lab.upatras.gr/, accessed on 15 May2023) constructed the most detailed habitat map existing of the NMPZ, accompanied by a detailed bathymetry of the area, both derived via marine geophysical techniques described in [31] and in Paragraph 2.2. The latter constitutes the baseline source of information used in the present work, representing the most up-to-date and complete mapping effort in the area.

Background Knowledge: Laganas Bay Bathymetry and Acoustic Habitat Mapping
Detailed bathymetry and seafloor habitat mapping in the NMPZ were derived via marine geoacoustical means, including MBES and SSS, respectively. The mapped seafloor components were ground truthed and validated by conducting a thorough supplemental visual seafloor inspection survey with a towed underwater camera (TUC). A dual-head MBES Elac Nautic Seabeam 1185, transmitting at 180 kHz, and an Edgetech 4200 SP dualfrequency SSS, transmitting simultaneously at 100 and 400 kHz, were used to survey a total area of 84 km 2 . Data acquisition, processing, and analysis methods are thoroughly described in [31], where the first high-detail habitat map of the NMPZ was released (Figure 1). The depth of the bay was found to range from 0 to 145 m ( Figure 1) while thorough interpretation of the geoacoustical and ground truth data revealed that "Posidonia beds" cover 27% (22.5 km 2 ) of the seafloor, C. nodosa beds ("Vegetated soft bottoms") make up 6% (4.8 km 2 ), "Reefs" 9% (7.8 km 2 ) and soft sediments account for 58% (48 km 2 ) of the seafloor in the NMPZ. The P. oceanica habitat type was further divided into two distinct meadow sub-types, following [68], namely: (a) the "continuous" one, characterized as homogenous, with low canopy height and developed on plain substrate, and (b) the "patchy" one, corresponding to meadows with higher canopy height, developed on wellformed "matte" (rhizomes' build-up).

Background Knowledge: Laganas Bay Bathymetry and Acoustic Habitat Mapping
Detailed bathymetry and seafloor habitat mapping in the NMPZ were derived via marine geoacoustical means, including MBES and SSS, respectively. The mapped seafloor components were ground truthed and validated by conducting a thorough supplemental visual seafloor inspection survey with a towed underwater camera (TUC). A dual-head MBES Elac Nautic Seabeam 1185, transmitting at 180 kHz, and an Edgetech 4200 SP dualfrequency SSS, transmitting simultaneously at 100 and 400 kHz, were used to survey a total area of 84 km 2 . Data acquisition, processing, and analysis methods are thoroughly described in [31], where the first high-detail habitat map of the NMPZ was released (Figure 1). The depth of the bay was found to range from 0 to 145 m ( Figure 1) while thorough interpretation of the geoacoustical and ground truth data revealed that "Posidonia beds" cover 27% (22.5 km 2 ) of the seafloor, C. nodosa beds ("Vegetated soft bottoms") make up 6% (4.8 km 2 ), "Reefs" 9% (7.8 km 2 ) and soft sediments account for 58% (48 km 2 ) of the seafloor in the NMPZ. The P. oceanica habitat type was further divided into two distinct meadow sub-types, following [68], namely: (a) the "continuous" one, characterized as homogenous, with low canopy height and developed on plain substrate, and (b) the "patchy" one, corresponding to meadows with higher canopy height, developed on well-formed "matte" (rhizomes' build-up).

Numerical Modeling of Wave Propagation
The effect of wind-and wave-generated currents on the seabed hydrodynamics was investigated numerically utilizing the MIKE21 software developed by the Danish Hydraulic Institute (DHI). In particular, the following modules were used: MIKE21 Spectral Waves (SW) for the simulation of wave propagation, and MIKE21 Flow Model (FM) for the simulation of wave-and wind-generated currents.
The MIKE21 SW module [69,70] has the ability to simulate the wind-induced wave growth and transformation, based on the numerical solution of the wave action conservation equation. Depth-induced wave breaking is included in the SW module to model wave-energy dissipation in shallow waters, following the formulation in [71]. In the present application, the values of the model parameters are Alpha = 1 and Gamma = 0.8. The spatial discretization is performed using a cell-centered finite volume method in an unstructured computational mesh, and the temporal discretization is based on a fractional time step approach. On all open boundaries of the computational domain, the parameters, i.e., the Significant Wave Height (H s ) and the Peak Wave Period (T p ), of the incident deep-water waves are imposed. Due to lack of existing wave measurements, a hindcast method was used to obtain the wave parameters, H s and T p . Mean annual wind data of the period 1990-2021 were used, which were obtained from the meteorological station of the Hellenic National Meteorological Service at Zakynthos. Based on the mean annual wind data, the dominant wind speed and direction with a return period of 1 year were found to be 14.86 m/s (U 10 ) and 180 • (South), respectively. Then the fetch length for this wind direction was estimated and the corresponding deep-water wave parameters were calculated to be H s = 2.9 m and T p = 6.54 s using the SMB-Wilson methodology [72].
The output of the SW module is the distribution of the radiation stresses, which is used as input to the FM module, and acts as the driving forcing for the wave-induced currents. The MIKE21 FM [69] module is based on the numerical solution of the two-dimensional shallow water equations, i.e., the depth-averaged continuity equation over depth and the depth-averaged horizontal momentum equations. A finite volume method is used for the spatial discretization of the governing equations, while a second order Runge-Kutta scheme is utilized for the temporal discretization [72]. The land boundary (coastline) is considered an impermeable boundary with zero velocity, while free outflow boundary conditions are applied on all open boundaries of the computational domain. An eddy viscosity formulation is followed for turbulence modeling, utilizing the Smagorinsky (1963) model for the horizontal diffusion terms. In the present application, the value of the model constant is C s = 0.28. The bed shear stress is calculated using a quadratic friction law, where the drag coefficient is determined from the Manning coefficient, and is considered in the present application as equal to 0.03125 s·m −1/3 .
The bathymetric data for the wave propagation modeling were comprised of both the aforementioned in-situ acquired bathymetric measurements and data from the digital database DHI C-MAP, which includes measurements of the Hellenic Navy's Hydrographic Service (Figure 2a The wave-generated bottom orbital velocity (Ubr) (Figure 3c) was estimated from the surface-wave parameters Hs and Tp and the MBES bathymetry, following the formula described in [73]. Ubr is a key wave energy parameter on the seabed, proven to be well related to seagrasses in the works of [1] and [59]. Both Ubr and Uwi (in m/s) were exported to 30 × 30 m raster files.   The duration of the numerical simulations was set to 8 h, which corresponds to a typical duration for intense weather phenomena. A steady current velocity field was reached after approximately 2 h of simulation. The simulated wave-induced currents (Uwi) in Laganas Bay for the above-mentioned wave parameters are presented in Figure 3b. The wave-generated bottom orbital velocity (Ubr) (Figure 3c) was estimated from th surface-wave parameters Hs and Tp and the MBES bathymetry, following the formula de scribed in [73]. Ubr is a key wave energy parameter on the seabed, proven to be well re lated to seagrasses in the works of [1] and [59]. Both Ubr and Uwi (in m/s) were exporte to 30 × 30 m raster files.   The wave-generated bottom orbital velocity (Ubr) (Figure 3c) was estimated from the surface-wave parameters H s and T p and the MBES bathymetry, following the formula described in [73]. Ubr is a key wave energy parameter on the seabed, proven to be well related to seagrasses in the works of [1] and [59]. Both Ubr and Uwi (in m/s) were exported to 30 × 30 m raster files.

Light at Seabed
The quantity of light on the seabed of the NMPZ was estimated using the photosynthetically active radiation (PAR) (i.e., light at 490 nm) at the sea surface (PARs) reduced at depth z by the diffuse attenuation coefficient for PAR (KdPAR), following Beer-Lambert Law [74], where PARz is the photosynthetically active radiation on the seafloor. While more complex light attenuation models are available, the simple Beer-Lambert model, assuming exponential attenuation of light with depth, has been proven adequate for shallow waters [75] and has also recently been used for assessing the mesophotic zone in the Mediterranean Sea [76].
Although the directionality of the light field affects the diffuse attenuation coefficient, Kd-PAR, this influence has been found to be weak at mesophotic depths [74,77] [78], fully covering the NMPZ (Figure 3d). EMODnet also offers an estimate of PARz for all European waters, but it is calculated based on the EMODnet Bathymetry. The latter is of very low accuracy in most parts of the Greek coastal waters and so the herein-generated PARz raster, based on the MBES bathymetry of the NMPZ, offers the best feasible accuracy.

Statistical Analysis and Variable Importance
Geospatial processing was performed in ArcGIS 10.8. The raster files of the environmental variables were upscaled to meet the Depth raster, i.e., to 30 × 30 m cell sizes, using a nearest neighbor interpolator. The polygon shapefiles regarding seafloor habitat types were converted to a point vector file (point shapefile), with points respecting the cell centers of the raster files. For each point, the raster values of the four environmental variables (Depth, Ubr, Uwi and PARz) were joined using the Spatial Join (Analysis) tool along with their x-y coordinates in metric UTM 34N geographic system. The resulting variable database was the basis for any statistical and predictive modeling tasks.
Statistical analysis and representation of the variable database were performed in GraphPad Prism 9, IBM SPSS, and PAST. Principal Component Analysis (PCA) with biplot interpretation was performed in the PAST statistical software to reveal any underlying mechanisms that link environmental variables to the presence of each seagrass type. Violin plots were created in GraphPad Prism 9 to compare the value distribution of the environmental variables between areas with or without each type of seagrass cover and to investigate any direct classification boundaries that each one offers. Further statistical analysis took place employing correspondence analysis (CA), performed in IBM SPSS, to investigate the correlation between categorical variables, namely seagrass types and the four environmental variables (Depth, Ubr, Uwi and PARz). To utilize CA, variables were converted to ordinal, with depth categorized into 5 m interval classes, while Uwi, Ubr and PARz were converted into 0.005 m/s, 0.1 m/s and 2 mol·photon/m 2 /day classes, respectively.
An information gain attribute evaluator was used to evaluate the worth of each environmental variable in separating the probability distribution of the different seafloor types by measuring the information gain (entropy) with respect to any bottom class. Entropy quantifies how much information exists in the probability distribution of a random variable: a skewed distribution having low entropy, whereas a distribution where events have equal probability has larger entropy [79]. Entry values vary from 0 (no information contribution) to 1 (maximum information contribution). Variables with a higher information gain value can be ranked higher, whereas those with a lower score can be removed.

Predictive Modeling
The random forests (RF) classification algorithm was used in this work for predictive model development. The considered seafloor classes were reduced to "P. oceanica" (aggregating "patchy" and "continuous" meadow types) and "C. nodosa" beds, as well as a collective class named "Other", which included "Reefs", "Nearshore sandy seafloor" and "Offshore soft sediments". RF has recently been used for habitat suitability modeling works [80][81][82] exhibiting significant results. This modeling technique is favored over others because it makes no underlying assumptions of the variables' distributions, it is robust to overfitting, it allows for nonlinear interactions between the response and environmental variables [83,84], and it can handle existence-absence data, in contrast to the most important habitat suitability modeling software candidate, Maxent (e.g., in [85]), which can handle presence-only 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 seagrasses in Laganas Bay. 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 majority vote. The major component of its decision tree is their 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 internal branch. The summary of the steps needed until the final decision is the depth of the tree. In our case, the number of trees for each forest was set to 100 and each tree included 2 variables, with a maximum depth of 10.
Two distinct approaches were employed for the purpose of training the model. In the first case, database points coinciding with the TUC transect points (manually classified into seafloor habitat classes) were used for model training (training set) and the rest points were used for cross validation (test set). In the second case, 20% of the observations (20% split) were used as a training set and the remaining 80% as a test set. Training the classifier on the TUC transects can be considered a realistic case, according to which one needs a detailed bathymetry and a visual census ground truthing dataset to proceed to habitat predictive modeling without any prior knowledge of the spatial distribution of the bottom habitats. The 20% split case was decided to test the optimal accuracy of a well-trained classifier for building generalized decision boundaries.
After the classifier was trained on each training set and their validation metrics were estimated in the Weka explorer environment, the exported serialized classifiers were applied to the predictors' vector file (environmental variables' points) of the entire region of the NMPZ. The final classification output contained the probability distribution (in 0-1 range) for "P. oceanica", "C. nodosa", and "Other" bottom classes, assigned to the predictors' point vector file along with each point's coordinates. Using a 0.5 distribution threshold, a final classification map was created for each training case, which, having been converted back to raster files of 30 × 30 m resolution, was suitable for further geostatistical compensation and geographic representation. Some elementary statistics on the above variables, regarding either the overall coverage of P. oceanica and C. nodosa seagrass species, or their shallow and deep limits, are given in Table 1. The importance of each variable, as assessed using the information gain attribute evaluator, also given in Table 1, ranks Ubr (0.42) first and then Depth (0.40), followed by PARz (0.32) and, clearly last, Uwi (0.19). The PCA biplot, as applied to the environmental variables, is presented in Figure 4.  The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub- attribute evaluator, also given in Table 1, ranks Ubr (0.42) first and then Depth (0.40), followed by PARz (0.32) and, clearly last, Uwi (0.19). The PCA biplot, as applied to the environmental variables, is presented in Figure 4. The overall variance explained by PC 1 and PC 2 was 72% and 14%, respectively, with PC 1 demonstrating close to +0.8 loadings for all four parameters, while PC 2 demonstrated approximately +0.4 for Ubr and Uwi and −0.4 for Depth and PARz, indicating a negative correlation between those two groups. The four variables are grouped in 2 perpendicular principal axes, according to the biplot, one parallel to the maximum variances of the wave energy variables (Ubr and Uwi) and the other to the ones of PARz and Depth. Seagrasses were well separated from the rest of the seafloor in the PCA biplot, mostly found on negative PC 1 values (−2-1 range) and around zero PC 2 ones (−1.5-0.75 range), but the individual seagrass species data point clouds were mostly not well separated. The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub- The PCA biplot, as applied to the environmental variables, is presented in Figure 4.  The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub- The PCA biplot, as applied to the environmental variables, is presented in Figure 4.  The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub- The PCA biplot, as applied to the environmental variables, is presented in Figure 4. attribute evaluator, also given in Table 1, ranks Ubr (0.42) first and then Depth (0.40), followed by PARz (0.32) and, clearly last, Uwi (0.19). The PCA biplot, as applied to the environmental variables, is presented in Figure 4. The overall variance explained by PC 1 and PC 2 was 72% and 14%, respectively, with PC 1 demonstrating close to +0.8 loadings for all four parameters, while PC 2 demonstrated approximately +0.4 for Ubr and Uwi and −0.4 for Depth and PARz, indicating a negative correlation between those two groups. The four variables are grouped in 2 perpendicular principal axes, according to the biplot, one parallel to the maximum variances of the wave energy variables (Ubr and Uwi) and the other to the ones of PARz and Depth. Seagrasses were well separated from the rest of the seafloor in the PCA biplot, mostly found on negative PC 1 values (−2-1 range) and around zero PC 2 ones (−1.5-0.75 range), but the individual seagrass species data point clouds were mostly not well separated. The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub- The detailed value distributions of the four environmental variables are given in the form of violin plots in Figure 5, regarding all considered bottom types, namely: P. oceanica, C. nodosa and non-seagrass ("other"), while P. oceanica is further analyzed in its two sub-types, namely "patchy" and "continuous". Although the ranges of values and the simple statistical descriptives (i.e., mean and std) of most variables are overlapping between different bottom types, their exact distributions have characteristic signatures that seem suitable for separating one from another. To better investigate the relationship between the environmental variables and each seagrass type, their values were converted to ordinal (using distinct intervals) and correspondence analysis (CA) was applied to each variable. The two-dimensional CA plots of Figure 6 suggest distinct environmental variable ranges as being the most suitable for each seagrass type, especially regarding Ubr and PARz. emote Sens. 2023, 15, x FOR PEER REVIEW types, namely "patchy" and "continuous". Although the ranges of values a statistical descriptives (i.e., mean and std) of most variables are overlapping ferent bottom types, their exact distributions have characteristic signatures t able for separating one from another. To better investigate the relationship environmental variables and each seagrass type, their values were conver (using distinct intervals) and correspondence analysis (CA) was applied to The two-dimensional CA plots of Figure 6 suggest distinct environmental v as being the most suitable for each seagrass type, especially regarding Ubr a Figure 5. Violin plots of Depth (a), Uwi (b), Ubr (c) and PARz (d), for seagrass (P nodosa) and non-seagrass (other) seafloors. P. oceanica is further analyzed in its namely "patchy" and "continuous". , Ubr (c) and PARz (d), for seagrass (P. oceanica or C. nodosa) and non-seagrass (other) seafloors. P. oceanica is further analyzed in its two sub-types, namely "patchy" and "continuous". Figure 7 compares the original to the predicted classification maps as well as to the spatial probability distributions of P. oceanica and C. nodosa in the NMPZ, as generated by the random forest classifier trained either by the 20% split subset of the data or by datapoints along the TUC. The model validation metrics (Figure 7), sorted by seagrass species and training set, exhibited remarkable performance for either training set or seagrass type, with C. nodosa having been predicted up to 10% better by most metrics. The overall accuracy for C. nodosa was between 92% and 92.6%, while for P. oceanica the accuracy was between 83.4% and 84.6%. The F-measure for the realistic case of training the model along the TUC was almost equal between P. oceanica and C. nodosa types, reaching 77.5% and 76.4, respectively. , Ubr (c) and PARz (d), for seagrass (P. oceanica or C. nodosa) and non-seagrass (other) seafloors. P. oceanica is further analyzed in its two sub-types, namely "patchy" and "continuous".  Figure 7 compares the original to the predicted classification maps as well as to the spatial probability distributions of P. oceanica and C. nodosa in the NMPZ, as generated by the random forest classifier trained either by the 20% split subset of the data or by datapoints along the TUC. The model validation metrics (Figure 7), sorted by seagrass species and training set, exhibited remarkable performance for either training set or seagrass type, with C. nodosa having been predicted up to 10% better by most metrics. The overall accuracy for C. nodosa was between 92% and 92.6%, while for P. oceanica the accuracy was between 83.4% and 84.6%. The F-measure for the realistic case of training the model along the TUC was almost equal between P. oceanica and C. nodosa types, reaching 77.5% and 76.4, respectively.

Discussion
Although no clear threshold values separating P. oceanica from C. nodosa could be drawn from the elementary statistics of the four parameters (Table 1), correspondence analysis drew distinct favorable ranges between the seagrass species, especially in the Ubr and PARz domains. P. oceanica seems to favor the low Ubr (0.1-0.2 m/s) and middle PARz

Discussion
Although no clear threshold values separating P. oceanica from C. nodosa could be drawn from the elementary statistics of the four parameters (Table 1), correspondence analysis drew distinct favorable ranges between the seagrass species, especially in the Ubr and PARz domains. P. oceanica seems to favor the low Ubr (0.1-0.2 m/s) and middle PARz value ranges (6-14 mol·photons/m 2 /day), while C. nodosa prefers middle Ubr values (0.25-0.3 m/s) and sharp PARz values (around 10 mol·photons/m 2 /day).
PCA hardly managed to separate the two seagrass species in the principal components space and specifically in the PC 1 versus PC 2 one (Figure 4). The P. oceanica cluster is located, more or less, along the PC 1 axis in a range of −2 to +1 with two dense sub-clusters around −1.5 and 0-1.0. C. nodosa cluster seems to be clear mainly in the negative PC 1 and PC 2 values, located between the two P. oceanica sub-clusters and forming an elongated cluster down to −3.0 in the PC 2 axis. It is noteworthy that the PCA biplot exhibited a well-defined separation of both seagrasses from the rest of the seafloor, in the mid-low wave energy (Ubr and Uwi) principal axes versus the mid-low depth and light ones, respectively, summarizing the findings of the variables' distributions, as revealed in the violin plots of Figure 5.
Closer analysis on the variables' distribution per seagrass type, as expressed by the violin plots ( Figure 5), reveal that C. nodosa favors a narrow range of depths (around 15-20 m) while P. oceanica appears to flourish in a wider range of depths (5-40 m). Patchy and continuous P. oceanica meadows favor shallower and deeper depths, respectively. A reverse distribution skewness between P. oceanica and C. nodosa for both wave derivatives (Ubr and Uwi) was also seen, with C. nodosa favoring higher wave energies and P. oceanica non-zero, middle to low ones. C. nodosa appears to flourish in areas where the Ubr fall within a narrow range of 0.35-0.5 m/s. P. oceanica's median Ubr is 0.3 m/s, and the "continuous" sub-type, which represents relatively uniform meadows, is mostly found in areas where the Ubr is below 0.3 m/s. Uwi didn't show such a sharp transition between seagrass types. Nonetheless, the data indicates a distinct preference of seagrasses towards Uwi lower than 0.3-0.4 m/s. This pattern can be explained by the fact that when the Uwi exceeds this threshold it causes sand transport, which in turn forms seasonal beach equilibrium profiles. However, it appears that the uniform sub-type of continuous P. oceanica has a strong preference for lower Uwi, specifically within the narrow range of 0.2 to 0.3 m/s. The above observations are further supported by the spatial distribution of seagrasses and the other seabed types. The non-seagrass ("other") bottom types showed distribution peaks in both very low and very high Ubr values, below and above the seagrasses' equivalents, supporting the assumption that seafloor wave energy is suitable for seagrass existence when it is neither too low nor too high. This is likely because sediment deposition may occur in low energy environments, especially on sandy floors, and sediments resuspension and/or transport may occur in high energy ones, not allowing for seagrasses to root and expand on mobile sandy substrates. For instance, the nearshore sandy seafloor of the NMPZ (see Figure 1) is a sandbank of changing thickness, following the equilibrium beach profile, as controlled by the local swash-backwash ratio and longshore drift dynamics. In this region, as indicated in Figure 3b,c, high wave dynamics correlate very well with the lack of seagrasses on the seafloor from 0 to 25 m depth. The correlation of high Ubr values with the lack of seagrasses is also exhibited in the area delineated in Figure 3c(III), where P. oceanica extents follow exactly the spatial distribution of moderate Ubr values.
In the N sides of the two islets in the NMPZ, where wave dynamics are minimized due to the isles' wave-masking, sediment depositional conditions form wide beaches whose terrestrial and underwater profiles change vastly over the year in a seasonal fashion. In these wave-protected areas, Ubr is close to zero and its spatial distribution perfectly matches the absence of P. oceanica (see Figure 3c(I,II)). In areas where there are strong bottom currents, the presence of P. oceanica is noticeably absent. Those wave-induced currents, caused by back-wash action, have been simulated by the wave propagation model in the area indicated in Figure 3b(IV) in the SE part of NMPZ. There, a well-shaped downslope current is formed, developing a wedge-shaped sandy area with total absence of P. oceanica. In that area, sand-waves have been witnessed through visual census (see photo next to Figure 3b(IV)), which are likely a local recurring bottom feature, developing after periods with relatively high wave energies. The instability of the substrate observed in this area, as in the case of the nearshore sandy seafloor (sandbanks), pose the risk of burial and uprooting of plants [61] not allowing for seagrasses to flourish.
P. oceanica is found in a broader light radiation range (6-14 mol·photons/m 2 /day) compared to C. nodosa, which is restricted to a narrower range around 10 mol·photons/m 2 /day. Those findings, although indicative of the Laganas Bay case, should not be generalized because seagrasses might be competing with each other, especially where they seem to share a similar border. In central Laganas Bay, C. nodosa shares its deep limits with the shallow limits of P. oceanica and if one of them was not present, the other might have had wider extents. The deep limits of P. oceanica though, seem to well respect light radiation, having a mean ± std value of 4.5 ± 1.5 mol·photons/m 2 /day, much lower than that of C. nodosa which is 8.1 ± 1 mol·photons/m 2 /day.
What we are also convinced of, is that the shallow limits of P. oceanica are well controlled by Ubr, judging by the NW and NE areas of Laganas Bay, where northward wave energy is halted by the islets and P. oceanica flourishes in much shallower depths (~3-5 m), but starting from similar Ubr values (~0.1 m/s) as in the central Bay, where the seafloor is exposed to the open-sea waves. There, the same Ubr limit is found at approximately 25 m depth, where the presence of P. oceanica starts sharply. Further supporting this hypothesis, C. nodosa seems to be absent in the NW and NE areas of Laganas Bay, where its favored Ubr ranges don't seem to occur at all.
The suitability of wave and light dynamics on the seafloor was further verified in this study by their successful incorporation into the predictive modeling processes. Ubr, PARz and Depth alone, as implied by the attribute evaluation process, were evidenced to be the most significant predictors for such modeling. Random forest performed remarkably well for either seagrass species, with overall accuracies exceeding 80%, as suggested by the validation metrics (Figure 7), when trained on either ground truthing sites or on samples spread all over the NMPZ (20% data split). The above implies that habitat predictive modeling can be an effective tool for mapping Mediterranean seagrasses when seafloor wave and light energies are measured or modelled in fine scales. A key point for such high-precision modeling exercises is the existence of high-detail bathymetric data so that wave and light can be precisely modelled on the seafloor.

Conclusions
The research findings demonstrate that the predicted levels of wave energy and light on the seafloor are critical in determining the fine scale spatial distribution patterns of the most important Mediterranean seagrasses, namely C. nodosa and P. oceanica. Specifically, Ubr exhibited a gradual and well-expressed increase as the seafloor moves from nearshore sand towards areas inhabited by C. nodosa and finally P. oceanica. This transition has a clear spatial match with Ubr decrease, as evidenced even when habitat boundaries are examined at finer scales. The deep limits of P. oceanica seem to be well controlled by light irradiance, although, yet again, local interruptions seem to occur when bottom hydrodynamics are locally anomalous.
The study area of Laganas Bay (NMPZ) is a unique area for exhaustive testing of seagrass interrelation to wave and light seafloor energies, allowing for the implementation of accurate predictive models, as its bathymetry and habitat extents are determined in detail through extensive application of hydracoustics and visual census mapping. Moreover, the NMPZ is a well-preserved and protected MPA, forced to minimum anthropogenic stressors, allowing the habitat extents to be controlled by natural forcings, free of any human intervention.
Given its success, this study is a stepping stone towards achieving more holistic habitat suitability models for finer-scale seagrass distribution, by incorporating light and hydrodynamics in the relevant modeling processes, extending the capabilities of more efficient monitoring and managing seagrass meadows in other MPAs in the Mediterranean Sea.