Predicting the distribution of threatened orbicellid corals in shallow and mesophotic reef ecosystems

Orbicellid corals are threatened primary reef-building corals throughout the Caribbean in shallow and mesophotic coral ecosystems (MCEs), yet a poor understanding of where they occur limits population monitoring and management. The goals of this study were to predict suitable habitat for orbicellid coral species and to identify how abiotic environmental factors constrain that habitat on the eastern Puerto Rico Shelf. The probability of occurrence for Orbicella annularis and O. faveolata/O. franksi (combined) from shallow to mesophotic depths on the eastern Puerto Rico Shelf was predicted using maximum entropy models. Contributions of abiotic predictors, including bathymetry, seafloor topography, temperature, wave exposure, and bottom velocity, were assessed. Model performance was assessed using area under the receiver operating characteristic curve, standard error of the replicate model runs, and mean absolute error. Both O. annularis and O. faveolata/O. franksi distributions were best predicted by rugosity, temperature, and wave exposure. O. faveolata/O. franksi occurred at shallow and mesophotic depths, and acted as a proxy for identifying the spatial extent of MCEs, contrary to O. annularis, which was predicted at shallow depths. Results for O. faveolata/O. franksi in mesophotic depths indicated potential for large areas of unexplored and unmonitored MCEs along the southeast shelf of St. Thomas, US Virgin Islands and within the Virgin Passage. These spatial predictions of potential mesophotic reef habitats will provide direction for future MCE exploration efforts.


INTRODUCTION
Caribbean coral reef-building is dominated by 2 genera, Acropora and Orbicella. These reef-building corals face multiple threats, including coastal development, land-based pollution, and climate change (Pandolfi et al. 2003). The branching acroporid corals have seen extreme population loss since the 1970s (Aronson & Precht 2001), leaving 3 massive orbicellid species (O. annularis, O. faveolata, and O. franksi) as the foundational species. However, orbicellid populations have also declined throughout the Caribbean (Jackson et al. 2014, Edmunds 2015, and their recent precipitous decline has led to all 3 species being classified as 'threatened' under the United States Endangered Species Act (79 FR 53852, 10 September 2014).
Orbicellids are primary reef-building corals (Goreau & Wells 1967, Reed 1985, Weil & Knowlton 1994 in shallow and mesophotic coral ecosystems (MCEs), which are found globally between 30 and 150 m depths. O. annularis primarily occupies shallow reef habitat (Pandolfi & Budd 2008) in depths less than 20 m (Weil & Knowlton 1994), whereas O. faveolata and O. franksi are common in depths up to 60 m throughout the Caribbean (Armstrong et al. 2006, Menza et al. 2008, Smith et al. 2008. In the United States Virgin Islands (USVI), mesophotic orbicellid banks are the dominant overall coral reef type, with cover and density surpassing shallow coral reef development (Smith et al. 2019). High cover orbicellid mesophotic reefs are primarily found along the southern insular shelf of the northern USVI (St. Thomas and St. John;Smith et al. 2010;see Fig. 1) and extend along the southern shelf of Puerto Rico. Mesophotic depths also extend out along the northern shelf area of St. Thomas, but these areas have sparse orbicellid and total coral cover (Groves 2016). Orbicellids are strongly structured by wave exposure and tend to form reefs in sheltered locations (Done 1983, Chollett & Mumby 2012. Due to their depth, mesophotic reefs are protected from strong wave action and provide this ideal sheltered habitat (Lesser et al. 2009), although Smith et al. (2016a) suggest that upper MCEs are also structured by storms. Other abiotic characteristics are associated with the structure of hardbottom mesophotic habitats including seafloor complexity as measured by metrics like rugosity and slope (Sheppard 1982, Sherman et al. 2010, temperature (Lesser et al. 2018), and light (Kahng et al. 2019. However, it is unknown if these characteristics are correlated with other factors that contribute to coral development. Given the importance of Orbicella spp. in the Caribbean, mapping their distribution in shallow and mesophotic depths is critical to guide effective management of dwindling populations. Characterization of MCEs is limited due to the logistical constraints of diver-based and remote observations at mesophotic depths (Hinderstein et al. 2010). In some areas of the Caribbean, there have been recent improvements in classifying well-developed mesophotic shelf areas (e.g. the insular shelf of the northern USVI; Costa et al. 2017), but the resulting maps were not designed to show the distribution of specific coral species, such as orbicellids.
Species distribution modeling (SDM) is a quantitative approach to predict the potential habitat of an organism across a landscape or seascape based on the abiotic variables that shape its distribution (Gui san & Zimmermann 2000, Franklin & Miller 2009. SDMs have been utilized to predict the distributions of threatened and endangered species (Godown & Peterson 2000, Engler et al. 2004, and can be used for monitoring declining populations, implementing restoration plans, managing and conserving critical habitat, or identifying unknown populations (Guisan et al. 2013). They have been applied to deep and shallow coral reef habitats (see Bridge et al. 2012, Chollett & Mumby 2012, Franklin et al. 2013, Georgian et al. 2014, Silva & MacDonald 2017. SDMs are also used to explore how abiotic factors influence species' distributions through variable contribution (i.e. how much does the abiotic variable contribute to the predicted distributions), and SDMs can be used to explore how variable contribution changes with different model spatial scales (resolution and extent) (Van Horn 2002).
Here we present a predictive species distribution modeling effort for orbicellid corals on the eastern Puerto Rico Shelf from shallow to mesophotic depths. The study objectives were to (1) develop SDMs for the threatened coral species O. annularis and O. faveolata/O. franksi (combined) using coral presence records and abiotic variables associated with orbicellid distributions to identify mesophotic coral habitat, (2) assess whether abiotic variable contributions vary within shallow and mesophotic depths or by geographic region (northeast and southeast shores), and (3) provide spatial predictive maps to help guide characterizations of potential mesophotic reef habitat for future field validation.

Study area
St. Thomas and St. John, USVI, are located on the eastern Puerto Rico Shelf at the interface of the northeastern Caribbean Sea and tropical western Atlantic. The shelf area also consists of the main island of Puerto Rico, Culebra, Vieques, and the British Virgin Islands. The shelf area encompassing the northern USVI (up to 65 m depth) is 1572 km 2 (Kadison et al. 2017). The study area included shallow to mesophotic depths (extending to 60 m) surrounding the islands of St. Thomas and St. John, USVI, and Culebra and Vieques, Puerto Rico (Fig. 1).

Study species
Orbicella spp. occurrence data were compiled from multiple shallow and mesophotic reef survey efforts conducted from 2007 to 2017. These datasets included (1) spatially stratified random drop camera and diver surveys characterizing the mesophotic depths of the north and south shores of St. Thomas (Smith et al. 2010, 2016a, Groves 2016depth > 30 m); (2) permanent transects from the USVI Territorial Coral Reef Monitoring Program (Smith et al. 2015;5−40 m); (3) (Herzlieb et al. 2006), occupy the same niche, and may hybridize (Szmant et al. 1997). Species occurrences outside the spatial extent of the predictor datasets were not included, and each cell (at a 50 m resolution) , and absences are displayed in addition to the number of presences and calculated prevalence for each coral group. Some areas of the study area have straight edges due to the size of some of the grids of the predictors. In particular, the regional ocean modeling system (ROMS) grid (bottom velocity), was square shaped and encompassed a large part of the area. The bathymetry layer was then clipped to the ROMS layer within the predictive framework was limited to a single occurrence (i.e. duplicate occurrences in each grid cell were removed) to reduce potential sampling biases for the model. For each coral group, prevalence was also calculated (number of presences / total number of presence and absence points combined) (see Figs. 1 & 2) to be used as a parameter in the models based on modeling methods conducted by . Using prevalence as a para meter in the models is further explained in Section 2.5.1. We tested the spatial autocorrelation of the presence points and found that they were spatially autocorrelated, which can bias model results (Dormann et al. 2007). However, after thinning the presence points, we found that we had an insufficient number of presence points to create accurate models. We discuss the limitations of spatially autocorrelated data in Section 4.2. Data compilation and cleaning were conducted in R software v. 3.3.3 (R Core Team 2017). Spatial autocorrelation analyses were conducted using the 'Spatial Analyst' toolbox in ArcGIS v. 10.5.1.

Abiotic variables
Abiotic variables tested during model development included bathymetry (depth), seafloor complexity, sea surface temperature (SST), benthic orbital velocity, wave height, wave period, wave direction, and bottom velocity ( Table 1). The extent of the study area was limited to 60 m depth to account for known depth limitations of Orbicella spp. (Smith et al. 2010). The extent of the study area was also limited by the available coverage of high-resolution bathymetry data and the extent of the bottom velocity data. All abiotic variables were resampled to a 50 m resolution. We modeled at a 50 m resolution to account for variability in GPS data and to match the resolution of the sampling grid used for the NOAA NCMRP data collection.
For the study area, bathymetry was mosaicked from multibeam sonar and light detection and ranging (LiDAR) data collected from 2008 to 2016 by NOAA (retrieved from https://www.ncei.noaa.gov/). The native spatial resolution of the bathymetry data ranged from 2 to 10 m, but we developed bathymetric variables at a resolution of 50 m. Seafloor complexity variables derived from bathymetry included aspect, standard deviation of the depth, rugosity, slope, change in slope, total curvature, planform curvature, and profile curvature at a 50 m resolution. Rugosity was calculated as the ratio of the surface area to the planar area across the neighborhood of a central pixel, derived from the compiled bathymetry data. Using this method, flat areas are equal to 1, with more rugose areas >1 (Jenness 2004, Wilson et al. 2007); this is closer to a measurement of relief rather than the rugosity produced by individual coral colonies. Bathymetric derivations were conducted using the 'Spatial Analyst' toolbox in ArcGIS v. 10.5.1, and the rugosity was derived using DEM Surface Tools (Jenness 2004(Jenness , 2013. We used the National Aeronautics and Space Administration Multi-Scale Ultra-High Resolution Sea Surface Temperature data (NASA MUR; https:// podaac. jpl.nasa.gov/MEaSUREs-MUR). Daily daytime SSTs were provided at 1 km resolution for 2011− 2017. Average, maximum, minimum, range, and summer average (June−September) temperatures were first calculated from the native 1 km resolution SST grids across the temporal range of the data. The resulting grids were then resampled to 50 m resolution using bilinear interpolation and clipped to the model extent.
Wave height, wave period, and wave direction were acquired from the Caribbean Coastal Ocean Observing System (CariCOOS; https://www.caricoos. org/ data-download) using the CariCOOS Nearshore Wave Model (Anselmi-Molina et al. 2012) and a Simulating Waves Nearshore spectral wave model (Booij et al. 1999) made twice daily from 2013 to 2016 at a 200 m native resolution. The average and maximum wave height and period (along with average wave direction) were calculated from the 200 m native resolution grids, and then resampled to a 50 m resolution using bilinear interpolation. The average and maximum benthic orbital velocities were also calculated across all grids (using the equation of Wiberg & Sherwood 2008). Bottom velocity was acquired from Regional Ocean Modeling System (ROMS) simulations based on in situ current measurements throughout the water column (see Cherubin et al. 2011). It was included to represent currents at depth driven by oceanographic variables other than surface waves. The ROMS model had a temporal range from 2007 to 2008 and native resolution of 330 m, and was resampled to a 50 m resolution using bilinear interpretation.

Model extents
Within our model domain, there is variability in environmental conditions between regions and depths ( Fig. 1). Environmental conditions differ be tween the northeast and southeast shores of this area (Groves 2016, Groves et al. 2018) as well as between shallow and mesophotic depths (Smith et al. 2016a (2) south shore only (referred to as '0− 60 m south shore model'); (3) shallow areas only (0− 30 m depth limit; referred to as 'shallow only model'); and (4) mesophotic areas only (30−60 m depth limit: referred to as 'mesophotic only model') ( Fig. 2). We removed the north shore of our study area because mesophotic coral development is sparse in this region (Groves 2016), and we hypothesized that the spatial predictions for Orbicella spp. and contributing environmental variables driving those predictions would change with the removal of this area. Additionally, the southeast shore of the study region contains highly developed and dense mesophotic reefs (Smith et al. 2010(Smith et al. , 2016b. We wanted to ensure that we were producing model results for the southeast shore region that were as ac curate as possible and would not be biased by lack of mesophotic reef development along the northeast shore. Similarly, we modeled at shallow depths and at mesophotic depths separately so those models would not be biased by including the full depth range. O. annularis was not included in the meso photic only model due to a low number of presences in the area. All abiotic variables were clipped to these 4 extents for modeling input, and all Orbicella spp. occurrence data that fell outside these extents were not included in the models.

Predicting orbicellid distributions
Maximum entropy modeling (MaxEnt v. 3.4.0; Phillips et al. 2017) was used to create spatially explicit predictions of O. annularis and O. faveolata/O. franksi distributions. MaxEnt is a presence/background sample SDM (sometimes referred to as presenceonly) that predicts habitat suitability based on correlation between species presences and environmental conditions and extends this correlation to other areas with similar conditions (Phillips et al. 2006(Phillips et al. , 2017. Maximum entropy modeling is a common method for estimating species distributions from presence data when absence records are unavailable or unreliable (Elith et al. 2011, Merow et al. 2013, and is commonly implemented using the Java software MaxEnt  (Phillips et al. 2006(Phillips et al. , 2017. MaxEnt uses a statistical machine learning algorithm to estimate the distribution (geographic range constrained by the extent of the environmental predictors) of a species by finding the distribution which has maximum entropy (Phillips et al. 2006(Phillips et al. , 2017. The MaxEnt formulation is used to maximize the likelihood of a parametric expo nential distribution (Phillips et al. 2004), and is mathematically equivalent to a Poisson point process model (Renner & Warton 2013). We chose to use MaxEnt because it is a common modeling technique for predicting the distribution of deep-water corals (Georgian et al. 2014, Kinlan et al. 2020) and mesophotic corals , Silva & MacDonald 2017. We also chose to use MaxEnt due to the abundance of presences for Orbicella spp. at mesophotic depths along the south shore, making it difficult to use other modeling methods. Because we were look-ing for areas of undiscovered mesophotic reef in this section of the Caribbean using Orbicella spp. as a proxy, we determined that MaxEnt models would be the best in identifying unsampled areas with the highest likelihood of finding orbicellids.

Model development
We developed our modeling methods by following the advice outlined by Merow et al. (2013) and the mesophotic modeling methods of Costa et al. (2015). Model development, performance, and evaluation were based on 10 replicate MaxEnt models for each coral group (O. annularis and O. faveolata/O. frank si). The parameters that were set in the MaxEnt models were as follows: random seed = on, replicated run type = subsample, replicates = 10, maximum iterations The random seed refers to the randomly chosen starting point used to subset the data for model training and cross validation. The regularization value (or beta multiplier) is a MaxEnt parameter used in machine learning to reduce overfitting of the model. The regularization value was determined using the 'Maxent-VariableSelection' package in R v. 3.3.3 (Jueterbock et al. 2016), which runs preliminary MaxEnt models and tests a series of beta multipliers. The beta multipliers tested in the preliminary model runs included β = 1, 2, 3, 5, 7, 11, and 13, and were tested for each of the 4 comparative models. The regularization value that yielded a minimum small-scale adjusted Akaike's information criterion (AICc) (Burnham & Anderson 2002) for each model can be found in Table S1 in the Supplement. We also ran preliminary models testing a different number of background samples (5000, 10 000, and 20 000), and found that 10 000 background samples yielded the smallest AICc. Of the 22 predictors, 8 created the best models for both coral groups. Images of the 8 predictors can be found in Fig. S1 in the Supplement. We used a pairwise correlation analysis to reduce the number of predictors. We re moved highly correlated predictors, and ran a series of preliminary MaxEnt models testing different combinations of predictors that created the best model based on the AICc and the area under the receiver operating characteristic (ROC) curve (i.e. area under the curve; AUC). We also chose predictors that were most ecologically relevant in driving Orbicella spp. distributions. The default prevalence calculated for each coral group based on the presence and absence data was 0.086 for O. annularis and 0.432 for O. faveolata/O. franksi for the 0−60 m model; however, the prevalence differed with changes in the spatial ex tent of the models (see Fig. 2 for the other prevalence calculations), and was included as a parameter in the MaxEnt models. An inherent issue with presence/ background modeling (e.g. MaxEnt) is that pre valence cannot be calculated from presence-only datasets. The prevalence parameter in MaxEnt models is arbitrarily set to 0.5 by default (Elith et al. 2011). Prevalence in this case is a measure of the sampling effort. It can be used for a more accurate interpretation of the logistic transformed spatial predicted outputs, and allows for a direct comparison between the coral groups that are modeled. The prevalence para meter used in these models is only estimated within the study area and does not represent prevalence throughout the entire distribution of a species (Elith et al. 2011). With a priori organism prevalence data, MaxEnt can be applied to specifically predict the probability of species occurrence rather than the habitat suitability (e.g. Elith et al. 2011. Data were randomly partitioned so that 70% of records were used for training and calibrating the model, and 30% were used to test model performance. For each replicate model run, different portions of the data were randomly subset for training and testing.
The resulting spatial predictive maps from MaxEnt depicted the complementary log−log (cloglog) transformed output, which represents the probability of occurrence on a scale of 0 to 100%, where 100% demonstrates that there is a 100% probability of the species occurring within that specified grid cell. The cloglog transformation gives a more appropriate estimate of the probability of species occurrence compared to the previous default output of MaxEnt, the logistic transformation (Phillips et al. 2017). The maximum training sensitivity plus specificity threshold (max SSS) calculated from the MaxEnt models was applied to the resulting continuous spatial predictions to create binary spatial predictive maps, which were used to identify the area the species is likely to occupy (Liu et al. 2013(Liu et al. , 2016.

Model performance
The discrimination capacities (the ability of the model to correctly distinguish between presences and background samples) of all models were assessed using ROC curves, which compare the sensitivity (true positive prediction rate) to 1 -the specificity (false positive prediction rate) over the range of predicted values. The AUC is the measure of the model's performance. AUC of 0.5 is indicative of the null model, meaning the model is no better at discriminating between presences and background points than random chance. Models with an AUC greater than 0.9 indicate superior model performance (Hosmer et al. 2013). ROC curves with test and training AUC values were generated by MaxEnt for each of the 10 model replicates and averaged.
Model performance was also assessed using reliability, the agreement between predicted and observed values (Liu et al. 2009). Reliability was represented by mean absolute error (MAE), a mea surement of the average magnitude of the predictive errors. MAE is defined as: ( 1) where n is the total number of species presences and true absences, y i is the predicted probability of , and x i is the observed species occurrence (0 or 1) based on a method used by Costa et al. (2015). The absence data were thinned similarly to the presence data, where 1 absence was recorded per grid cell to account for sampling biases. If a presence and an absence fell within the same grid cell, that cell was counted as a presence. The occurrences used to calculate MAE were the test presence points (30% of the data) set aside by MaxEnt for each of the 10-replicate model runs.
The true absence data available for the full modeling region were used (not the background points generated by MaxEnt

Model uncertainty
Model uncertainty was quantified to provide some indication of how sampling variability and model structure may have influenced the predicted distributions (Barry & Elith 2006, Lobo et al. 2008, Elith & Leathwick 2009). Uncertainty maps for the model replicates of each coral group were created by taking the standard error among the 10 replicate cloglog spatial distribution maps, and were assessed visually for areas of high uncertainty.

Contribution of predictor variables and response curves
Contributions of predictor variables to species distribution predictions were measured using the percent contribution of each predictor as calculated by MaxEnt. When the MaxEnt model is trained, it records which environmental variables make the greatest contributions to the model's predictions and converts these into a percentage. Percent contribution of variables was assessed across the 4 different model extents to identify how the abiotic variable contribution changed with depth or geographic region. Response curves were used to assess how probability of occurrence predictions changed in the context of the environmental variables.

Spatial predictions of occurrence
According to the average test AUC values, both the Orbicella annularis (mean ± SE: 0.955 ± 0.002) and O. faveolata/O. franksi (0.873 ± 0.002) models exceeded the random model prediction (AUC > 0.5). O. annularis had a higher probability of occurrence in shallow depths (<15 m) in the full depth range model (0−60 m; Fig. 3A), whereas O. faveolata/O. franksi occurred more widely in both shallow and mesophotic depths (Fig. 4A). O. faveolata/O. franksi had the highest probability of occurrence in the Virgin Passage, which is an area of reef habitat between the islands of St. Thomas, USVI, and Culebra, Puerto Rico. The max SSS threshold identified more potential habitat area on the southeast shelf compared to the northeast shelf (Fig. 4B). Based on the max SSS threshold area calculations, O. annularis was likely present in 474 km 2 (13.6% of the 3497 km 2 total survey area), and O. faveolata/O. franksi was likely present in 604 km 2 (17.3% of the total survey area) ( Table 2). Modeling results for the other model extents can be found in 'Model extent analyses' and Table S2, both in the Supplement.

Spatial prediction uncertainty and error
The MAE results and uncertainty maps demonstrated greater discrepancies between test presences/ true absences and predicted probabilities for O. faveolata/O. franksi. Model uncertainty was high based on visual assessments of the model uncertainty maps (Figs. 3C & 4C), indicating high variability be tween the model runs. Based on visual assessments, the highest uncertainties for each coral group were co-located with the highest probability of occurrence values: in the Virgin Passage and along portions of the southern insular shelf for O. faveolata/O. franksi (Fig. 3C), and in shallow areas for O. annularis (Fig. 4C) Fig. 6). Depth (32.96%) and SST range (6.01%) both followed rugosity as important contributors to O. annularis predictions. Orbital velo city (14.24%) and average SST (8.70%) followed rugosity as contributing variables for O. faveolata/O. franksi.

Orbicella annularis
Varying spatial extents of abiotic predictors along the eastern Puerto Rico Shelf did not affect predictions of O. annularis. The comparative models showed little change in percent contribution of all variables, with rugosity, depth, and SST range re maining the top 3 predictors for the 0−60 m, south shore and shallow only models (Fig. 5). The rugosity response curves for O. annularis indicated that probability of occurrence increased with increasing ru gosity (Fig. 7). The probability of occurrence tapered off when a rugosity value of 1 was surpassed. Rugosity values >1 are indicative of high-relief habitat, with 1 representing flatter areas. The variation in rugosity no longer mattered in predicting probability of occurrence after this 'threshold' was reached. In addition, the probability of O. annularis occurrence decreased with increasing depth and decreased with variability in SST range (Fig. 7).

Orbicella faveolata/O. franksi
In contrast to O. annularis, the variables that best predicted occurrence of O. faveolata/O. franksi changed with the extent of the model. The mesophotic only model for O. faveolata/O. franksi was markedly different from the other models. For the 3 models that included shallow water extents (0−60 m full, 0−60 m south shore, and shallow only), rugosity was by far the most important variable. In the mesophotic model (30−60 m), wave direction became the most important variable, followed by rugosity and average SST (Fig. 6). The contribution of the second and third most important variables in the models that included shallow water were minor relative to rugosity, but varied slightly between models. Rugosity, wave direction, average SST were the most important predictors in the 0−60 m full model. Rugosity, orbital velocity, and average SST were the 3 predictors contributing most to the probability of occurrence in the 0−60 m south shore model. Lastly, rugosity, orbital velocity, and SST range were the top 3 contributing variables in the shallow only model.
As with O. annularis, the response curves indicated an increase in rugosity above a value of 1, which was associated with increased probability of occurrence for O. faveolata/O. franksi (Fig. 7). The response curves for wave direction showed a strong shift from low probability at 45°t o a peak in probability be tween 80  at the lower SST range, but dropped off markedly above 27.7°C. The orbital velocity response curve showed that high probability of occurrence peaked with very low orbital velocity, approaching 0.

DISCUSSION
The resulting spatial predictions indicated that unmonitored and unmanaged areas of Orbicella spp. are widespread in both shallow and mesophotic coral reef habitat along the eastern Puerto Rico Shelf. In particular, there is the potential for large areas of previously unidentified O. faveolata/O. franksi reef habitat at mesophotic depths. However, some of the model performance measures indicated areas of uncertainty and disagreement between spatial predictions and true occurrences. The abiotic variables that shaped these predicted distributions included rugosity, depth, wave exposure, and temperature.

Rugosity
Rugosity was the strongest predictor for O. annularis and O. faveolata/O. franksi across several models, with increasing relief (rugosity >1) indicating an increased probability of occurrence. The response curves indicated that corals were excluded from flat areas. These flat areas are indicative of soft substrate, sediment, or depositional environments, and are not ideal habitat for coral settlement and growth. Corals could also be responding to high relief substrate in regards to settlement and subsequent survival. High relief substrate is a driver of coral distributions in Hawaii because antecedent reef structure (i.e. not created by recent reef building) creates ideal substrate (Jokiel et al. 2004), and similarly, antecedent reef buttresses also create ideal habitat for  (Sherman et al. 2010). In seafloor mapping studies, high rugosity values were indicative of hard coral occurrences in a shallow reef area of Puerto Rico (Prada et al. 2008) and in shallow and mesophotic coral reefs in Bonaire (Trembanis et al. 2017). Shallow coral occurrence was also predicted by both slope and rugosity in Hawaii (Franklin et al. 2013). Spatial resolution is a key consideration for seafloor topography inclusion in SDMs (Pittman et al. 2009, Lecours et al. 2015. In our study, the rugosity measurement used in the models is more similar to a measurement of relief, and not a measurement of individual coral colonies that create the rugosity. Corals create rugosity (Alvarez-Filip et al. 2009, Graham & Nash 2013, and Orbicella spp. are primary reef-building corals in shallow and mesophotic reef habitats (Weil & Knowlton 1994, Newman et al. 2015. Orbicella spp. distributions were modeled at a 50 m resolution (i.e. much greater than the sizes of individual coral colonies), and the spatial scale sug-gests other types of rugose substrate, such as antecedent reef, are creating this ideal habitat. The resolution of the rugosity layer was derived from the 50 m bathymetric grid. The importance of rugosity in these models reflects the overall macro-structure and general topography of the reef, which could suggest that the antecedent reef structure is driving the model re sults. Shallow and mesophotic reefs are built upon drowned reefs submerged from sea level rise during the Pleistocene epoch (Lesser et al. 2018), and this antecedent reef could be providing the ideal topographic complexity and high relief habitat for Orbicella spp., similar to conditions in Hawaii (Jokiel et al. 2004) and Puerto Rico (Sherman et al. 2010).
Rugosity could also be acting as a proxy for the angle of available substrate, or slope, as they are positively correlated (dependent on spatial scale) (Pittman et al. 2009). Higher relief or more complex habitat is indicative of more angled surfaces. Although slope was tested in preliminary SDMs, it was not included in final model runs because the inclusion of 74 Fig. 6. As in Fig. 5, but for Orbicella faveolata/O. franksi (but includes the mesophotic only model [M]). Note differences in the y-axes for the variables, as they do not contribute equally to models rugosity in the models produced higher AUC values. Seafloor slope affects sedimentation, which also structures MCEs (Sherman et al. 2016). High-relief areas favor the downslope transportation of sediment away from corals, improving coral development (Sher man et al. 2010). Slope and rugosity are also important for coral recruitment (Lesser et al. 2018). We suggest that rugosity is a strong predictor because it is a positive feedback loop for coral reef growth, i.e. rugosity encourages coral recruitment, which creates more rugosity. In shallow reef habitat, coral recruits prefer downward-facing substrate, and in deep areas, upper horizontal surfaces are preferred (Harborne et al. 2006). In studies using settlement tiles to observe coral recruitment, recruits also prefer rugose texture and refuges (Edmunds et al. 2014), although our models did not pick up on such fine-scale rugose features. The preference for more complex habitat reflects the preference of a coral for low light levels, low sedimentation rates, and protec-tion from competition and predation (van Woesik et al. 2014). Areas without relief or topographic complexity can limit coral recruitment due to lack of shelter, excessive algal growth, and the presence of sediment (Edwards & Clark 1999). Rugose surfaces also lead to 'sticky water,' an oceanographic process where areas of high reef density are poorly flushed because the prevailing currents are directed away from these regions, leading to high levels of natal recruitment at well-developed reefs (Andutta et al. 2012).

Wave direction
Wave direction was a strong predictor of O. faveolata/O. franksi in the shallow and mesophotic models that included the north shelf habitats. This likely represents the strong spatial structuring of mean wave climate between the north and south shelves of the northern USVI. The north shelf is exposed to the Atlantic swell and sees a dominant wave climate from the northeast (45°), whereas the south shelf is sheltered from the Atlantic swell and has a wave climate dominated by regional trade winds from the east (90°). These dominant wave directions correspond to the low and high probability occurrence predictions in the models for O. faveolata/O. franksi. The wave direction importance in the models diminishes when the north is removed; thus, in the northinclusive model, the wave climates are likely separating the low occurrence mesophotic orbicellid areas of the north with the high occurrence mesophotic orbicellid areas of the south ( Fig. 1 Wave direction could also be a proxy for food availability. Species-specific changes from autotrophy to heterotrophy occur from shallow to mesophotic depths (Lesser et al. 2018). Montastraea cavernosa becomes more heterotrophic at mesophotic depths, relying more on zooplankton and inorganic nutrients (Lesser et al. 2010). Potential food sources are transported onto the reef through upwelling or other sources. There is evidence from the Indo-Pacific to suggest that deep upwelling and wave direction drive food availability, which further drives coral trophic zonation (Williams et al. 2018, Radice et al. 2019).

Temperature
SST was also a contributing abiotic variable across all model extents. Generally, the response curves showed a pattern of decreased probability of occurrence with increasing temperatures or more variable temperature ranges. These environmental variable preferences are in close agreement with the findings of Castillo & Helmuth (2005), who reported that O. annularis did not respond well physiologically to acute changes in temperature, and tended to prefer cooler waters in shallower areas. The SST range predictor did not capture acute changes in temperature, because it was averaged across several years; however, the model results provided insight into the habitat preferences for O. annularis based on low variability in temperature fluctuations. Changes in temperature over long periods of time have also contributed to the decline of O. annularis cover in Florida and the Caribbean (St. John, USVI, Edmunds & Elahi 2007;Florida Keys, Ruzicka et al. 2013), providing further evidence that O. annularis requires a certain ideal temperature range to thrive. Temperature plays a key role in differentiating between shallow and mesophotic reef habitats as well, and is a limiting factor in determining the lower depth limits of MCEs (Bongaerts et al. 2015). MCEs are cooler than shallow reefs (Lesser et al. 2009, Smith et al. 2016b, Lesser et al. 2018) most likely from internal oceanographic processes (Bak et al. 2005). Smith et al. (2016b) identified that benthic temperatures in Orbicella spp. mesophotic sites were cooler than shallow coral reefs in the USVI specifically, which could explain the difference in temperature contributions for the shallow and mesophotic only models. Orbicella spp. at meso photic depths are more sensitive to elevated temperatures, and have likely acclimated to temperatures at these depths, thus making cooler temperatures ideal for Orbicella spp. (Smith et al. 2016b). We used SST data from 2011−2017 to closely match the data collected for presences of Orbicella spp., which is common in the coral modeling literature (see Georgian et al. 2014, Chen et al. 2020). These current Orbicella spp. presences may be driven by current temperature conditions, but increases in temperature due to climate change are already causing, and will continue to cause, negative impacts to coral populations.
Benthic temperature data were not used in the MaxEnt models because they were not available for the entire model extent. Winship et al. (2020) recommended using in situ temperature data to create an interpolated temperature grid at depth for deep coral SDMs. While some in situ temperature data were available for the study region, there were not enough data to accurately create a bottom temperature grid. The use of SST could be a proxy for benthic temperatures instead. Costa et al. (2015) identified that SST in Hawaii was a proxy for benthic temperature at mesophotic depths, and similar results were seen in other parts of the Pacific (Kahng et al. 2012). In the Caribbean, comparisons between satellite-measured temperature and in situ temperature measurements have only been documented within the first 20 m of water (Castillo & Lima 2010, Neal et al. 2014, Gomez et al. 2020. For the purposes of this study, correlations with benthic and satellite-measured temperatures would potentially make SST an ideal predictor to use in modeling species distributions at shallower depths. It is unknown if this would work similarly at mesophotic depths in the Caribbean and is something to consider for future study.

Model error and uncertainty
The distributions of both coral groups are potentially inaccurate in some areas, given the high MAE values. The MAE values were high when calculated separately for the presences and absences (rather than combined) for O. faveolata/O. franksi. High MAE for presences is indicative of underpredictions from the model, whereas high MAE for absences is indicative of overprediction. The high MAE exhibited may be due to missing or inaccurate environmental predictors such as backscatter and light-attenuation grids (Kahng et al. 2019. Environmental data are sparse for marine environments compared to terrestrial environments (Robinson et al. 2011). Most datasets are available for the surface and subsurface levels, and rarely for deeper depths (Lecours et al. 2015). Several predictors used in these SDMs were surface values at coarse resolutions (1 km), and may not be indicative of coral distributions at mesophotic depths. Several of the predictors were also resampled from coarser (1 km) to finer resolutions (50 m), further introducing error.
High-resolution multibeam or LiDAR provides bathymetry and rugosity data that are key habitat descriptors to include in marine SDMs (Valavanis et al. 2008, Pittman et al. 2009). Acoustic backscatter data differentiate between hard-and softbottom habitats (Anderson et al. 2008), but were not available for the entirety of the region. However, combining backscatter data from different multibeam datasets could prove problematic due to the subjective way backscatter is processed. It is important to differentiate between hard-and softbottom, but not all hardbottom habitat is ideal for coral settlement and growth. The upper mesophotic zone along the southern insular shelf contains a variety of hardbottom habitat including horizontal banks, steep slopes, and walls. Orbicella spp. occur on horizontal banks, whereas Agaricia spp. are common on steep slopes and walls. The upper mesophotic zone also supports large areas of softbottom habitat, which is largely composed of rhodolith beds. These beds contain small stony coral patches, including the sparse occurrence of orbicellid colonies (Smith et al. 2019). Parsing out different kinds of ideal habitat for orbicellids would be advantageous in future modeling efforts and would also contribute greatly to a regionwide benthic habitat map. High-quality benthic habitat maps cover most of the modeling area (see Kendall et al. 2001, Kågesten et al. 2015, Costa et al. 2017); however, a benthic habitat map encompassing the entire model extent has not been created due to lack of multibeam and acoustic backscatter data for parts of the northeast shore of the northern USVI.
Another source of error and uncertainty in our models comes from the spatially autocorrelated Orbicella spp. presences, which can greatly bias model results (Dormann et al. 2007). Spatial autocorrelation of the data is most likely due to the heterogeneous nature of the compiled datasets, using data from multiple different sources with varying survey designs as outlined in Section 2.2. Future modeling ef forts could utilize the methods of Georgian et al. (2019), who calculated a residual autocovariate (RAC) variable for each taxonomic group they modeled. The RAC method derives a term representing the spatial autocorrelation present in the residuals of a preliminary model. This term is then added as an explanatory variable in the subsequent models and reduces modeling error.

Conclusion
The results of this study provide predicted distributions of Orbicella spp. and mesophotic reef locations that are located along the eastern Puerto Rico Shelf on a broad scale under current environmental conditions. Shallow and mesophotic reef ecosystems have the potential for large areas of unexplored, unmonitored, and unmanaged areas of orbicellid-dominated reef habitat. The spatial predictive maps now need to be ground-truthed, in conjunction with uncertainty maps, to further characterize mesophotic orbicellid communities, which would lead to more effective monitoring of these declining species. The SDMs could be refined, with the collection of environmental data at higher spatial resolutions, larger spatial extents, the addition of new environmental predictors, and through the use of more refined modeling techniques (i.e. presence/absence or density model-ing). Multibeam, acoustic backscatter, benthic habitat maps, and light attenuation are missing for some of the region, and collection of these data will greatly refine the orbicellid models. Gridded benthic temperature data are needed for the entire modeled region to assess whether satellite-derived temperature is correlated with benthic temperature. This in turn could provide more information about temperature as a driver in orbicellid distributions. Refined models would theoretically lead to more accurate spatial predictive maps, which would further lead to a better understanding of these unexplored orbicelliddominated reefs.