Habitat Suitability Modeling to Inform Seascape Connectivity Conservation and Management

Abstract

Coastal habitats have experienced significant degradation and fragmentation in recent decades under the strain of interacting ecosystem stressors. To maintain biodiversity and ecosystem functioning, coastal managers and restoration practitioners face the urgent tasks of identifying priority areas for protection and developing innovative, scalable approaches to habitat restoration. Facilitating these efforts are models of seascape connectivity, which represent ecological linkages across heterogeneous marine environments by predicting species-specific dispersal between suitable habitat patches. However, defining the suitable habitat patches and migratory pathways required to construct ecologically realistic connectivity models remains challenging. Focusing on two reef-associated fish species of the Florida Keys, United States of America (USA), we compared two methods for constructing species- and life stage-specific spatial models of habitat suitability—penalized logistic regression and maximum entropy (MaxEnt). The goal of the model comparison was to identify the modeling algorithm that produced the most realistic and detailed products for use in subsequent connectivity assessments. Regardless of species, MaxEnt’s ability to distinguish between suitable and unsuitable locations exceeded that of the penalized regressions. Furthermore, MaxEnt’s habitat suitability predictions more closely aligned with the known ecology of the study species, revealing the environmental conditions and spatial patterns that best support each species across the seascape, with implications for predicting connectivity pathways and the distribution of key ecological processes. Our research demonstrates MaxEnt’s promise as a scalable, species-specific, and spatially explicit tool for informing models of seascape connectivity and guiding coastal conservation efforts.

1. Introduction

Understanding the spatial, temporal, and environmental drivers of marine species distributions is paramount to developing ecologically sound conservation and place-based management strategies [1,2]. These efforts are increasingly urgent as global climate change interacts with local and regional stressors (e.g., pollution, overexploitation) to degrade and fragment marine ecosystems [3,4]. Coastal habitats, including mangroves [5], seagrasses [6], corals [7], and salt marshes [8], have experienced precipitous declines worldwide over the last several decades, thereby vastly reducing the suitable habitat space for many marine species. Biodiversity loss stemming from environmental degradation enhances the risk of functional collapse by reducing ecosystem resilience to interacting stressors [9,10]. Identifying priority locations for protection and restoration of coastal habitats at both local and global scales is therefore one of the most important challenges in marine conservation and spatial planning [11,12]. Seascape connectivity is increasingly recognized as a dynamic and spatially explicit process regulating biodiversity patterns and ecosystem functions, and playing an important role in guiding coastal conservation efforts [13,14,15,16].

Seascape connectivity represents ecological linkages and flows across heterogeneous environments [17,18], typically related to the flow of organisms, energy, nutrients, or genetic material. The magnitude and location of these exchanges are shaped by patterns of intra- and inter-habitat connectivity, with the latter being especially important for species whose movements span multiple habitat types. In tropical coastal seascapes, for example, multi-habitat use is widespread with almost half of all coral reef associated fish species having also been recorded in two or more non-reef habitats [19]. Seascape connectivity can be categorized into two major types—structural connectivity and functional connectivity [13]. Structural connectivity describes the spatial patterns and relationships of the seascape itself (e.g., patch areas, inter-patch distances, habitat corridors), whereas functional connectivity describes the degree to which animals move among resource patches in response to those structural patterns [18,20]. Functional connectivity is therefore inherently species- and life stage-specific, as it depends on the behavioral and life history traits of the organism, and on the spatiotemporal scales of their movements [18,21]. Thus, management programs that maximize functional connectivity across scales and communities, including through place-based conservation measures such as marine protected areas, are expected to achieve greater ecological outcomes [13,14,15]. Facilitating these efforts are models of potential connectivity, in which limited data on species behavior or dispersal are related to metrics of seascape structure [20]. These include graph-theoretic approaches, in which the seascape is represented by a spatial graph constructed of suitable habitat patches (nodes) connected by a series of dispersal links (edges) [22]. In this regard, potential connectivity models can be used to estimate connectivity thresholds, movement corridors, and stepping-stones unique to individual seascape residents, which can then serve as targets in conservation and restoration planning [23,24].

To successfully inform conservation efforts, connectivity models must achieve a high level of ecological realism, which relies on knowledge of the current distribution of suitable habitats and species–environment relationships. One way to access this information is through habitat suitability modeling, also referred to as species distribution modeling, a method that uses known species locations and their associated environmental conditions to predict habitat suitability over space and time [25]. When used as a precursor to connectivity modeling, predictive habitat suitability models (HSMs) can reveal the locations and sizes of suitable habitat patches, which can then serve as nodes in spatial graphs [26]. Furthermore, assuming that animals travel along pathways that minimize their ecological costs, HSMs can be inverted to produce cost surfaces over which connectivity is predicted (e.g., least-cost path models), potentially increasing the realism and precision of predicted edges [27]. As the integration of three-dimensional surfaces alongside standard two-dimensional predictors becomes more common in spatial modeling [28,29], the value of HSMs and HSM-derived cost surfaces for connectivity analysis will likely increase. Although realistic HSMs can help bridge the gap between connectivity modeling and conservation planning, applications of this multi-step approach have been applied only to terrestrial landscapes [26,30,31].

The limited availability and quality of species occurrence and environmental data represent major barriers to habitat suitability modeling in complex coastal seascapes. However, recent advancements in remote sensing technology have increased the accessibility of data on marine habitat mosaics, three-dimensional seascape terrain structure, and oceanographic conditions, which can serve as ecologically relevant spatial predictors in HSMs [28,29,32]. When possible, systematic surveys (e.g., SCUBA censuses) can be performed to collect detailed presence–absence data for focal species, which can then be used to predict the probability of species presence via generalized linear or additive modeling (e.g., logistic regression) [25,33]. The cost of systematic surveys, however, may limit their coverage and usefulness for seascape-wide modeling. In contrast, presence-only species data may be compiled from a variety of sources covering broad geographic and environmental space, including online biodiversity databases (e.g., gbif.org, obis.org), fishery-dependent and fishery-independent programs, and citizen science initiatives. When coupled with information on environmental conditions at a set of background points (i.e., presence-background data), these data can be used to estimate the relative likelihood of species occurrence [34]. MaxEnt, an open-source machine learning software that uses the principle of maximum entropy to model species distributions, is a popular presence-background technique [35,36]. Indeed, MaxEnt and standard regression techniques have displayed similar predictive performance in several terrestrial studies [37,38,39,40], but whether this outcome holds in complex seascapes and under the scrutiny of connectivity modeling remains unclear.

We compared penalized logistic regression and MaxEnt models of habitat suitability for two common fish species across a spatially heterogeneous coastal seascape in the Florida Keys, United States of America (USA). Our primary objective was to determine which modeling method produced the more realistic HSM for use in subsequent connectivity assessments. To meet this objective, we examined each model’s ability to discriminate between suitable and unsuitable locations using both threshold-independent and threshold-dependent assessments. Our secondary objective was to identify the most influential environmental and spatial predictors of habitat suitability for each species to better understand the species-seascape interactions that shape patterns of connectivity.

2. Materials and Methods

2.1. Study Area

The Florida Keys—a string of islands located off the southern tip of Florida between the Atlantic Ocean and the Gulf of Mexico—were selected as a case study for this research (Figure 1). The seascape in this region has been described as a mangrove–seagrass–reef continuum featuring patches of shoreline mangroves, dense seagrass beds interspersed with patch reefs, and finally the Florida Reef Tract located 5–15 km offshore [41]. Despite protection from the Florida Keys National Marine Sanctuary (FKNMS) and Biscayne National Park (BNP), this region has experienced significant declines in biogenic habitat over the last 50 years. A range of local, state, and federal organizations are now focusing substantial resources on habitat restoration to combat loss and fragmentation, with the offshore reef tract being a major target [42]. These restoration efforts will benefit from a thorough understanding of the daily, seasonal, and ontogenetic migrations of reef fishes between spatially isolated habitats across the seascape, as these movements increase reef productivity and resilience through the transfer of nutrients, the maintenance of top-down control on coral predators, and the enhancement of grazing pressure on harmful epiphytes and macroalgae that could otherwise shift reefs into a macroalgal state [43,44,45]. In fact, the transfer of nitrogen and phosphorus from sheltering resident and migratory reef fishes to their host corals has been demonstrated to significantly increase zooxanthellae abundance and coral growth rates [46,47]. Fish-derived nutrient hotspots thus offer important supplements to oligotrophic tropical and subtropical coral reefs, including those of the Florida Keys [48]. Therefore, we focused on modeling habitat suitability for reef fishes in the coastal region (≤50 m depth) from Key Biscayne in the Upper Keys to Cudjoe Key in the Lower Keys.

2.2. Focal Species

We selected gray snapper (Lutjanus griseus) and bluestriped grunt (Haemulon sciurus) as focal species due to their complex, multi-habitat life histories. Though generally considered reef-associated as adults, these species occupy spatially discrete patches of varying habitat type through ontogeny, resulting in them being categorized as seascape nursery species [49,50]. In Southern Florida, larval L. griseus and H. sciurus settle in seagrass beds around September and February, respectively, before expanding to mangroves several months later [51]. Juveniles often remain in nearshore nurseries for months to years, then, as sub-adults, undertake a cross-shelf migration to join adults on offshore coral reefs. Although it is clear that their ontogenetic migrations play an important role in maintaining seascape-wide ecological connectivity, there remains a paucity of information on habitat suitability across the seascape for L. griseus and H. sciurus. Additionally, there may be considerable inter-species variation in habitat suitability, and consequently, functional connectivity, stemming from unique preferences for nearshore strata [51,52], tolerances to salinity fluctuation [53], and motivations for movement [54,55]. Addressing these knowledge gaps is of critical importance as L. griseus and H. sciurus play key ecological and economic roles in the region as abundant members of the fish assemblage, highly mobile predators and vectors of nutrient transfer, and valued sport and commercial fishery targets [56,57].

Data describing the size, abundance, and distributions of L. griseus and H. sciurus were obtained from two multi-agency monitoring programs coordinated by the NOAA Southeast Fisheries Science Center: the South Florida Reef Visual Census (RVC) and the Mangrove Visual Survey (MVS). Using a two-stage stratified random sampling design, the RVC program surveys fish communities on coral reefs and other hard bottom habitats biennially using a stationary visual survey method [58]. Similarly, the MVS program conducts annual belt transect surveys alongside randomly selected mangrove shoreline sites in Biscayne Bay, Card Sound, Barnes Sound, and northeastern Florida Bay [59,60]. At each site, a trained diver records the identity, abundance, and size structure of all fishes encountered; RVC surveys occur within a 7.5 m-radius imaginary cylinder extending vertically from the seafloor to the surface, while MVS surveys take place over 30 × 2 m transects established parallel to the shore. These programs were designed to collect data that enables the detection of spatiotemporal changes in reef fish species composition, size, and abundance using statistical analyses. Thus, we compiled RVC and MVS data collected at unique sampling sites in 2014, 2016, and 2018 to maximize the spatial coverage of georeferenced L. griseus and H. sciurus records. This process also ensured a temporal match between the two reef fish monitoring programs and the spatial predictors described below, which are based largely on data collected between 2014 and 2018. Considering the dynamic nature of the Florida Keys, where seafloor features change over time, we used temporal data alignment to prevent model inaccuracies.

2.3. Spatial Predictors

To explore the relationship between spatial predictors and L. griseus and H. sciurus distributions, we constructed raster data layers to map the following environmental categories: benthic habitat, bathymetry and seafloor morphology, and water conditions (

Table 1. Spatial predictors available for habitat suitability modeling. Grayed-out predictors were removed due to collinearity issues. BTM: Benthic Terrain Modeler ArcGIS extension (v3.0). Calculations performed using a 3 × 3 moving window of cells, unless otherwise noted.

	Predictor	Units	Description	Tool Used
Habitat	Benthic Habitat	Categorical (12)	Bottom habitat type	Not Applicable
	Distance to Mangrove	Meters	Euclidean distance to the nearest mangrove habitat	gridDistance function in the raster R package
Bathymetry & Seafloor Morphology	Depth	Meters	Water depth in each cell	Not Applicable
	Depth (Standard Deviation)	Meters	Local dispersion	Calculate Metrics tool in BTM
	Slope	Degrees	Rate of maximum change in depth	Slope tool in BTM
	Curvature	$\frac{1}{100}$th of a meter, convex (+) or concave (−)	Second derivative of the bathymetric surface	Curvature tool in ArcGIS Spatial Analyst
	Plan Curvature	$\frac{1}{100}$ th of a meter, convex (+) or concave (−)	Curvature perpendicular to the direction of maximum slope	Curvature tool in ArcGIS Spatial Analyst
	Rugosity	Ratio	Local surface roughness calculated as the ratio of surface area to planar area	Surface Area to Planar Area (slope-corrected) tool in BTM
	Broad-Scale Bathymetric Position Index (BPI)	Ridge (+), Flat (0), or Valley (−)	Depth of a cell relative to its surroundings, evaluated using concentric rings of 125 m and 1250 m	Broad-Scale BPI tool in BTM (standardized)
	Fine-Scale BPI	Ridge (+), Flat (0), or Valley (−)	Depth of a cell relative to its surroundings, evaluated using concentric rings of 5 m and 125 m	Fine-Scale BPI tool in BTM (standardized)
Water Conditions	Winter Temperature	Degrees Celsius	Mean winter (January–March) temperature	Krige function in the gstat R package
	Winter Salinity	Practical Salinity Units	Mean winter (January–March) salinity	Krige function in the gstat R package
	Winter Dissolved Oxygen	Milligrams per Liter	Mean winter (January–March) dissolved oxygen	Krige function in the gstat R package
	Summer Temperature	Degrees Celsius	Mean summer (July–September) temperature	Krige function in the gstat R package
	Summer Salinity	Practical Salinity Units	Mean summer (July–September) salinity	Krige function in the gstat R package
	Summer Dissolved Oxygen	Milligrams per Liter	Mean summer (July–September) dissolved oxygen	Krige function in the gstat R package

, Figure 2). All rasters had an interpolated resolution of 5 × 5 m and were referenced to the Florida East projected coordinate system and the NAVD vertical coordinate system where applicable.

2.3.1. Benthic Habitat

Benthic habitat data were obtained from Florida’s Unified Reef Map v2.0 [61], a seamless map of benthic habitats from Martin County to the Dry Tortugas derived from remote sensing imagery, high-resolution bathymetric data, and in situ observations using a five-tier hierarchical classification system. We supplemented the Reef Map with a separate GIS dataset of shoreline mangroves to fully capture the extent of this potentially important nursery habitat [62]. For this research, we used Level 1 map classifications, which delineate the major benthic habitat types while maintaining agreement between the Reef Map’s contributing agencies. Habitat data were rasterized, producing cells that reflected the IDs of 12 possible benthic habitats: individual or aggregated patch reef, scattered coral and rock in unconsolidated sediment, continuous seagrass, discontinuous seagrass, unconsolidated sediment, aggregate reef, pavement, reef rubble, mangrove, artificial, dredged and excavated, and ridge. We also assessed the importance of mangrove nursery proximity by constructing a raster of Euclidean distances to the nearest mangrove habitat using the raster R package (v3.3–13) [63] in R version 4.0.2 [64].

2.3.2. Bathymetry and Seafloor Morphology

NOAA’s 1/9th ArcSecond Resolution Continuously Updated Digital Elevation Model served as an initial bathymetric surface for this research [65]. Seafloor morphology was quantified from the bathymetric surface by applying metrics of slope, curvature, rugosity, bathymetric position index (BPI), and standard deviation of depth using the Benthic Terrain Modeler ArcGIS extension (BTM v3.0) [66] in ArcGIS v10.7.1. These rasters captured detailed information on the structure and complexity of the seafloor, including the locations of crests, flats, and valleys, the direction of benthic flow, and the roughness of the local surface [67]. The influence of seascape surface morphology on the distribution of tropical reef fishes has been demonstrated in coral reef ecosystems in the U.S. Caribbean [68,69], Hawaii [29,70,71], and elsewhere [72]; however, this approach has not been applied at a seascape scale to fishes of the Florida Keys.

2.3.3. Water Conditions

Seasonal water conditions, in addition to habitat availability and seafloor morphology, may play a role in shaping the distributions of L. griseus and H. sciurus. We investigated this possibility using raster data layers of mean bottom temperature, salinity, and dissolved oxygen. Using the gstat R package (v2.0–6) [73], we interpolated water quality data from sites that were sampled regularly over the 2014–2018 period by the Southeast Environmental Research Center’s Water Quality Monitoring Network and the Miami-Dade County Surface and Groundwater Quality Viewer. Our research focused on two critical seasons—winter (January–March) and summer (July–September)—to capture not only annual extremes, but also important periods in the life histories of L. griseus and H. sciurus, including habitat expansion and spawning activity [51,74,75].

2.4. Model Development

2.4.1. Filtering and Partitioning of Occurrence Records

As functional connectivity for L. griseus and H. sciurus across the Florida Keys seascape is maintained primarily by the cross-shelf (5–15 km) ontogenetic migrations of their sub-adult life stage, we focused our habitat suitability modeling efforts specifically on this subpopulation. Therefore, we restricted L. griseus and H. sciurus records using the size cut-offs defined in a previous Florida study [51], where sub-adults are those between the size at year 1 and the size at maturation (

Table 2. Sub-adult gray snapper (Lutjanus griseus) and bluestriped grunt (Haemulon sciurus) occurrence records in the southern Florida study area. Prior to habitat suitability modeling, data were randomly partitioned into calibration (70%) and evaluation (30%) subsets.

	Size (cm TL)	Presence Sites	Absence Sites	Total
Calibration Data
Lutjanus griseus	9.51–24.71	378	1129	1507
Haemulon sciurus	11.90–25.33	447	1060	1507
Evaluation Data
Lutjanus griseus	9.51–24.71	167	479	646
Haemulon sciurus	11.90–25.33	216	430	646

). Prior to habitat suitability modeling, species occurrence records were partitioned into calibration (70%) and evaluation (30%) subsets following a random split approach (Table 2).

2.4.2. Selection of Spatial Predictors

Pearson correlation coefficients (r) and variance inflation factors (VIF) were used to assess collinearity among the environmental raster data layers using thresholds of |0.7| and 5 for r and VIF scores, respectively [76]. Four spatial predictors were removed due to multicollinearity issues—standard deviation of depth, plan curvature, and dissolved oxygen across both seasons. Of the 12 predictors retained for modeling, the highest correlation was between summer and winter salinities (r = 0.61) and the largest VIF score, belonging to depth, was 2.63 (Supplementary Table S1 and Supplementary Figure S2).

2.4.3. Penalized Logistic Regressions

Penalized logistic regressions predicting habitat suitability for sub-adult L. griseus and H. sciurus were constructed in R using the glmnet (v4.1.1) [77] and caret (v6.0.86) [78] packages, in R version 4.0.2. [64]. For each species, two logistic regressions were fit via penalized maximum likelihood. The first set of models applied the lasso penalty (i.e., L1-regularization, alpha = 1), a method that minimizes the absolute magnitude of the regression coefficients [79]. The lasso penalty reduces variance by shrinking or assigning a value of zero to some coefficients, thereby finding the optimal balance between model fit and complexity. The second set of models applied the ridge penalty (i.e., L2-regularization, alpha = 0), a method that minimizes the sum of the squared coefficients [80]. Unlike lasso, the coefficients in ridge regression can only asymptotically approach a value of zero. We applied these penalty terms independently to determine whether predictive performance varies based on regularization strategy.

The appropriate shrinkage parameters (lambda) for lasso and ridge were determined based on 10-fold cross validation using the cv.glmnet function of the glmnet package [77]. Penalized logistic regressions were then fit for each species separately using caret, again using 10-fold cross validation for model calibration. Thus, each model was fit using the 70% of presence–absence records set aside for model calibration and the 12 spatially explicit environmental predictors. Finally, the fitted training models were used to extrapolate predictions across the study area via the predict function of the R package raster [63].

2.4.4. MaxEnt Models

Maximum entropy models predicting relative habitat suitability for sub-adult L. griseus and H. sciurus were constructed using MaxEnt version 3.4.1. [81]. MaxEnt automatically applies L1-regularization to find the most parsimonious model [35,82]. The default regularization multiplier is 1.0, however, we used species-specific tuning to identify the regularization value and feature classes (i.e., functions of continuous environmental covariates) that enhanced predictive ability while minimizing overfitting [83,84]. After comparing five potential regularization multipliers (0.25, 0.5, 1.0, 2.0, and 5.0) using the ENMEval R package (v0.3.1) [85], it was determined that a value of 5.0, in combination with linear, quadratic, and hinge features, was optimal for constructing presence-background HSMs for sub-adult L. griseus and H. sciurus. Additionally, to prevent environmental bias stemming from spatially biased occurrence data [86,87], we created a Gaussian kernel density surface to capture the distribution of RVC and MVS sampling effort. This bias grid was fed into MaxEnt via the “bias file” option, enabling the sampling distribution to be factored out during construction of the training algorithm.

MaxEnt models were constructed using the 12 spatially explicit environmental predictors and the 70% of presence-only records designated for model calibration. Initial tuning and subsequent modeling were conducted using 10-fold cross validation and a set of 10,000 background points selected according to the bias file described above. The complementary log–log (cloglog) transformation was used, producing surfaces that reflected the predicted relative habitat suitability (or relative likelihood of occurrence) on a scale of 0 to 1.

2.5. Model Assessment

2.5.1. Discriminatory Ability

Predictive performance was compared between penalized logistic regression and MaxEnt models using the area under the receiver-operator curve (AUC). AUC is a threshold-independent, rank-based statistic that indicates a model’s ability to discriminate between a random absence (or background) point and a random presence point [35]. By assessing model performance over a variety of thresholds, the AUC test statistic provides an indication of discriminatory ability on a continuous scale and enables comparisons between modeling algorithms. AUC values range from 0 to 1, with the latter representing perfect discrimination. AUCs ≤ 0.5 suggest random or worse than random performance [88].

2.5.2. Binary Predictive Performance

We assessed the binary classification accuracy of each model using the 30% of presence–absence records set aside during the initial train–test split. The predicted suitability surfaces for each species were first discretized to a binary scale using a standard threshold of 0.5 for the predicted probability (penalized regression) or relative likelihood (MaxEnt) of presence. A confusion matrix was then calculated for each species–model combination, and the accuracy, sensitivity (i.e., percentage of correctly predicted presences), and specificity (i.e., percentage of correctly predicted absences) were examined. Though this standard threshold enables a general comparison of map accuracy between modeling algorithms, it may not represent the suitability cut-off at which the models optimally distinguish between suitable and unsuitable locations, which is an essential goal when defining nodes for subsequent connectivity modeling. Thus, we also identified the suitability thresholds at which each model achieved a maximum sum of sensitivity plus specificity (max SSS). Max SSS provides an indication of how conservative of a suitability cut-off must be used to maximize discrimination between the presence and absence (or suitable and unsuitable) locations [89].

2.5.3. Variable Importance

Variable importance scores were also calculated for each species across the three modeling techniques. For penalized logistic regressions built using glmnet and caret, variable importance was assessed using varImp, a function of the caret package that scales variables from 0 to 100 according to the absolute value of their standardized coefficients. For MaxEnt models, jackknife resampling was used to assess the influence of each predictor, this procedure sums the change in regularized gain (i.e., a goodness of fit measure based on a variable’s ability to distinguish species presence sites) across the ten cross validation folds. Regardless of model type, larger values indicate a higher level of importance.

3. Results

3.1. Discriminatory Ability

We used the AUC test statistic to determine whether penalized logistic regression and MaxEnt modeling techniques differ in their ability to discriminate between suitable and unsuitable locations across a variety of thresholds. According to the AUC statistic, regularization strategy had little effect on the discriminatory ability of regression HSMs, as lasso- and ridge-penalized models displayed similar performance within species. For sub-adult L. griseus, both penalized regressions achieved an AUC value of 0.74, indicating a good model fit. Discriminatory ability improved slightly for penalized regressions of sub-adult H. sciurus suitability, with lasso and ridge regressions producing AUCs of 0.76 and 0.75, respectively. However, regardless of the regularization strategy, the penalized logistic models were outperformed by MaxEnt, which yielded AUC values of 0.88 for L. griseus and 0.86 for H. sciurus (

Table 3. Performance assessment for the various species–model combinations. Discriminatory ability on a continuous scale was assessed using the area under the receiver–operator curve (AUC). Binary predictive performance was assessed using confusion matrices following discretization at two suitability thresholds—a standard threshold of 0.5 and the threshold at which each model achieved a maximum sum of sensitivity and specificity (max SSS). Accuracy, sensitivity, and specificity values are displayed as percentages and thresholds represent predicted suitability levels on a scale of 0 to 1.

	Discrimination	Binary Performance (Standard)				Binary Performance (Max SSS)
	AUC	Cut-Off	Accuracy	Sensitivity	Specificity	Cut-Off	Accuracy	Sensitivity	Specificity
Lutjanus griseus
Lasso Regression	0.74	0.5	77.9	32.7	93.3	0.28	69.4	73.9	67.8
Ridge Regression	0.74	0.5	77.7	32.7	93.1	0.24	65.5	80.0	60.5
MaxEnt	0.88	0.5	50.5	84.9	38.7	0.65	59.0	80.0	51.8
Haemulon sciurus
Lasso Regression	0.76	0.5	73.2	37.5	91.2	0.32	72.8	71.3	73.5
Ridge Regression	0.75	0.5	74.2	35.7	93.5	0.32	72.8	70.8	73.7
MaxEnt	0.86	0.5	41.3	79.1	24.3	0.65	51.7	72.6	42.3

).

3.2. Binary Predictive Performance

To produce the discrete patches of suitable habitat (i.e., nodes) required for modeling potential connectivity, the continuous habitat suitability surfaces must be discretized to a binary scale. Therefore, we used confusion matrices to assess each model’s binary predictive performance following discretization at two thresholds. When first discretized using a standard suitability cut-off of 0.5 and compared to the withheld validation data, lasso and ridge HSMs for sub-adult L. griseus achieved classification accuracies of 77.9% and 77.7%, respectively. Map accuracy was slightly lower for sub-adult H. sciurus, with values of 73.2% and 74.2% for lasso and ridge regressions, respectively. Relative to the penalized regressions, overall map accuracy for MaxEnt was low when assessed at this suitability cut-off, ranging from roughly 40–50% (Table 3). The models also varied in terms of sensitivity and specificity, with penalized regressions successfully identifying known absence (i.e., unsuitable) locations more frequently than known presence (i.e., suitable) locations and MaxEnt following the opposite trend.

The max SSS threshold selection strategy, which optimizes discrimination between known presence and absence localities, identified substantially different suitability cut-offs for each of the modeling algorithms (Table 3). In general, the penalized logistic models would have to lower their suitability thresholds to roughly 0.30 to maximize discrimination, whereas the MaxEnt models optimized discrimination between known presence and absence sites at a threshold of 0.65, a far less conservative cut-off. Although the overall map accuracy of the MaxEnt models still trailed those of the penalized regressions, the percentage of correctly predicted presence locations (i.e., sensitivity) was higher for MaxEnt than nearly all regressions. Additionally, discretization using max SSS thresholds increased the balance between sensitivity and specificity across all species–model combinations.

3.3. Variable Importance

To better understand the species–seascape interactions driving patterns of habitat suitability and, consequently, potential connectivity, we examined the continuous habitat suitability surfaces and variable importance plots produced by each model. Within species, there was a high level of agreement between the habitat suitability predictions of the three modeling techniques; however, lasso- and ridge-penalized regressions predicted smooth, gradual patterns of decreasing suitability as distance from shore increased, whereas MaxEnt captured patchy distributions of suitable habitat with noticeable fine-scale differences (Supplementary Figure S3). Furthermore, the models revealed species-specific responses to the various spatial predictors, resulting in considerable inter-species variation in predicted suitability across the seascape (Figure 3). For brevity, we focus here only on the variable contributions of the top-performing continuous HSM for each species (i.e., the MaxEnt models); variable importance information for the penalized regressions is provided in Supplementary Figure S4 and response curves for the top five MaxEnt predictors are provided in Supplementary Figure S5.

For sub-adult L. griseus, MaxEnt’s jackknife procedure identified benthic habitat type as the single most influential predictor of habitat suitability, followed by Euclidean distance to the nearest mangrove, slope, depth, and broad-scale BPI (Figure 4, Supplementary Figure S5). Predicted habitat suitability for L. griseus was high in shallow waters (<5 m depth) within roughly 200 m from mangroves, with the highest values predicted at the interface of shoreline mangroves and seagrass meadows (i.e., the seagrass fringe). Patches of high suitability were also identified along the shoreward side of the barrier reef tract, primarily in areas of increasing slope (>5 degrees) and over aggregate reefs and isolated patch reefs with shallow peaks surrounded by dense seagrass (Figure 3). In contrast, unconsolidated sediment and discontinuous seagrass were predicted to have the lowest relative likelihood of sub-adult L. griseus presence. Broad-scale BPI positively influenced predicted suitability levels, with broad peaks and ridges favored over flats and valleys. Though their roles were negligible relative to the top five predictors, mean winter salinity and mean summer temperature both showed negative relationships with predicted suitability when used in isolation. The remaining predictors—curvature, rugosity, fine-scale BPI, mean summer salinity, and mean winter temperature—were assigned a variable importance score of zero.

The five most influential predictors regulating the distribution of suitable habitat for sub-adult H. sciurus were benthic habitat type, slope, Euclidean distance to the nearest mangrove, depth, and broad-scale BPI (Figure 4, Supplementary Figure S5). Similar to the L. griseus model, predicted habitat suitability for sub-adult H. sciurus was high along the seagrass fringe and over individual patch reefs and aggregate coral reefs, especially those with broad peaks (Figure 3). However, relative to the L. griseus model, MaxEnt predicted higher suitability levels for H. sciurus over ridges and in patches of pavement, scattered rock, and reef rubble, and lower suitability levels in seagrass with the exception of those in shallow areas immediately adjacent to mangroves. Although depths shallower than 5 m were predicted to be the most suitable, the relative likelihood of H. sciurus presence remained above 50% at depths down to roughly 25 m. Predicted habitat suitability for H. sciurus also declined rapidly as Euclidean distance from the nearest mangrove approached 200 m, however, suitability levels for this species began gradually increasing again at a distance of around 1 km rather than continuing to decline. Furthermore, the mangrove shorelines on the windward and leeward sides of the Florida Keys had consistently higher estimates of habitat suitability than those along the mainland in Biscayne Bay (Figure 3). Although this trend was also visible in the L. griseus model, it was much more pronounced for H. sciurus. Water conditions contributed only weakly to overall model fit, however, there were positive relationships between the relative likelihood of sub-adult H. sciurus presence and mean winter salinity and temperature. In contrast, curvature, rugosity, fine-scale BPI, and mean summer salinity and temperature had contribution scores of zero.

4. Discussion

4.1. Model Performance

We compared penalized logistic regression and MaxEnt models of habitat suitability for two economically and ecologically critical reef fish species in the Florida Keys, USA, with the goal of identifying which modeling algorithm produces the most realistic and detailed products for use in subsequent connectivity modeling. MaxEnt’s AUC values were consistently higher than those of either the lasso- or ridge-penalized logistic regressions, suggesting that MaxEnt was better able to distinguish between suitable and unsuitable locations for sub-adult L. griseus and H. sciurus when evaluated across a range of suitability thresholds. Although the overall accuracy of the MaxEnt models fell below those of the penalized regressions when discretized using standard and max SSS suitability thresholds, MaxEnt produced similar or improved sensitivity estimates relative to the penalized regressions. MaxEnt’s consistently high sensitivity values suggest that these models were able to reliably identify known presence locations, which is essential for delineating the suitable habitat nodes required to produce spatial graphs for connectivity assessment. The high-resolution suitability maps produced by MaxEnt also appear to be better suited for conversion to resistance surfaces, as these models predicted patchy distributions of suitable habitat across the seascape that more closely aligned with the known ecology of the study species, as described below in Section 4.2.

Consistent with our findings, MaxEnt’s predictive performance has paralleled or exceeded that of other machine learning and regression techniques in several comparative studies. For instance, MaxEnt and penalized logistic regression techniques yielded similar AUC values when used to model the distributions of several tree species in Spain, with both models outperforming standard logistic regressions [90]. Similarly, MaxEnt achieved the highest predictive performance out of five modeling methods when used to model habitat suitability for the invasive Argentine ant across the Iberian Peninsula, producing predictive distributional maps that highlighted areas susceptible to invasion [39]. MaxEnt has even performed well in spatially and topographically complex seascapes, as demonstrated by a comparative study of ten presence-only modeling algorithms applied to demersal fish species of Australia [91]. Our research, therefore, joins these and other examples from the literature that illustrate MaxEnt’s usefulness as a tool for mapping species distributions across a range of taxa and environmental settings, especially in scenarios where distributional patterns are thought to be driven by complex species–environment relationships. These results can likely be attributed to MaxEnt’s ability to harness categorical data and linear, quadratic, hinge, and threshold functions of continuous environmental variables, while simultaneously maintaining a balance between model fit and complexity using regularization [82]. Despite its growing promise and popularity among ecologists, attempts to leverage MaxEnt products for connectivity modeling remain scarce.

4.2. Variable Importance

Of the 12 environmental predictors assessed in our study, habitat type was identified as the main driver of suitability for sub-adult L. griseus and H. sciurus, with dense seagrass beds, shoreline mangroves, patch reefs, and shoreward aggregate reefs being especially important. According to our models, these habitats play a variable role in supporting sub-adult L. griseus and H. sciurus, depending on their geographic location. This trend was especially apparent for shoreline mangroves, with suitability predicted to be highest in mangroves along the leeward and windward sides of the Keys and surrounding Biscayne Bay’s southern islands. Previous research in the region revealed a similar pattern of habitat use whereby sub-adult L. griseus and H. sciurus selected for easterly mangroves along the Keys, whereas larger-bodied adults were more common in Biscayne Bay’s expansive mainland forests [51,52]. The high suitability levels predicted along the seagrass fringe and in areas within several hundred meters of a nearby mangrove shoreline likely reflect the regular diel migrations of grunts and snappers between daytime resting sites in mangroves and nocturnal foraging grounds in seagrass meadows, which cover distances of up to 1 km [51,92,93]. Furthermore, patch reefs and sections of the reef tract adjacent to lush seagrass beds and mangrove shorelines, particularly those along the Keys, were predicted to have a higher relative likelihood of presence for both species compared to those located alongside unvegetated substrates. Together, these results suggest that patch reefs serve as critical stepping-stones connecting mangrove and seagrass nurseries to adult habitat on offshore reefs. However, sub-adult dispersal between back reef habitats may be limited to only immediately accessible patches of topographically complex habitat, with implications for the replenishment of adult populations across the barrier reef tract [94,95]. These findings also suggest that nutrient supplementation and top-down control of coral predators by migratory reef fishes may be greatest in patch reef, back reef, and reef crest communities neighboring seagrasses and mangroves, with consequences for coral growth and the placement of coral restoration efforts [44,48].

Beyond the type and spatial arrangement of benthic habitats, suitability varied as a function of depth and seafloor surface morphology, as quantified by slope and bathymetric position index (BPI). Generally, the likelihood of L. griseus and H. sciurus presence decreased with increasing depth, with this pattern being especially abrupt for L. griseus. The predicted suitability maps also displayed within-patch variation, with both species responding positively to even small increases in slope (<5 degrees). Additionally, our models revealed that the influence of BPI on habitat suitability is both scale- and species-dependent, with broad-scale BPI having a stronger influence over predicted suitability levels for L. griseus than H. sciurus. Nonetheless, fine-scale BPI and two other fine-scale metrics of seafloor surface complexity—curvature and rugosity—were dropped entirely from both MaxEnt models. These results indicate that the distributions of sub-adult L. griseus and H. sciurus across the seascape are primarily driven by broad-scale habitat features and geographic location, a finding that is consistent with previous research on predatory reef fishes and which may be related to the high vagility and large home range sizes of these species (1–5 km) [96,97,98,99]. However, previous work in the Caribbean revealed that seafloor morphology and geographic location interact to drive the distributional patterns of herbivores, in addition to invertivores and piscivores, suggesting that this trend is likely common across coastal reef fish communities rather than being restricted to mobile predators [68].

Although neither salinity nor temperature acted as a major determinant of suitability, these variables acted as filters to mediate the relative suitability of otherwise similar habitats. Previous research has demonstrated significant inter-species variation in salinity tolerances, with L. griseus being abundant in low-to-intermediate salinities and H. sciurus being abundant in stable, high salinities [53,100]. These relationships were reflected in MaxEnt’s predictions and were especially noticeable in Biscayne Bay, where salinity fluctuates significantly as a result of both freshwater discharge and tidal exchange. As such, the western and southwestern mainland coasts, which are characterized by extreme salinity fluctuations and lower overall means [59], were predicted to have higher suitability levels for L. griseus than H. sciurus. In contrast, habitats along the leeward and windward sides of the Upper Keys that have narrow salinity ranges dominated by seawater were predicted to have higher suitability levels for H. sciurus. Additionally, the present finding that relative suitability for sub-adult L. griseus decreases rapidly as summer temperatures approach 31 °C is in agreement with previous laboratory experiments that estimated 33 °C as being close to the maximum for juvenile gray snapper feeding [101]. Considering that settlement and grow-out occur from summer through early fall [51], this temperature constraint may be an artifact of juvenile habitat selection. On the other hand, the positive relationship between mean winter temperature and predicted suitability for sub-adult H. sciurus may reflect the winter spawning and settlement behavior of this species, as warmer winter temperatures decrease the overwinter mortality of juveniles and increase the chances of successful sub-adult dispersal to the offshore reefs [102,103].

4.3. Implications for Seascape Connectivity Modeling and Conservation

Based on our case study, habitat composition and arrangement, depth, and broad-scale bathymetric features are among the most important factors to consider when planning conservation efforts for reef fish species with complex, multi-habitat life histories. In particular, our models highlight the value of targeting mosaics of interconnected habitats, rather than single biotopes, when planning marine protected areas, reserve networks, and resource management [13,14,104]. The importance of considering surrounding seascape context and connectivity is not limited to the protection and conservation of existing ecosystems, but also to the restoration of degraded and fragmented habitats. By strategically placing restoration activities within the seascape, restoration practitioners can enhance intra- and inter-habitat connectivity, thereby increasing dispersal and harnessing key ecological processes including herbivory and primary production, predation and secondary production, and nutrient turnover [44]. Although the application of seascape connectivity as a spatially explicit metric in restoration planning has been limited thus far, projects that have incorporated surrounding seascape context in their site selection process have seen positive outcomes [105]. Integrative approaches that combine the strengths of habitat suitability modeling and connectivity modeling are becoming increasingly common and accessible thanks to improvements in data availability and the development of decision support tools like Marxan Connect and Zonation [15,24].

In conclusion, our research demonstrates that MaxEnt, a presence-background machine learning approach, outperforms traditional presence–absence techniques in terms of predictive performance and ability to produce habitat suitability maps that reflect the known ecology of sub-adult L. griseus and H. sciurus in the Florida Keys. These results are consistent with previous terrestrial studies that have found similar or improved predictive performance of MaxEnt relative to standard or penalized logistic regressions [37,90,106]. Furthermore, our research, coupled with previous work on warm-water kelps [107], demersal fishes [91], and shallow and deep-sea corals [108,109], demonstrate that MaxEnt’s promise as a fast, open-source tool for mapping species distributions extends beyond the land–sea interface. Although analogs exist from the terrestrial literature (e.g., [26,30]), the application of habitat suitability modeling as a precursor to marine connectivity assessments and conservation planning remains scarce. We anticipate that MaxEnt-derived nodes and linkages can be combined with data on species biological traits to produce detailed and ecologically realistic spatial graphs for connectivity assessment. As a next step, we plan to construct and operationalize these graphs for habitat restoration planning through scenario testing, including the iterative addition or removal of nodes and linkages (i.e., restoration and fragmentation scenarios, respectively). As habitat restoration efforts ramp up across the spatially and topographically complex seascape of the Florida Keys, we encourage restoration practitioners and coastal managers to adopt a multi-step site selection strategy that harnesses the strengths of both habitat suitability and connectivity modeling.