Ecosystem simulation models are valuable tools for understanding the impacts of environmental and anthropogenic stressors on marine ecosystems and for informing resource management (Fulton, 2010; Christensen and Walters, 2011; O'Farrell et al., 2017). A wide diversity of ecosystem modeling platforms is now available, including the spatially-explicit modeling frameworks Atlantis (Fulton et al., 2004, 2007, 2011), Ecopath with Ecosim with Ecospace (Pauly et al., 2000; Christensen and Walters, 2004; Walters et al., 2010) and OSMOSE (Shin and Cury, 2001, 2004; Grüss et al., 2016c). The simulation of predator-prey interactions is a key feature of ecosystem models (Plaganyi, 2007; Christensen and Walters, 2011; Grüss et al., 2016a). To simulate predator-prey interactions, spatially-explicit ecosystem models need to represent patterns of spatial overlap between predators and prey. To be able to represent patterns of spatial overlap between predators and prey, many spatially-explicit ecosystem modeling platforms, including Atlantis and OSMOSE, need to be provided with distribution maps defining the annual or seasonal spatial distributions of functional groups¹ and life stages.

The extrapolation of spatial distribution patterns predicted by statistical models integrating environmental covariates (commonly called “species distribution models”) is a useful means to produce distribution maps for spatially-explicit ecosystem models (Drexler and Ainsworth, 2013; Grüss et al., 2014, 2016d; Hattab et al., 2014). For example, Drexler and Ainsworth (2013) fitted generalized additive models (GAMs) to monitoring data collected in the U.S. Gulf of Mexico (GOM). The authors then extrapolated the predictions made by their GAMs to the entire GOM Large Marine Ecosystem (LME) to obtain distribution maps for some of the functional groups (all life stages combined) represented in the Atlantis model of the GOM LME (“Atlantis-GOM”; Figure 1). While Drexler and Ainsworth's (2013) methodology is more sophisticated than usual methodologies for generating distribution maps for ecosystem models, it has some limitations. Firstly, Drexler and Ainsworth (2013) employed only one groundfish trawl survey dataset to produce all their distribution maps, which resulted in unreliable predictions of spatial distributions for those functional groups that are poorly sampled by groundfish trawl (e.g., some of the pelagic functional groups represented in Atlantis-GOM). Secondly, Drexler and Ainsworth (2013) used the same set of environmental covariates (dominant sediment type, bottom depth, bottom temperature, bottom dissolved oxygen (DO) concentration, and surface chlorophyll-a concentration) in all their GAMs; however, other environmental parameters may have a stronger influence of the spatial distribution patterns of the functional groups represented in Atlantis-GOM (e.g., sea surface temperature (SST) in the case of most pelagic functional groups). Lastly, Drexler and Ainsworth's (2013) GAMs did not account for spatial autocorrelation (spatial structure), because their purpose was to be employed for performing extrapolation over the entire GOM LME; this led to significant, unmodeled spatial patterns in Drexler and Ainsworth's (2013) GAM residuals for functional groups associated with small-scale habitat features, such as red grouper (Epinephelus morio) and gag (Mycteroperca microlepis) (unpublished data).

Interpolation from the predictions of geostatistical models is a means to account for spatial structure in spatial distribution patterns (Guisan and Zimmermann, 2000; Saul et al., 2013; Thorson et al., 2015; Grüss et al., 2017b). For instance, Grüss et al. (2017b) fitted geostatistical binomial generalized linear mixed models (GLMMs) to a blending of fisheries-dependent and fisheries-independent monitoring data collected in the U.S. GOM. The authors then interpolated geostatistical GLMM predictions to a spatial grid covering the entire U.S. GOM to be able to generate distribution maps for an OSMOSE model of the West Florida Shelf. Geostatistical binomial GLMMs accurately reflect the probability of encountering functional groups/life stages in a marine ecosystem. However, these models cannot be used to interpolate the spatial distributions of functional groups/life stages over the entire GOM LME, due to the scarcity or limited access to monitoring data collected in the southern GOM (i.e., the Mexican and Cuban GOM and the international waters of the GOM).

To address the above mentioned issues, we developed a novel methodology combining extrapolation and interpolation of the predictions made by statistical habitat models to produce distribution maps for the fish and invertebrates represented in the Atlantis-GOM ecosystem model. This methodology consists of fitting binomial GAMs to a large monitoring database for the U.S. GOM and a large environmental database for the GOM LME, and geostatistical binomial GLMMs to the large monitoring database for the U.S. GOM; and then employing GAM predictions to infer fish and invertebrate spatial distributions in the southern GOM, and geostatistical GLMM predictions to infer fish and invertebrate spatial distributions in the U.S. GOM. The large monitoring database we compiled for the U.S. GOM gathers all the fisheries-dependent and fisheries-independent monitoring data collected in the region between 2000 and 2016 using random sampling schemes. Our large environmental database gathers all the environmental parameters known to influence the spatial distributions of fish and invertebrates of the GOM that are available for both the U.S. and the southern GOM. Thus, our methodology allows for reasonable extrapolation in the southern GOM based on a large amount of monitoring and environmental data, and for interpolation in the U.S. GOM accurately reflecting the probability of encountering fish and invertebrates in that region. Though initially designed for the GOM, our methodology can serve other regions of the world characterized by a large surface area and only partially covered by available monitoring data, such as LMEs bordered by several countries (e.g., the Benguela Current LME).

We present our methodology, which allowed us to construct a total of 116 annual and seasonal distribution maps for the fish and invertebrate functional groups and life stages represented in the Atlantis-GOM ecosystem model. We first briefly present the GOM LME and the Atlantis-GOM model. Then, we provide an overview of our methodology. Next, we describe the compilation of the large monitoring database for the U.S GOM and of the large environmental database for the GOM LME. Afterwards, we explain how we developed and evaluated binomial GAMs for the entire GOM LME and geostatistical binomial GLMMs for the U.S. GOM. Finally, we combine GAM predictions and geostatistical GLMM predictions according to our methodology to generate distribution maps for the fish and invertebrates represented in Atlantis-GOM.

In the present study, we developed a novel methodology for producing distribution maps for a spatially-explicit ecosystem model using large monitoring and environmental databases and a combination of interpolation and extrapolation. This novel methodology has a number of advantages. First, it makes use of all the monitoring data collected in the ecosystem of interest using random sampling schemes, thereby allowing the development of more reliable statistical models. Then, our methodology performs extrapolation using GAMs integrating only environmental covariates that are pertinent for the functional groups and life stages of interest. Finally, our methodology performs interpolation using geostatistical models in the regions where monitoring data are available, thereby accounting for spatial structure in spatial distribution patterns in those regions. This is particularly desirable because the residuals from GAMs fitted to monitoring data typically show spatial structure (as was the case in the present study; Supplementary Figure 1).

In order to construct reliable distribution maps for the Atlantis-GOM ecosystem model, we decided not to strictly follow our methodology in some cases. First, when a GAM did not pass the validation test, we used a GAM for a related functional group/life stage to generate distribution maps for the GOM LME. Second, we did not fit geostatistical GLMMs for those functional groups and life stages whose depth, longitudinal and latitudinal ranges are not entirely covered by the monitoring data from the large monitoring database for the U.S. GOM. For those, only GAM predictions were employed to obtain distribution maps for Atlantis-GOM. From the large monitoring and environmental databases, we were able to produce 53 annual maps and 64 seasonal maps (for 32 different functional groups/life stages) for the Atlantis-GOM model. We had no data for filter feeding sharks [basking shark (Cetorhinus maximus), and whale shark (Rhincodon typus)] in the large monitoring database for the U.S. GOM. To be able to generate a distribution map for this functional group, we fitted a MaxEnt model (Phillips et al., 2006, 2017) to presence-only data for basking and whale sharks obtained from the Ocean Biogeographic Information System⁴ (Supplementary Data Sheet 3). Thus, ultimately, we produced a total of 117 annual and seasonal maps for the fish and invertebrate functional groups and life stages represented in the Atlantis-GOM ecosystem model (Supplementary Figure 4).

The usefulness of our methodology combining GAM and geostatistical GLMM predictions is exemplified by the distribution maps generated for pink shrimp and adult red grouper (Figure 4). The distribution maps produced for juvenile and adult pink shrimps from GAM and geostatistical GLMM predictions concurred with previous studies in that pink shrimp hotspots are found on the West Florida Shelf, around the Dry Tortugas and the Florida Keys, in the region off Alabama and Mississippi where depth ranges between 20 and 40 m, on the Texas-Mexico border, and on the Campeche Bank (Costello and Allen, 1970; Bielsa et al., 1983; Ramirez-Rodriguez et al., 2003; Drexler and Ainsworth, 2013; Grüss et al., 2014). However, reflecting monitoring data for the U.S. GOM, the geostatistical GLMM of adult pink shrimp predicted that the life stage has a very high probability to be encountered in the region off Alabama and Mississippi where depth ranges between 20 and 40 m, contrary to the GAM of adult pink shrimp. More importantly, the GAM of adult red grouper was unable to correctly predict the spatial distribution of the life stage in the U.S. GOM. Consistent with the literature (Coleman et al., 1996, 2011; Lombardi-Carlson et al., 2008; Albañez-Lucero and Arreguin-Sanchez, 2009), the GAM and the geostatistical GLMM of adult red grouper predicted that the life stage has a high probability to be encountered over the entire West Florida and on the Campeche Bank. However, like the GAM fitted in Drexler and Ainsworth (2013), the GAM of adult red grouper fitted in this study also wrongly predicted that the life stage has a high probability to be encountered over the rest of the U.S. GOM shelf, while there is a consensus in the literature that adult red grouper is very rarely encountered west of the Florida-Alabama border (Coleman et al., 1996, 2011; Franks, 2005; Burns and Robbins, 2006; Lombardi-Carlson et al., 2008); it has been suggested that the presence of red grouper west of the Florida-Alabama border relates to displacement of individuals by storms and hurricanes or anomalous oceanic current patterns resulting in the settlement of grouper larvae west of the De Soto Canyon (Franks, 2005; Burns and Robbins, 2006; Dance et al., 2011; Karnauskas et al., 2013).

A major limitation of the methodology introduced in this study is that it relies on GAMs that do not account for spatial structure. One can account for spatial structure in GAMs through the integration of a spatial autocorrelation term (Guisan and Thuiller, 2005; Zuur et al., 2014), or of an interaction term between eastings and northings (Wood, 2006; Grüss et al., 2016d); however, it was not possible to implement either of these options in the present study, since we needed to perform extrapolation to the entire GOM LME while we had monitoring data only for the U.S. GOM. The autocorrelated residual variance found for the great majority of the functional groups/life stages in this study (Supplementary Figure 1), combined with the broad spatial scale of the maps of probability of encounter we generated, mean that these maps are not particularly useful for localized predictions. By necessity, mapping the general spatial distribution of functional groups/life stages usually requires some abstraction of the populations. Hence, the predictions we made with GAMs in the present study would not be useful for informing whether one would be likely to capture functional groups/life stages at a particular location of the GOM LME. This is evident, for example, if we compare the ranges of autocorrelation between the residuals of the model fits for gag and red grouper life stages found in this study (Supplementary Figure 1) with the smaller ones found when modeling individual catch per unit effort information in Saul et al. (2013); the ranges of autocorrelation for gag and red grouper found in Saul et al. (2013) were much closer to the sizes of individual habitat patches. Thus, it is important to emphasize that the modeled products that we reported in this study are valuable to capture larger spatial patterns necessary for running simulations with ecosystem models, but not to predict the spatial distribution patterns of fish and invertebrate populations with high precision.

Another limitation of the present study is that, for each functional group, we relied on one reference species to determine whether to model ontogenetic differences in spatial distribution patterns and environmental influences on spatial distribution patterns (Supplementary Table 3). In other words, we made the assumption that all the species of a given functional group share the same spatial distribution patterns, namely those of the reference species of the functional group (although data extraction and statistical model fitting considered all the species making up the study functional groups). This assumption is reasonable given the large number of fish and invertebrates we had to deal with (63) and the very low spatial resolution of the Atlantis-GOM ecosystem model (Figure 1B). Yet, studies using our methodology to construct distribution maps for a spatially-explicit ecosystem model with the higher resolution than that of Atlantis-GOM may want to proceed another way; for example, these studies may want to rely on several reference species for deciding whether the spatial distribution patterns of the juveniles and adults of each functional group differ and which environmental parameters should be integrated in the GAM of each functional group.

Some modifications could be introduced in the structure and validation process of our geostatistical models in future studies. For simplicity, we did not integrate environmental parameters in our geostatistical GLMMs, because preliminary analyses conducted for previous studies using geostatistical models to infer fish and invertebrate spatial distributions in the U.S. GOM found that spatial residuals in probability of encounter have a much larger influence on probabilities of encounter than environmental variables (Grüss et al., 2017b,c). Also, typically, relationships between the probability of encounter of fish and invertebrates and environmental parameters are non-linear, which is something that the geostatistical GLMMs employed in this study cannot represent at present. Thus, it appeared more appropriate not to integrate environmental parameters in our geostatistical models, and we recommend that the future studies wishing to integrate environmental parameters in our geostatistical models should first incorporate non-linearities into these models. Moreover, to validate our geostatistical models, we used standard procedures, i.e., convergence diagnostics and analyses of Pearson residuals (Thorson et al., 2015; Grüss et al., 2017b,c). Future studies could implement the Leave Group Out Cross Validation procedure for our geostatistical models and estimate their adjusted R² and AUC; it would then be possible to compare the AUCs of binomial GAMs with those of binomial geostatistical GLMMs fitted to the same monitoring data and, thus, to compare the ability of binomial GAMs and binomial geostatistical GLMMs to discriminate between encounters and non-encounters.

We generated maps of relative abundance for the Atlantis-GOM ecosystem model from probability of encounter maps. Thus, as is often done in the mapping literature (e.g., Maxwell et al., 2009; Hattab et al., 2013; Cormon et al., 2014; Grüss et al., 2014, 2017b,c), we assumed a linear relationship between probability of encounter and abundance. This assumption is appropriate for many of the functional groups represented in Atlantis-GOM that are structure-oriented; the abundance of structure-oriented functional groups, such as red grouper and gag, may be primarily determined by the encounter/non-encounter of a suitable structure (Saul et al., 2013; Grüss et al., 2017c). Yet, spatial patterns of probability of encounter and spatial patterns of abundance often rely on slightly different determinants (Nielsen et al., 2005; Koubbi et al., 2006; Aarts et al., 2012; Grüss et al., 2014); thus, if we had produced distribution maps from abundance predictions (e.g., if we had had monitoring data for the entire GOM LME), some of those distribution maps may have been slightly different from the distribution maps constructed in the present study. In the future, if we had access to a reasonable amount of monitoring data for the southern GOM, it would be interesting to evaluate the consequences of employing abundance vs. probability of encounter maps to produce maps of relative abundance for Atlantis-GOM with our methodology.

The goal of this study was to obtain the maximum of data possible to generate distribution maps for a spatially-explicit ecosystem model. However, while our methodology represents a an improvement over previous ones (Grüss et al., 2016a), further improvements could be made. First, for the sake of statistical rigor, only monitoring data collected using random sampling schemes were employed to fit statistical models in the present study. Developing methods to use data from both random and fixed sampling stations would allow us to fit geostatistical models for coastal functional groups/life stages such as juvenile red drum and juvenile and adult snooks (Table 4) for which a large amount of monitoring data is collected using fixed sampling schemes in Louisiana, Mississippi and Alabama (Grüss et al., 2016a). Moreover, rather than combining interpolation and extrapolation as we did, future studies could employ only interpolation and combine the predictions made by different geostatistical models, each fitted to the monitoring data collected in a given region (e.g., in the waters of the different countries bordering a LME). Such an endeavor would be highly challenging in our case, due to the scarcity or limited access to monitoring data collected in the southern GOM.

Our inability to fit geostatistical binomial GLMMs for some of the functional groups and life stages represented in Atlantis-GOM revealed data gaps for two groups of fish and invertebrates of the U.S. GOM: (1) coastal functional groups and life stages, such as juvenile red drum, juvenile and adult snooks, ladyfish and oysters; and (2) large pelagic fish, such as bluefin tuna and white marlin (Table 4). Data collected at random stations are missing for Louisiana, Mississippi and Alabama for the first group, as mentioned above. In parallel to the development of new statistical methodologies employing monitoring data collected using fixed sampling schemes, we encourage the initiation of new monitoring programs using random sampling schemes in Louisiana, Mississippi and Alabama, for the sake of having a set of monitoring programs using similar sampling protocols over the entire U.S. GOM. Data for species such as bluefin tuna and white marlin are collected almost exclusively by the National Marine Fisheries Service Pelagic Observer Program, which operates in the offshore waters of the GOM LME, while these species are found both inshore and offshore in the GOM. We recommend the use of fisheries-independent longline surveys to collect data for these species in inshore areas of the GOM and the development of statistical models that can be fit to a mix of data collected using fixed and random samplings schemes.

The methodology developed in the present study will allow for the use of more reliable inputs in spatially-explicit ecosystem models. Spatially-explicit ecosystem models, especially highly sophisticated models such as applications of the Atlantis modeling platforms, rely on numerous data inputs, many of which can have a significant impact on ecosystem model predictions (Fulton et al., 2007; Shin et al., 2010; Walters et al., 2010; Steele et al., 2013). Many spatially-explicit ecosystem models have been designed worldwide (Fulton, 2010; Espinoza-Tenorio et al., 2012; Steele et al., 2013), and significant attention has been devoted to enhance the structure and assumptions of ecosystem modeling platforms (e.g., Shin et al., 2010; Steenbeek et al., 2013, 2016; Travers-Trolet et al., 2014; Grüss et al., 2016b,c). Therefore, it is now high time to develop and employ methodologies such as ours for ensuring that the spatially-explicit ecosystem models available around the world are provided with reliable inputs generated using the best available data and information. We hope that the methodology developed in this study will serve other regions of the world, particularly LMEs bordered by several countries (e.g., the Benguela Current LME, and the Humboldt Current LME).

Designed and analyzed the models: AG, MD. Conceived the models: AG, MD, CA, and EB. Conducted the literature review: AG and JT. Queried the data: AG and ML. All authors contributed to writing the paper.