This study compiled a database of cold-water corals, focusing on both Hexocorallia and Octocorellia at predominantly the species level, with the data obtained from public databases, including the ICES Vulnerable Marine Ecosystems data portal, the NOAA Deep Sea Coral Data Portal, and the Ocean Biogeographic Information System portal (OBIS), as well as peer-reviewed scientific outputs. To ensure the accuracy of the presence records used for subsequent modeling, the data underwent three steps of filtering. Firstly, records with position accuracy > 5 km were removed. Secondly, records without position accuracy or depth information were removed. Thirdly, records with directly reported depths that differed by more than 200 m from the inferred depth based on the spatial position were removed. Finally, geographic filtering was applied to keep only one occurrence record in a single grid cell (25 km²) to reduce the influence of sampling bias.

The selection of background data for modeling in this study used commonly employed prevalence ratios of 1:1, 1:5, and 1:10, as well as a fixed number of 10,000 points (Phillips et al., 2006; Barbet-Massin et al., 2012; Hysen et al., 2022). Two methods were used to select background points: the target-group background sampling method and kernel density sampling method (Elith et al., 2010; Fitzpatrick et al., 2013; Cerasoli et al., 2017; Georgian et al., 2019; Burgos et al., 2020; Finucci et al., 2021; Georgian et al., 2021; Robinson et al., 2021; Stephenson et al., 2021; Anderson et al., 2022). To select target-group backgrounds for each species, random points were selected from the remaining cold-water coral presences with the selected points within 5 km of occurrences excluded. For kernel density backgrounds, sampling effort was modeled by fitting a kernel density estimate to presence locations of the six species studied (3001 locations remained with only one presence in each cell). The created two-dimensional estimated kernel density was used as a probability grid to select background points for each species, with points within 5 km of occurrences of the species excluded from selection using ArcGIS software.

The SRTM15+ V2.0 bathymetric grid with a resolution of 15 arc-sec (Tozer et al., 2019) was projected into the WGS_1984_EASE_Grid_2.0_Global coordinate system with a cell size of 5 km. The projected SRTM15+ V2.0 was further utilized to calculate terrain variables, including slope, curvature, and BPI9 (bathymetry position index with an analysis window size of 9×9) using ArcGIS software (Wilson et al., 2007) ( Table 1 ).

In this study, nine global environmental variables of the bottom layer, such as temperature and salinity, were chosen from the Bio-ORACLE v2.0 dataset with a grid size of 5′ (Assis et al., 2018) ( Table 1 ). These variables were calculated by utilizing pre-processed global ocean re-analyses, which combined satellite and in situ observations at regular two- and three-dimensional spatial grids with monthly averages over the period of 2000-2014 (Assis et al., 2018). Additionally, particulate organic carbon (POC) of the bottom layer was adopted as a candidate predictor with a spatial range (-75°31’30″S, 83°58″N) (after Davies and Guinotte, 2011). Furthermore, the variables aragonite saturation state (omega aragonite), calcite saturation state (omega calcite), and alkalinity, with a geographic range of approximately (-78.5°S, 68°N), were interpolated into seabed grids using three-dimensional gridded datasets GLODAP and WOA09, utilizing the cookie-cutter upscaling method first described by Davies and Guinotte (2011) and implemented by Yesson et al. (2017). All 13 environmental variables were resampled using the bilinear algorithm and projected into WGS_1984_EASE_Grid_2.0_Global with a cell size of 5 km and a global range of (-90S, 90N, -180W, 180E).

To investigate the potential influence of the collinearity of the environmental variables on predictive performance, the variance inflation factor (VIF) and Pearson’s correlation coefficient were utilized in this study. The VIF was used to estimate the correlation between each variable and all the other variables iteratively. Each time, the variable with the highest VIF was excluded until no VIF was larger than a set threshold (Yesson et al., 2015; Yesson et al., 2017; Khosravifard et al., 2020). The thresholds of VIF and Pearson’s r for this study were set as 5 and 0.8, respectively, with all VIFs <5 and Pearson’s r <0.8 indicating a low collinearity of the variables (Heiberger and Holland, 2004; Yesson et al., 2017). The R package HH was utilized for VIF calculation (Heiberger, 2022).

Four modeling algorithms, including MAXENT, MARS, BRT, and RF, were tested in this study. MAXENT predicts species’ potential distribution based on the maximum entropy principle, constrained by the features that the expected value of each environmental variable matches its empirical average (Phillips et al., 2006; Elith et al., 2011). BRT uses an iterative boosting algorithm to call the regression tree algorithm and build a combination of trees (Ridgeway, 1999) while stepwise adding modifications to the modeled regression trees to fit the data better (Friedman, 2001). RF uses the bagging method (bootstrap aggregation) to randomly select a number of bootstrap samples from the training data to fit the trees, then average all the fitted trees for a prediction (Breiman, 2001). RF is less sensitive to the tuning of model parameters than other classification and regression tree methods, such as BRT (Freeman et al., 2016). MARS is a flexible regression method that is used for fitting non-linear relationships or interactions, applying a piecewise linear function instead of smoothing (Friedman, 2001; Elith et al., 2006).

In this study, predictions were made using the R package Biomod2 4.1-2 (Thuiller, 2003; Thuiller et al., 2009; Thuiller et al., 2022). To achieve better model performance, an equal weight between presences and background points was set in this study, with the weighted sum of presences equaling the weighted sum of background points, as recommended by previous studies (Barbet-Massin et al., 2012; Liu et al., 2019). The default settings of Biomod2 4.1-2 were implemented for the additional arguments in all four models used in the study.

In this study, model evaluation was conducted by segmenting the data based on spatial blocks rather than resorting to random subsampling of presence and background data. The latter approach tends to lead to an overestimation of model performance, whereas the former is considered to yield a more robust estimate of accuracy (Santini et al., 2021). To achieve this, the presence records of each species were divided longitudinally into eight numerically equal folds, with the corresponding background points divided accordingly. For each replicate, one of the eight folds was withheld from the model calibration and reserved for evaluation (Valavi et al., 2019; Winship et al., 2020).

The accuracy of model predictions was evaluated using four metrics: sensitivity, specificity, AUC, and TSS. Sensitivity and specificity were computed as the proportion of correctly predicted as presences and absences, respectively (Barbet-Massin et al., 2012). Maximizing the sum of sensitivity and specificity has demonstrated promise as a method for selecting thresholds in presence-only predictions (Liu et al., 2013; Liu et al., 2015). In this study, the sensitivity and specificity yielding the highest combined value were utilized for assessing model performance.

The AUC statistic, a widely employed nonparametric method for evaluating species distribution predictions, was employed. AUC is insensitive to prevalence and is calculated as the area under the receiver operating characteristic curve (DeLong et al., 1988). It measures the model’s true positive rate (sensitivity) against the false positive rate (1 - specificity), with values ranging from 0.5 (indicating performance no better than random) to 1 (indicating perfect discriminatory capacity) (Jiménez-Valverde, 2012). Another commonly used metric, the TSS (true skill statistic) was used to gauge the agreement between predicted habitat suitability and the known presences/backgrounds of the validation dataset. TSS is calculated as the sum of sensitivity and specificity minus 1 (or true negative rate minus false negative rate) and ranges from -1 (indicating very poor performance) to 1 (indicating perfect agreement) (Allouche et al., 2006).

Pearson’s correlation was employed to assess the consistency of predicted habitat suitability across each run of algorithm pairs for each species (Grimmett et al., 2020). The Fleiss’ kappa statistic (Fleiss, 1971) was used to measure the agreement among predictions from 8 runs of each algorithm for each species, providing insights into the stability of predictions for each algorithm. This statistic measures agreement among multiple raters, whereas Cohen’s kappa estimates agreement between two raters (Grimmett et al., 2020).

Using M. oculata as a case study, further investigations were conducted into the smoothed spatial residuals (Renner et al., 2015) and the uncertainty of predicted habitat suitability. Spatial residuals were calculated using the following formula: observed species presence/backgrounds (1/0) - normalized predicted habitat suitability. The results were then interpolated into a grid using Ordinary Kriging. To assess the uncertainty of predicted habitat suitability, the standard deviation of predictions from eight runs was computed. This measure offers insights into the spatial sensitivity of the model to the sampling of occurrences and backgrounds in block cross-validation modeling.

The response patterns of TSS are generally consistent with AUC ( Figures 2A, B ). All predictions performed well, with mean AUC values > 0.75. Models using kernel density backgrounds have relatively stable higher mean values of AUC and TSS statistics than corresponding models using target-group backgrounds, for all algorithms, all presence prevalence (the ratio of presences to backgrounds), and both environmental datasets ( Figures 2 , 3 ). BRT and RF have better discriminative performance than the other two algorithms, with relatively stable high mean AUC values for almost all presence prevalence and both environmental datasets ( Figures 2 , 3 ).

The predictive performances of BRT, MARS, and RF using kernel density backgrounds fluctuate, with mean AUC values decreasing slightly as background prevalence decreases, whereas the mean AUC values of MAXENT increase as the prevalence decreases ( Figures 2A , 3A ). Models based on target-group backgrounds performed poorly compared to kernel density backgrounds, with more variability in mean AUC values ( Figure 3 ), and lower mean discriminative performance of MAXENT than BRT, MARS, and RF, particularly for higher prevalence ( Figures 2 , 3 ).

Runs of all four algorithms are more stable using kernel density backgrounds than target-group backgrounds ( Figure 3 ). Prediction outputs of all four algorithms are most stable for the geographically widespread M. oculata, followed by E. rostrata and D. pertusum, using kernel density backgrounds, with narrow ranges of AUC values; meanwhile, the predictions of G. dumosa consistently show the worst performance across all algorithms and backgrounds ( Figure 4 ).

The mean sensitivity and specificity of predictions using kernel density backgrounds are both higher than predictions using target-group backgrounds. The mean sensitivity and specificity of RF, BRT, and MARS are mostly higher than that of MAXENT, regardless of species prevalence and environmental datasets used. Additionally, the mean sensitivity of each algorithm is higher than the corresponding specificity ( Figures 2C, D ).

Greater agreement is shown between predictions of each pair of algorithms using kernel density backgrounds than using target-group backgrounds with both environmental datasets among all four cases of species prevalence ( Figure 5 ). Similar characteristics were found among the predictions for the six species ( Figure 6 ). The highest agreement is evident between BRT and RF predictions using kernel density backgrounds, exhibiting a mean Pearson’s correlation coefficient ranging from 0.834 to 0.938. This is followed by BRT and MARS, with Pearson’s correlation coefficient ranging from 0.816 to 0.872, and then MARS and RF ( Figure 5 ). In contrast, MAXENT predictions exhibit the lowest level of agreement with predictions using the other modeling techniques ( Figure 5 ).

The greatest mean agreement between BRT and RF was observed in predictions using a prevalence ratio of 1:1 and the filtered environmental dataset (12 environmental variables), yielding a Pearson’s correlation coefficient of 0.938. This is followed by predictions using a prevalence ratio of 1:1 and the whole environmental dataset (17 variables), resulting in a Pearson’s correlation coefficient of 0.925. The level of agreement decreases as prevalence decreases ( Figure 5 ).

When using kernel density backgrounds, the mean agreement between BRT and RF predictions ranks among the top for five out of six species, except for G. dumosa. Madrepora oculata and E. rostrata show the highest mean agreement in predictions for each pair of algorithms, while the least agreement is found in G. dumosa predictions for each pair of algorithms ( Figure 6 ).

Fleiss’ kappa statistics highlight a stronger mean agreement of predictions for each algorithm when using kernel density backgrounds compared to target-group backgrounds. The mean Fleiss’ kappa ranges between 0.46 and 0.641 in predictions using kernel density backgrounds ( Figure 7 ).

The strongest agreement of predictions for each algorithm is achieved in cases where the filtered environmental dataset is used, and an equal number of background points to presences is used (prevalence 1:1) for BRT, MARS, and RF predictions. For MAXENT predictions, a prevalence ratio of 1:10 is employed. The highest mean Fliess’ kappa values are recorded for RF (0.641, 12KD1), followed by BRT (0.613, 12KD1), MARS (0.596, 12KD1), and MAXENT (0.606, 12KD10) ( Figure 7 ).

The lowest agreement for each algorithm is observed in predictions of G. dumosa, using both kernel density backgrounds and target-group backgrounds, as well as both environmental datasets. For M. oculata, high Fleiss’ kappa values are noted in BRT, MARS, and RF predictions, particularly when using the filtered environmental dataset (12KDMO), with corresponding values of 0.729, 0.756, and 0.689, respectively. Conversely, the Fleiss’ kappa value for E. rostrata is high in MAXENT prediction, particularly when using the filtered environmental dataset (12KDER) at 0.685. Notably, high agreement for each algorithm is also observed in predictions for S. varibilis using target-group backgrounds ( Figure 8 ).

Figure 9 ; Supplementary Figure 1 illustrate that the spatial residuals of RF predictions are lower compared to predictions of BRT, MARS, and MAXENT. When using a 1:1 ratio of presence to backgrounds, BRT predictions exhibit fewer areas of high residuals compared to using backgrounds of 1:5, 1:10, and 10K. A similar trend is observed in MARS predictions. However, in the case of MAXENT predictions, using a 1:1 ratio of presence to pseudo-absence leads to more areas of high residuals compared to backgrounds of 1:5, 1:10, and 10K.

Figure 10 ; Supplementary Figure 2 illustrate that the uncertainty of BRT and RF predictions is lower compared to MARS and MAXENT predictions. Among them, MAXENT predictions display the highest level of uncertainty, with the 1:10 ratio prediction showing fewer areas of high uncertainty compared to other ratio cases.

Supplementary Figure 3 provides visual assessments of the predicted habitat suitability maps for M. oculata. It reveals a similar spatial pattern among the predicted habitat suitability of the four modeling algorithms, with high habitat suitability generally aligning with the spatial locations of M. oculata presences.

The study revealed a pronounced spatial bias in species presence records, with the highest sampling intensity concentrated in the North Atlantic and SW Pacific regions ( Figure 1 ). Spatial sampling biases in biodiversity data emerge from intricate interactions among geography, species traits, and human behavior (e.g., preferences or avoidance of certain species or habitats) (Baker et al., 2022). Addressing this sampling bias poses a significant challenge for presence-only and presence-background species distribution models. Biased data can obscure the actual correlation between occurrences and environmental predictors, irrespective of the model used, leading to predictive distributions that predominantly reflect sampling effort rather than actual habitat suitability (Baker et al., 2022; Barber et al., 2022). In this study, the sampling strategy of background points was identified as the most influential factor on predictive performance, surpassing the effects of modeling algorithms, species prevalence (the ratio of species presences to background points), or collinearity of environmental datasets ( Figure 2 ).

Kernel density background data demonstrated significantly superior performance compared to target-group backgrounds in achieving higher discriminative accuracy and functional performance across all four algorithms, all species prevalence cases, both environmental datasets, and different niche breadths ( Figures 2 , 7 , 8 ). Although direct performance accuracy comparisons between predictions using kernel density backgrounds and target-group backgrounds have been sparse in previous studies, these methods have been favorably compared to other techniques aimed at reducing sampling bias. Target-group backgrounds outperformed radius-restricted backgrounds and covariates of sampling effort (e.g., maps of human population density and road networks) used as alternative means to estimate sampling effort, effectively mitigating the impact of sampling bias (Barber et al., 2022). Moreover, modeling with target-group backgrounds exhibited superior performance compared to random backgrounds (Phillips et al., 2009) and environmental-filtering backgrounds (Iturbide et al., 2015). Modeling using kernel density backgrounds exhibited a greater improvement in recreating known distributions compared to presences-filtering methods, which encompass spatially thinning presences to a fixed sampling density using geographic filtering or maintaining equal sampling intensity across environmental space using environmental filtering (Inman et al., 2021).

We speculate that the comparable distribution pattern of the six species with target-group backgrounds, often found on the continental shelf margin, seamounts, and steep slopes of sea islands ( Figure 1 ), could heighten the probability that background points (cells) indeed contain species occurrences. This increased likelihood could lead to greater difficulty in distinguishing between presences and backgrounds, ultimately diminishing the predictive performance of models using target-group backgrounds.

The study revealed that different algorithms performed optimally with distinct background sampling sizes, a relationship also observed in previous studies (Barbet-Massin et al., 2012; Liu et al., 2019). The findings regarding machine learning methods BRT and RF are consistent with previous research. Grimmett et al. (2020); Barbet-Massin et al. (2012); Jimenex-Valverde (2021); Hysen et al. (2022), and Barker & MacIsaac (2022) all support an equal prevalence of presences and absences. Liu et al. suggested a small multiplier in background points relative to the number of presences (Barbet-Massin et al., 2012; Liu et al., 2019; Grimmett et al., 2020; Jimenez-Valverde, 2021; Barker and MacIsaac, 2022; Hysen et al., 2022).

The conclusions regarding the optimal size of background samples for the MARS algorithm have been conflicting in previous studies. Hysen et al. (2022) found that MARS performed optimally with 10,000 background points for a sample size of 248 (Hysen et al., 2022), while Barbet-Massin et al. (2012) suggested fewer pseudo-absences (e.g., 100) for sample sizes of 30, 100, 300, or 1000 presences. However, our findings align with Barker and MacIsaac (2022), who recommended equal random pseudo-absences to centroids (presences). The reduced areas of high residuals in the predicted habitat suitability of MARS using a 1:1 ratio of presence to backgrounds, compared to other background ratios ( Supplementary Figures 1, 2 ), also signify the higher accuracy of predictions using a 1:1 ratio.

Likewise, conclusions regarding the optimal size of background samples for the MAXENT algorithm have been conflicting in previous studies. Liu et al. (2019) found that MAXENT performed optimally using background points as a small multiple of presences, and Barker and MacIsaac (2022) suggested equal random pseudo-absences to centroids (presences). However, our findings align with Grimmett et al. (2020) (e.g., 1900 background points for 100 presences), Phillips and Dudík (2008) (e.g., 10,000 for a few to thousands of presences) (Phillips and Dudík, 2008), and Hysen et al. (2022) (prevalence 1:10 is better than prevalence 1:1 and 10,000 sample size). The reduced areas of high residuals and lower uncertainty in the predicted habitat suitability of MAXENT using a 1:10 ratio of presence to backgrounds, compared to other background ratios, also indicate the higher accuracy of predictions using a 1:10 ratio.

Concerning prediction consistency, our results correspond with Grimmett et al. (2020), who identified the strongest agreement (Fleiss’ kappa statistic) between RF predictions using equal presence backgrounds and between MAXENT predictions using a large number of background points. Moreover, our results align with other studies that advocate equal presence-backgrounds in BRT, RF, and MARS predictions, while recommending the use of a substantial number of background points for MAXENT predictions.

No single modeling algorithm outperformed all others under all circumstances. Nonetheless, in this study, BRT & RF performed best overall using equal prevalence. However, MAXENT and BRT exhibited better performance than RF in predictions involving 225 species with dozens to thousands of presences (Elith et al., 2006; Valavi et al., 2021). Similarly, MAXENT outperformed other algorithms, such as a variant of RF, BRT, and MARS, in predictions involving 171 species in both random and spatial partitioning (Valavi et al., 2023). However, a number of previous studies have reported different outcomes, which are consistent with our results. For instance, Barker and MacIsaac (2022) found BRT to outperform RF, MARS, and MAXENT in predictions involving virtual species; Romero-Sanchez et al. (2022) found RF to achieve better performance than MAXENT in predictions for plus trees; and Hysen et al. (2022) found RF to perform better than MARS and MAXENT in nest predictions (Barker and MacIsaac, 2022; Hysen et al., 2022; Romero-Sanchez et al., 2022). Additionally, the smaller number of areas with high residuals and lower uncertainty in the predicted habitat suitability of RF and BRT, compared to those of MARS and MAXENT in M. oculata prediction, also indicates the higher accuracy of predictions using RF and BRT.

Our finding that the greatest agreement of MAXENT showed lower than BRT or RF is in conflict with the result of Grimmett et al. (2020), which found that MAXENT performed more consistently in terms of discriminative accuracy and spatial prediction stability than RF using 20 to 100 presences, and a background sample size determined by presences/(presences + background points) = 0.5, 0.1, 0.05, 0.01, and 0.005 (Grimmett et al., 2020). Additionally, Grimmett et al. (2020) found that Pearson’s correlation coefficient between RF and MAXENT increased with prevalence (~0.3-0.9). The inconsistent result from our study could be related to different sample sizes (20-100 presences in their study, 316-1223 presences in this study) and validation metrics between studies (random cross-validation in their study and spatial block cross-validation in this study). In general, increasing sample size has a positive effect on performance (measured by AUC and Schoener´s D), with model accuracy decreasing with a sample size below 300 presences (Gábor et al., 2020). Both Fleiss’ kappa of MAXENT and RF predictions increased with increasing sample size in Grimmett et al. (2020), which may indicate that the conflicting results may be related to differences in sample size.

The AUC values were found to be similar between predictions using the 12 filtered (uncorrelated) environmental variables and the full set of 17 correlated variables for each of the four modeling techniques, particularly when using kernel density backgrounds, indicating that the collinearity of the environmental dataset did not influence the discriminative accuracy of the four modeling methods significantly. Dormann et al. (2013) also found that the modeling methods significantly, GAM, MARS, BRT, and RF, worked reasonably well under moderate collinearity (Dormann et al., 2013). de Marco Junior and Nobrega (2018) found that the intensity of the effect of collinearity varied according to the algorithm characteristics, with more complex models, such as MAXENT, performing better than simple envelope ones based on PCA-derived variables (de Marco Junior and Nobrega, 2018). However, the collinearity of the environmental dataset did influence the functional performance of algorithms, with stronger agreement of each algorithm found in predictions using the VIF-filtered environmental dataset compared to when using the whole dataset. Furthermore, spatial autocorrelation in the predictors can inflate the variable importance estimates when the response of species to the environmental gradients is linear (Harisena et al., 2021). Therefore, it is recommended to address the collinearity between environmental variables before modeling.

It is worth noting that the methodologies employed in our study, namely, the use of different modeling algorithms, background sampling strategies, number of background points, and collinearity of environmental variables, are not exclusive to the study of cold-water corals, although specific cold-water coral species within the Hexocorallia and Octocorellia groups were used as a case study. These approaches can be widely applicable in species distribution modeling across various marine benthic taxa. Therefore, the insights gained from our investigation into the impact of these factors on predictive performance could potentially extend to other benthic groups and marine species. However, different species with different habitat preferences, ecological requirements, and dispersal abilities may exhibit unique responses to these factors. Therefore, while the specific results of our study pertain to cold-water corals, the underlying principles and considerations can serve as a foundation for researchers working with other marine species.