Spatio-temporal distribution of the spinetail devil ray Mobula mobular in the eastern tropical Atlantic Ocean

The distribution of the spinetail devil ray Mobula mobular in the eastern tropical Atlantic remains poorly known compared to the Pacific and Indian Oceans. We used fisherydependent data and generalized additive models to examine the environmental characteristics associated with the presence of M. mobular in the eastern Atlantic Ocean. Results revealed that the distribution of M. mobular is significantly associated with seasonal upwelling systems in coastal and pelagic areas. Our model predicted the presence of the species in areas where there is evidence of its occurrence, such as the Angolan upwelling system and the coast of Ghana. In addition, our model predicted new hotspot areas, including locations around the Mauritanian upwelling system, the Guinea coast, offshore Ghana and the south coast of Angola and Brazil, where sample sizes are limited. Those areas, as well as the environmental preferences depicted by the model, provide valuable information about the habitat and ecology of the spinetail devil ray. Future research lines derived from this study, as well as its limitations, are discussed. Furthermore, in light of our results we discuss the improvements that are needed to contribute to the conservation and management of this vulnerable species.


INTRODUCTION
Fishing activities have a direct effect on vulnerable large pelagic species such as elasmobranchs (Moyes et al. 2006). Mobulid rays are subject to fishing pressure, being caught either as target species in smallscale fisheries using a variety of gears (e.g. gillnets, driftnets, harpoons, gaffs, traps, trawls and longlines) or as bycatch in large-scale fisheries with longlines, trawls and purse seines (Ward-Paige et al. 2013, Croll et al. 2016). Due to their life-history characteristics (i.e. they can live more than 20 yr, e.g. Mobula birostris, but their reproductive cycle is likely to last 1 yr, giving birth to a single pup) (Couturier et al. 2012, Dulvy et al. 2014) and the increase in catches that have occurred in recent decades (from 200− 5000 t yr −1 between 1998 andWard-Paige et al. 2013), mobulid rays are considered highly vulnerable to overexploitation (Lawson et al. 2017) and have recently been added to Appendix II of the Convention on International Trade in Endangered Species (CITES) and to Appendices I and II of the Convention of Migratory Species (CMS) to regulate their international trade (Lawson et al. 2017). *Corresponding author: nlezamaochoa@gmail.com

Spatio-temporal distribution of the spinetail devil ray Mobula mobular in the eastern tropical Atlantic Ocean
Information on mobulid biology, population structure, ecology, distribution and habitat is crucial to the design of effective management measures to conserve their populations. In the specific case of manta rays (M. alfredi, M. birostris), there is evidence of a tendency to aggregate at specific locations to feed or clean, and this has led to the design of spatial management measures which aim to protect them (Croll et al. 2016, Barr & Abelson 2019. However, for datapoor mobulid species, specifically for mobulids that inhabit pelagic waters, information about their distribution, seasonal patterns and movements are still scarce, which makes it difficult to develop effective conservation measures. The tropical tuna purse-seine fishery is among the fisheries with the highest number of bycaught mobulid rays in the eastern Pacific Ocean, with an annual average of 2664 individuals (as per observations between 1993 and 2014) (Hall & Roman 2013, Lezama-Ochoa et al. 2019a. Although there is no information on population abundance of mobulid species in the eastern Pacific Ocean, these bycatch numbers are considered significant, given the lack of evidence of their survival when released back into the sea. The new handling and release techniques promoted by the European tropical tuna purse-seine fishery and the development of the Code of Good Practices (CGP) (Poisson et al. 2014, Goñi et al. 2015, Grande et al. 2020 to reduce mortality of vulnerable species (including mobulids) set the foundations by the Inter-American Tropical Tuna Commission (IATTC) to adopt Resolution C-15-04, promoting the correct management of mobulid ray species.
In the eastern Atlantic Ocean, tropical tuna purseseine fisheries mainly capture 4 species of mobulid rays: M. mobular, M. tarapacana, M. thurstoni and M. birostris. Some of these species primarily inhabit oceanic waters, whereas others are more restricted to coastal and shelf waters (Ward-Paige et al. 2013). The spinetail devil ray (Müller & Henle 1841) M. mobular is a circumglobal pelagic species that is found in tropical and subtropical waters and can inhabit both shallow coastal and offshore waters. Until recently, the spinetail devil ray was thought to consist of 2 separate species (M. mobular and M. japanica), one being endemic to the Mediterranean Sea and the other having a circumglobal distribution (Couturier et al. 2012). However, a recent taxonomic revision (White et al. 2018) reclassified it as a single species, M. mobular. M. mobular is listed as Endangered globally by the IUCN Red List of Threatened Species (Marshall et al. 2019).
Few studies have provided information about the distribution of M. mobular in the Atlantic Ocean, and most of these have been based on very localized visual (aerial or underwater) censuses in areas of the Mediterranean Sea (Notarbartolo di Sciara et al. 2015) and the Azores (Sobral 2013, Sobral & Afonso 2014. With regard to the eastern tropical Atlantic, the distribution of mobulids in general (and M. mobular in particular) has been poorly documented compared to the Pacific and Indian oceans (Weir et al. 2012). The first sightings of M. mobular in this region were recorded by local fishermen and documented by Cadenat (1959) in Abidjan (Ivory Coast). Some studies have reported mobulid rays as bycatch in longliners (Mas et al. 2015) and purse seiners (Ariz & Gaertner 1999, Amandè et al. 2010, Torres-Irineo et al. 2014, Grande et al. 2020, including M. mobular catches in most of them. It is difficult to study such highly mobile and oceanic species due to lack of easy access to the animals and their aggregation sites, making it even more important to get information to fill the knowledge gaps that still exist about their ecology, distribution and conservation status in the eastern Atlantic Ocean (Weir et al. 2012).
Compared to other devil ray species, few electronic tagging studies have been performed to study the habitat use of M. mobular (Canese et al. 2011, Croll et al. 2012, Francis & Jones 2017. For example, Croll et al. (2012) reported the species' preferences for both shallow (< 50 m depth) and warm waters (> 20°C) in the Gulf of California. Lezama-Ochoa et al. (2019b) revealed that the distribution of M. mobular is directly associated with productivity systems in the eastern Pacific Ocean, yet there are no studies describing these possible relationships in the Atlantic Ocean. Establishment of effective conservation measures that could protect the habitat -and hence the population -of mobulid rays depends on a clear understanding of the spatio-temporal distribution of the species and its relationship with environmental variables (Graham et al. 2012).
Species distribution models (SDMs) have been widely used in ecology to understand the relationships between the presence and/or abundance of a species and the environment (Whittaker et al. 2005, Araujo & Guisan 2006, Albouy et al. 2019. These models can be used to inform management options, design protected areas or forecast climate change effects for the development of future mitigation measures, among other uses (Graham et al. 2012). Regarding the tropical tuna purse-seine fishery, different SDMs have been developed to understand the distribution of target and non-target species (Martínez-Rincón et al. 2012, 2015, Sequeira et al. 2014, Lezama-Ochoa et al. 2016, Montero et al. 2016, Druon et al. 2017, with generalized additive models (GAMs) being the most used tool for modeling the distribution and environmental preferences of large marine species in the eastern Atlantic Ocean (Escalle et al. 2015, Lezama-Ochoa et al. 2018). However, these methods have not been widely used to study mobulid rays, and specifically M. mobular, in the Atlantic Ocean.
Based on records from various European tropical tuna purse-seine observer programs, we developed GAMs to explore the relationship between environmental variables and M. mobular in the eastern Atlantic tuna purse-seine fishery. We hypothesized that the distribution of M. mobular is associated with seasonal productive systems in coastal and pelagic areas of the eastern Atlantic Ocean. This study provides a valuable species case study to develop statistical SDMs that will help understand the ecology of M. mobular and ultimately inform conservation and management of this Endangered species and other mobulid species in the future.

Study area
The study area was located in the eastern tropical Atlantic Ocean, between 40°W, 15°E and 20°N,20°S (Fig. S1 in Supplement 1 at www.int-res.com/ articles/ suppl/n043p447_supp/). This region is the main fishing ground of European (French and Spanish) tropical tuna purse-seine vessels operating in the Atlantic Ocean (Escalle et al. 2016). The area is influenced by different oceanographic processes, such as permanent (occurring during the whole year with different intensity, i.e. Mauritanian upwelling) and seasonal (occurring during a specific period of the year, i.e. equatorial upwelling) upwelling systems. For example, the Benguela and Canary Currents generate permanent and seasonal upwellings from South Africa to Gabon (July−September) and from Mauritania to Senegal (November−May) (Hardman-Mountford et al. 2003, Tomczak & Godfrey 2003. There is a seasonal equatorial upwelling system that originates from the trade winds, begins in spring and reaches maximum chlorophyll values during summer (Binet & Marchal 1993). Meanwhile, the Congo River provides freshwater input and dissolved organic matter to the Cape Lopez area, which en hances primary production during the wet season (April− September) (Hardman-Mountford et al. 2003). Other mesoscale processes, such as the seasonal Guinea (summer) and Angola (winter) thermal domes located around 12°N, 11°W and 10°S, 9°E, respectively, increase productivity in those waters at dif ferent times throughout the year (Binet & Marchal 1993).

Data collection and environmental variables
Mobula mobular bycatch data was collected by the European Union (EU) scientific observer program on tropical tuna purse-seine fisheries operating in the eastern Atlantic Ocean. Since 2003, Spain and France have coordinated an observer program as part of the Spanish and French National Programs for Fisheries Data Collection, established through European Regulations (Commission Regulation [EC] No. 665/2008) (Gondra et al. 2017). These programs are run by the French Institut de Recherche pour le Développement (IRD) and the Spanish Instituto Español de Oceanografía (IEO) and AZTI (Lezama-Ochoa et al. 2018).
In 2012, the observer coverage increased to 100% during January and February in order to monitor the fishing operations of all purse seiners targeting tropical tunas and supply vessels 1 during a 2 month fishing aggregating device (FAD) fishing closure in the area off western Africa (ICCAT 2016, Gondra et al. 2017). Moreover, the Spanish tuna purse seiner associations Asociación Nacional de Armadores de Buques Atuneros Congeladores (ANABAC) and Or ganización de Productores Asociados de Grandes Atuneros Congeladores (OPAGAC), operating in the Atlantic Ocean, established a voluntary agreement for the application of a CGP for responsible tuna fishing activities -including the requirement for 100% coverage -which has been applied since 2013. Observer coverage rate has varied over the years (Amandè et al. 2010, Gondra et al. 2017, with higher observer coverage since 2012, reaching up to 61% in 2015 (Table S1).
Data recorded by the observer programs include information about all fishing activities conducted during the trip: set day and hour, set type, set position and catch of target (tunas) and bycatch species (in number or biomass). Bycatch groups include billfish, sharks, bony fish, rays, turtles and mammals (Lezama-Ochoa et al. 2018). Mobulid rays used to be identified by the observer program simply as 'Mantas' or to family level (Mobulidae). Over the years, however, the number of 1 Supply vessel: vessels used to search for, assess, maintain, retrieve, monitor and deploy drifting FADs for the benefit of the purse seiners individuals identified to genus level (Mobula sp.) has increased progressively with the improvement of observer experience, and thanks to the CGP and the development of improved identification guides. Currently, species are identified to genus and species level (Poisson et al. 2014, Grande et al. 2020).
In the present study, a fishing set was considered to be the sample unit and, depending on the fishing strategy used, each set was categorized as freeswimming ag gregations (free school sets) or associated (FADs) (Lezama-Ochoa et al. 2018). Sets associated with whales (mysticetes) were considered free school sets. Natural (logs) and artificial (man-made) objects, as well as sets associated with whale sharks, were considered FADs (Amandè et al. 2011, Escalle et al. 2016, Lezama-Ochoa et al. 2018. A total of 719 fishing trips and 17 127 fishing sets were observed between 2005 and 2015, of which 10 576 (61.7%) were FAD sets and 6551 (38.2%) were free school sets ( Fig. 1). A total of 433 sets with the presence of M. mobular were observed during the study period, of which 343 were FAD sets (79.3%) and 90 were free school sets (20.7%). Most of the fishing effort (number of sets) was located in the equatorial area and in coastal locations of Mauritania, Ghana or Gabon (Fig. 1).
The database contained information about M. japanica, M. mobular and M. rancurelli. Currently, M. japanica and M. rancurelli are considered to be junior synonyms of M. mobular (White et al. 2018, G. Notarbartolo di Sciara pers. comm.), and therefore all Mobula individuals were considered to be M. mobular. A total of 1032 individuals of M. mobular were identified; 828 (80.23%) in FAD sets and 204 (19.77%) in free school sets. The spatial-temporal distribution of fishing sets with the presence of M. mobular showed that the species occurrence appears to be variable in space and time (see Fig. S1b−d).
All environmental variables were collected from the EU Copernicus Marine Environment Monitoring Service (CMEMS) (http:// marine. copernicus. eu). This platform aggregates information from different databases on the physical state, variability and dynamics of the ocean and marine ecosystems worldwide (Lezama-Ochoa et al. 2019b). The following environmental variables (for the years 2005−2015) were downloaded: daily oxygen concentration (O 2 ; mg l −1 ), daily nitrate (Ni; mg l −1 ), daily chlorophyll (chl; mg m −3 ), daily sea surface temperature (SST; in°C ), daily salinity (Sal; in PSU) and daily sea surface height (SSH; in m) ( Table 1). All variables were extracted at 1/4° spatial resolution and downloaded from the 'Ocean product' section of the Copernicus database (http://marine.copernicus.eu/services-portfolio/ access -to-products/). The Copernicus Service provides guides to help the user write Python scripts and download the data in 'NetCDF' format. Once the 'NetCDF' file was downloaded for each variable, it was processed in R (R Core Team 2019). For more details on the catch− environmental data matching process, see Lezama-Ochoa et al. (2019b), and for the detailed information for each environmental variable, see Table 1. Information on statical variables, such as depth and distance from the coast, were obtained from Global Marine Environmental Datasets (GMED) (http:// gmed. auckland.ac.nz/download.html). Specifically, bathymetry values were obtained from the General Bathymetric Chart of the Oceans (GEBCO, 08 Digital Atlas; British Oceanographic Data Centre; www.gebco.net) (depth; m; 30 arc-s resolution). Land distance, or distance to the nearest land cell (water cells only), was calculated using Eu clidean distance (distance; km × 100; 5 arc-min resolution). Both variables were obtained as rasters (ASCII format) and processed in R (R Core Team 2019) using the 'rasters' package (Hijmans & van Etten 2012). Each corresponding environmental value ex tracted from the position data of species bycatch was included in the final dataframe.

Regime shifts identification
Lopez et al. (2020) identified 3 different environmental regimes in the eastern tropical Atlantic Ocean that are separated by calendar days, Day 155, Day 240 and Day 290 and correspond to (1) a stable period and reduction of the summertime productivity warm season, from October (i.e. Day 290) to June (i.e. Day 155); (2) a cool season, characterized by high productivity, high dissolved oxygen concentrations and low SST values, from early June (i.e. Day 155) to late August (i.e. Day 240) and (3) a warm season, from late August (i.e. Day 240) to mid-October (i.e. Day 290), with lower productivity and warmer sea surface waters (for more details see Lopez et al. 2020). Given that our study focused on a species bycaught in the same fishery and in the same areas as the ones investigated by Lopez et al. (2020), the same environmental regimes proposed by Lopez et al. (2020) were considered in this study. These 3 environmental regimes (i.e. stable, cool and warm) were used here to characterize the distribution of M. mobular in different seasonal oceanographic regimes, representing different activities: (1) transition or migration between warm− cool season (stable season), (2) feeding activities (cool season) and (3) reproductive or breeding activities (warm season).

Modeling
Relationships between the explanatory variables (oceanographic variables) were examined to avoid correlation and collinearity and to reduce overfitting. Firstly, covariates were examined for correlation using pair plots and Pearson's rank correlations in R (R Core Team 2019). Pairs of variables with high correlation values (Pearson correlation r > 0.6) were identified and only one of the correlated pairs was included in the modeling process, based on the literature and the ecology of the species (Zuur et al. 2010, Passadore et al. 2015, Lezama-Ochoa et al. 2019b). Secondly, a variance inflation factor (VIF) analysis was conducted using a cut-off value of 3 (Zuur et al. 2009) to avoid collinearity. This was done by using the 'vifstep' function in the 'usdm' package in R (Naimi 2015). Based on these tests, distance and SST were removed from the analysis. All other covariates had low correlation and collinearity values. The correlation and collinearity analyses are shown in Fig. S2.
GAMs were used to model presence−absence of the species (response variable) in relation to environmental characteristics in the eastern Atlantic, assuming a binomial distribution and a logit link function. The general structure of the GAM is represented as: where g is the link function (logit for binomial family), μ i is the expected response variable (presence− absence), α is the intercept, f n are smooth functions (thin plate regression splines), and X n are the covariates (Guisan et al. 2002).
Other possible modeling algorithms, such as zero-inflated negative binomial models, were considered in previous steps in this work due to the high number of zeros in several sets. However, mobulid rays are a data-poor bycatch species, and their abundance is difficult to model when attempting to identify specific hotspot areas, leading to difficulties in identifying strong relationships with the environment. A manual forward step-wise variable selection procedure was applied to select the final model covariates (Wood 2006). As an additional measure, to avoid model overfitting, the degrees of freedom of the smooth functions were restricted for each explanatory variable (measured as number of knots, k) to k = 4 for the main effects (Giannoulaki et al. 2013, Jones et al. 2014, Lezama-Ochoa et al. 2018b). The final model was chosen based on the lowest Akaike's information criterion (AIC) value, which represents the most parsimonious model (Akaike 1974).

Model validation and evaluation
Ideally, model performance is assessed using an independent data set. When this is not available, the performance of the model is usually validated by splitting the data into 2 sets: one to model the relationship between presence−absence data and the en vironmental variables, called the training data, and the other to evaluate the quality of the predictions, called the testing data. In this study, a k-fold cross-validation method (with k = 5) was used to construct the testing (20%) and training (80%) data (Elith & Leathwick 2009).
Model predictions were assessed through a confusion matrix using the 'PresenceAbsence' package in R (Freeman & Freeman 2012). The area under the receiver-operating curve (AUC) and the mean true skill statistic (TSS) indices were calculated from the confusion matrix (Pearson 2010). The AUC provides an overall measure of accuracy that is thresholdindependent and ranges from 0 to 1. Values around 0.5 indicate that the prediction is as good as random and values close to 1 indicate perfect prediction (Fielding & Bell 1997). This index measures the ability of the model to correctly predict where a species is present or absent (Elith et al. 2006). From the confusion matrix we also obtained 'specificity' (proportion of absences correctly predicted), and 'sensitivity' (proportion of observed occurrences correctly predicted).
The TSS is another accuracy index but it is threshold-dependent and not affected by the size of the validation set (Allouche et al. 2006). The TSS index, which is calculated as sensitivity plus specificity minus 1, ranges from −1 to +1, where 0 indicates no predictive skill, +1 indicates perfect agreement and values of zero or less indicate a performance no better than random (Brodie et al. 2015). Usually both performance scores are used in combination to evaluate the predictive performance of a distribution model.

Environmental preferences of M. mobular
Daily spatial prediction maps of the probability of M. mobular presence and their average and standard deviations at a 1 × 1° resolution for the whole period (2005−2015) and for each previously identified regime (i.e. stable, cool and warming seasons) were ob tained from the final GAM model using the 'predict.gam' function of the 'mgcv' package (Wood 2014). All statistical analyses were performed using R free software (R Core Team 2019).

RESULTS
The final model, which included depth, chlorophyll, SSH and oxygen concentration, explained 10.2% of the total deviance, with an AIC value of 4480.383 (the most parsimonious model based on the lowest AIC value) (Tables 2 and S2). The individual contribution of each variable is shown in Table 2. Depth (5.69%), chlorophyll (3.85%) and SSH (3.85%) were the most significant covariables to predict the presence of Mobula mobular. Relationships between the response variable and the explanatory covariates are shown by smooth functions in Fig. 2. The final explanatory variables included in the model suggest that species presence is associated with productive systems, particularly in coastal areas of < 2000 m depth. M. mobular is more likely to be found in areas with relatively high concentrations of chlorophyll (0.2−1 mg m −3 ), low concentrations of oxygen (<220 mmol l −1 ) and negative SSH values (−0.2 m) (Fig. 2).
AUC, sensitivity, specificity and TSS values for cross-validated models are shown in Table 3 to identify the most important areas for the species despite the low number of presences in the data set. The model predicted different areas with higher species presence probability in the eastern tropical Atlantic Ocean, particularly around Cape Verde, Cameroon and Angola (around 0.03−0.04) and along the coast of Senegal, Ghana and Angola (around 0.02) (Fig. 3). Average predictions by environmental regime showed a preference of M. mobular for temporary and permanent productive systems throughout the year. During the stable season (October− June), M. mobular presence probability peaked off Guinea and the Angola coast, whereas in the cold season (June−August) the probability of presence increased along the coasts of Cape Verde, Gabon and Angola, the Cape Lopez area and the eastern coast of Brazil (Fig. 4, Animation S1 in Supplement 2 at www.int-res.com/ articles/ suppl/ n043p447 _ supp/). The southwestern portion of the prediction map is an area of seamounts crossed by the Mid-Atlantic ridge that could benefit the aggregation of mobulid rays by providing opportunities for cleaning or feeding, and the environmental conditions of the region could be most beneficial for the presence of the species (Stevens et al. 2018). During the warm season (August− October), however, the probability of M. mobular presence decreased in general but remained in some parts of the Mauritania, Cameroon and Angola coasts (Fig. 4, Animation S1). Daily predictions provided more accurate information by months (Animation S1).

DISCUSSION
Modeling the distribution of datapoor marine species is challenging, particularly for large, mobile species with low observer coverage that have potentially been significantly impacted by fisheries (Temple et al. 2018). This study explored Mobula mobular habitat preferences in the eastern tropical Atlantic Ocean using bycatch data from the tropical tuna purse-seine fishery. SDMs based on observer bycatch data can provide new information on the distribution of this Endangered species.
Model outputs suggest a relationship between M. mobular occurrence and the most productive coastal waters of the eastern tropical Atlantic Ocean. Areas of mixed waters usually concentrate plankton and generate areas of high productivity, where filterfeeding large pelagic species, such as mobulid rays, tend to aggregate (Luiz et al. 2009, Stewart et al. 2017, Barr & Abelson 2019. Food availability is only one of several factors that could drive the distributions of elasmobranchs; for example, species are known to make seasonal migrations related to sexual reproduction and parturition. Two of the 4 most influential variables in the final model were depth and SSH. Depth explains the species' preference for coastal waters (areas < 2000 m depth), whereas SSH is related to mesoscale processes that are usually linked to productivity (Kahru et al. 2007). Oxygen con centration was also an important variable in the SDM. In the present study, the variables selected for , which suggests that this species can be found in areas with the same environmental variables across oceans but with different environmental preferences in each case (different chlorophyll and/or temperature values). The Atlantic and Pacific oceans share similarities in terms of oceanographic characteristics, such as the seasonal divergence that occurs along the equator as a consequence of the trade winds that blow from the northeast and southeast, which create upwelling and bring denser, nutrient-rich water up from the bottom of the ocean (Picaut et al. 1984). Based on these similarities, we consider that environmental preferences of M. mobular may be similar across regions.
The distribution of M. mobular seems to change seasonally in concordance with the seasonal oceanographic process of the eastern Atlantic Ocean, likely related to food availability. Chlo ro phyll shows the highest concentrations during specific times of the year and in places where mobulid rays ag gregate for feeding. These increases in productivity are due to upwelling systems related to these oceanographic processes. M. mobular may show a preference for permanent and seasonal mesoscale processes such as eddies (which retain productive bodies of water) or domes (Jaine et al. 2014). The mean of the predictions shows a preference for coastal areas. The predictions by regime and day, however, provide more details than quarterly analyses about the species' daily and seasonal distribution throughout the year (see Animation S1). For example, M. mobular may be associated with the permanent upwelling system of the Mauritania and Angola coasts all year long, but with different in tensity levels depending on the season.
Average predictions by regime and day suggest that the species had some association with the Guinea and Angola domes during winter and summer, respectively. These domes are meso scale structures which concentrate high productivity habitat. They are characterized by upward displacement of iso therms in the thermocline layers down to depths of more than 300 m (Siedler et al. 1992). These small cyclonic gyres (in the Northern Hemisphere) promote upwelling due to the net influence of Ekman transport, which drives colder, nutrient-rich water up into surface water (i.e. promote mixing; Wyrtki 1981). The Guinea Dome shows similar oceanographic characteristics to the eastern Pacific Costa Rica Dome (Wyrtki 1981), where M. mobular appears to aggregate during the summer (Lezama-Ochoa et al. 2019a,b).
In the Equatorial area during May−September, the presence of M. mobular can be partially explained by both the equatorial upwelling and the thermal elevation structure close to Cape Lopez (1.5°S, 8°W−3°E). The upwelling in those areas is associated with wind-   Table 3. Accuracy indices to evaluate the performance of the model for each interaction and the mean. AUC: Area Under the Curve; TSS: True Skill Statistic induced divergence of the surface waters and is considered to be different than the Guinea and Angola Domes (Wyrtki 1981). Specifically, the productivity system that originates around the Cape Lopez area during the beginning of the cold season is the result of both warm and less saline waters in the north that are separated from cold, more saline waters coming from the south (Lopez et al. 2020). In addition, during the spring and summer, discharge from the Congo River can increase nutrients in the water and, therefore, the productivity of these waters during this time of year. The coastal upwelling along the Ghanaian coast during summer and around the Cameroon coast during winter could also be associated with wind induced-upwelling processes, which explains the presence of the species (Nieto & Mélin 2017). The probability of presence of M. mobular in the southwestern region of the eastern Atlantic Ocean may be explained by the presence of seamounts and the mid-Atlantic ridge. These elevated areas are a perfect habitat for this species, since they can be used for protection, reproduction, cleaning or breeding (Stevens et al. 2018). Finally, we did not have any observations for the Brazilian coast, but the model predicted the possible presence of M. mobular in this area throughout the whole year, with highest predictions during late summer. Up -welling events along the Brazilian coast originate during spring and summer (Luiz et al. 2009, Coelho-Souza et al. 2012. Relationships between oceanographic processes and the distribution of pelagic/semi-pelagic elasmobranchs and bycaught species have been previously studied in the eastern tropical Atlantic Ocean (Escalle et al. 2016, Lezama-Ochoa et al. 2016, Gaube et al. 2018, Braun et al. 2019. Lopez et al. (2020) found relationships between silky sharks, mesoscale features and seasonal upwelling events in the eastern tropical Atlantic Ocean. In that sense, remote sensing data offers the opportunity to study in detail the oceanographic processes that explain the distribution of large marine species (Hazen et al. 2013). However, to date no studies have an alyzed environmental drivers of M. mobular distribution in the Atlantic Ocean, or used fisheries-dependent data or alternative independent data sets. Moreover, only a few studies (Sobral 2013, Thorrold et al. 2014) have tried to relate environmental conditions to occurrences of M. mobular, and none of these studies were in the tropical Atlantic Ocean. Information on M. mobular in the eastern Atlantic Ocean is scarce. Future telemetry, stable isotope and genetic studies are crucial to validate current results and models but also to finally shed some light on the basic biology and ecology of these poorly known Endangered species, which will be the only way to actively work towards sound management and conservation measures. The hotspot areas predicted by the SDM developed in this study could be important in helping researchers to implement tagging programs or small-scale initiatives, increasing the probability of successful tagging missions. Data obtained through those tagging initiatives could also be used to validate current models and develop further SDMs. As the main areas of importance are coastal regions, local management measures for conservation (economic and social) should be taken into account in future management plans. These measures will likely not be sufficient to ensure the conservation of the species, but they will mark the beginning for future panels of discussion about correct management procedures, for example, regarding local small-scale fisheries. Collaborations between stakeholders, scientists, non-governmental organizations (NGOs) and industry to monitor and improve bycatch mitigation measures have increased in recent years. Observer programs have become necessary, and are valuable as they are sometimes the only source of information available to study vulnerable species (Gondra et al. 2017).
One of the main problems with mobulid ray bycatch data is the lack of correct species identification, since members of the family share various morphological similarities, thus making it difficult to obtain reliable information at the species level (Ward-Paige et al. 2013). Moreover, some species are caught by multiple fisheries, which highlights the need to analyze different fisheries data to assess the overall impact on their populations (Croll et al. 2016). Ideally, genetic studies should also be included in future projects so that population structure can be studied; these types of studies will ultimately help with the correct identification of mobulid species (Couturier et al. 2012, White et al. 2018, Hosegood et al. 2019. Small-scale artisanal fisheries also capture mobulid rays along the African coast (Diop & Dossa 2011), however this information is limited and usually not published. It would be important for future studies to include this information in order to improve, compare and develop more SDM models and predictions.

Limitations of the study and future perspectives
The SDM developed here was able to explain only a limited percentage of deviance. Model prediction depends entirely on the quality and quantity of the data, the processes to be modeled and the availability of environmental variables included in the model. Although the results of our study are in line with previous resources and literature (Cadenat 1959, Couturier et al. 2012, Weir et al. 2012, Kaschner et al. 2013, Lezama-Ochoa et al. 2019b, they should be interpreted with caution due to the low level of explained deviance. However, these results can still be useful and, particularly for data-poor bycatch species, can provide new information about species distribution as here for M. mobular (Scales et al. 2017). When bycatch species data are modeled, there is usually only low deviance explained by the model . In contrast to target species, such as tunas, where 1 or 2 environmental factors (i.e. temperature) may explain their spatial distribution, the presence of bycatch species may be explained by diverse factors not considered, known or included in the models. Despite this data gap, the information obtained from the observer programs is useful to improve our knowledge about the most vulnerable bycatch species, with the aim of reducing their mortality and protecting them, as here for M. mobular. Thus, based on the different calculated statistics indices (AUC, sensitivity, specificity, etc.), we can consider the performance of the model as good, based on the available data and previous knowledge about this species in the studied area (which can be confirmed and improved).
Mobulid ray distributions all overlap with the main fishing areas in free school sets in very productive waters in the Pacific tropical area (Lezama-Ochoa et al. 2019a). However, in the present study most of the data came from FAD fisheries, which dominate the eastern Atlantic purse-seine fishery. Although mobulids do not exhibit strong aggregation behavior around floating objects (Lezama-Ochoa et al. 2019a,b), the model was able to predict M. mobular presence along the equatorial area, where a greater distribution of FADs is usually found, and where the equatorial upwelling system could be sufficiently productive to provide habitat for this species.
Some recent studies have analyzed the post-release survival of elasmobranchs in purse seiners and longlines, with most of the effort focused on sharks (Musyl et al. 2011, Hutchinson et al. 2015, Escalle et al. 2018, Musyl & Gilman 2018. In the case of mobulid rays, Francis & Jones (2017) reported the post-release survival of 7 out of 9 individuals of M. mobular that were tagged in the tuna purse-seine fishery in New Zealand waters. In the eastern Pacific Ocean, a pilot project is currently in progress with the objective of estimating the post-release survival rates of mobulids using similar handling and releasing practices as in purse seine vessels (J. D. Stewart, pers. comm.). Future studies in the eastern Atlantic Ocean should also be directed toward the study of post-release survival of mobulid rays to contribute to the identification of the most vulnerable species and to improve and/or modify the existing methods in the purseseine fleet.
There is no stock assessment of M. mobular in the Atlantic Ocean. As a data-poor bycatch species, there is no information available to develop robust stock assessment models to support effective management measures. Alternative methods to quantify the impacts of purse-seine fisheries and prioritize studies that focus on species of potential concern are needed. The use of ecological risk assessments (ERA) and specifically, productivity−susceptibility analyses (PSA), is an alternative solution to assess the relative vulnerability of data-limited bycatch species such as mobulid rays (Griffiths et al. 2017). Recently, Griffiths et al. (2018) developed a flexible quantitative approach in the Eastern Pacific Ocean that uses fewer input parameters than PSA to quantify the cumulative impacts of multiple fisheries on M. mobular and uses predicted distribution maps provided by SDMs (Griffiths et al. 2019). The application of this new method based on SDMs, using fishery-dependent data from the eastern Atlantic Ocean, could provide a better assessment of the impact of the global purseseine fishery on M. mobular.

CONCLUSIONS
This study improves our understanding of the environmental characteristics associated with the presence of M. mobular in the tropical eastern Atlantic Ocean, based on purse-seine fishery observer data and GAMs. Results suggest that M. mobular has a preference for productive systems in coastal habitats of Africa. More specifically, the species' distribution seems to be directly related to the permanent upwelling systems of Mauritania and Angola throughout the whole year and, in a seasonal mesoscale process, in areas such as Guinea, the Angola Dome and the coastal upwelling off the coast of Ghana. The distribution model provides valuable information and reasonable prediction capabilities of the distribution of M. mobular, even when the presence sample size was limited. This new information is important to predict effects of overlap between the main fisheries and habitat important for pelagic species in this region, to help provide informed advice for effective fisheries management and contribute to the protection of this vulnerable species.