Habitat suitability of the Atlantic bluefin tuna by size class: An ecological niche approach

An ecological niche modelling (ENM) approach was used to predict the potential feeding and spawning habitats of small (5–25 kg, only feeding) and large (>25 kg) Atlantic blueﬁn tuna (ABFT), Thunnus thynnus , in the Mediterranean Sea, the North Atlantic and the Gulf of Mexico. The ENM was built bridging knowl- edge on ecological traits of ABFT (e.g. temperature tolerance, mobility, feeding and spawning strategy) with patterns of selected environmental variables (chlorophyll-a fronts and concentration, sea surface current and temperature, sea surface height anomaly) that were identiﬁed using an extensive set of pre-cisely geo-located presence data. The results highlight a wider temperature tolerance for larger ﬁsh allowing them to feed in the northern – high chlorophyll levels – latitudes up to the Norwegian Sea in the eastern Atlantic and to the Gulf of Saint Lawrence in the western basin. Permanent suitable feeding habitat for small ABFT was predicted to be mostly located in temperate latitudes in the North Atlantic and in the Mediterranean Sea, as well as in


a b s t r a c t
An ecological niche modelling (ENM) approach was used to predict the potential feeding and spawning habitats of small (5-25 kg, only feeding) and large (>25 kg) Atlantic bluefin tuna (ABFT), Thunnus thynnus, in the Mediterranean Sea, the North Atlantic and the Gulf of Mexico. The ENM was built bridging knowledge on ecological traits of ABFT (e.g. temperature tolerance, mobility, feeding and spawning strategy) with patterns of selected environmental variables (chlorophyll-a fronts and concentration, sea surface current and temperature, sea surface height anomaly) that were identified using an extensive set of precisely geo-located presence data. The results highlight a wider temperature tolerance for larger fish allowing them to feed in the northern -high chlorophyll levels -latitudes up to the Norwegian Sea in the eastern Atlantic and to the Gulf of Saint Lawrence in the western basin. Permanent suitable feeding habitat for small ABFT was predicted to be mostly located in temperate latitudes in the North Atlantic and in the Mediterranean Sea, as well as in subtropical waters off north-west Africa, while summer potential habitat in the Gulf of Mexico was found to be unsuitable for both small and large ABFTs. Potential spawning grounds were found to occur in the Gulf of Mexico from March-April in the south-east to April-May in the north, while favourable conditions evolve in the Mediterranean Sea from mid-May in the eastern to mid-July in the western basin. Other secondary potential spawning grounds not supported by observations were predicted in the Azores area and off Morocco to Senegal during July and August when extrapolating the model settings from the Gulf of Mexico into the North Atlantic. The presence of large ABFT off Florida and the Bahamas in spring was not explained by the model as is, however the environmental variables other than the sea surface height anomaly appeared to be favourable for spawning in part of this area. Defining key spatial and temporal habitats should further help in building spatially-explicit stock assessment models, thus improving the spatial management of bluefin tuna fisheries.

Introduction
Atlantic bluefin tuna, Thunnus thynnus (Linnaeus, 1785), is a highly migratory species able to tolerate wide ranges of environmental conditions (Arrizabalaga et al., 2015) in tropical and temperate waters of the Atlantic Ocean and Mediterranean Sea. The earliest scientific reference on that species comes from Aristotle in his treatise History of Animals, written in 350 B.C., describing the migratory and reproductive habits of tuna in the Aegean and Black Sea (D'Arcy Wentworth Tompson, 1947). The two main ABFT populations are known to spawn in the Gulf of Mexico (western stock) and the Mediterranean Sea (eastern stock) respectively while the productive waters of the North Atlantic are the main feeding grounds. Adolescent fish are known to feed on the shelf of the North West Atlantic, in the North East Atlantic and in the Mediterranean Sea (e.g. Fromentin and Powers, 2005;Katavić et al., 2013;Rooker et al., 2008b). The habitat at the scale of the entire distribution of the species remains relatively unknown despite being heavily exploited recently. Total ABFT catches increased dramatically after the mid 80s, from less than 20,000 t up to more than 53,000 t in the mid 90s. Introduction of annual quota regulations by the International Commission of the Conservation of Atlantic Tunas (ICCAT) gradually reduced these numbers below 15,000 t in the mid 2000s (ICCAT, 2014). The recent overexploitation patterns of fish stocks led managers to more frequently require information on the distribution of marine resources, enabling them to identify the areas of most suitable habitat (European Commission, 2008;Rubec et al., 1999). Understanding the dynamics and spatial distribution of species is crucial for management, as spatial variability governs the definition of management units, stocks and boundaries (Fromentin and Powers, 2005). Herein we provide an approach for accomplishing such a task by delivering an indirect identification of the Atlantic bluefin tuna habitats based on the association between environmental traits and presence data. The relatively recent introduction of remote sensing and access to relevant data from the scientific community allowed for incorporating environmental data into distribution and abundance analyses. This has confirmed that bluefin tuna distribution is significantly affected by spatial and temporal variations of environmental conditions (see e.g. . During their seasonal migrations, ABFTs seem to track changes in water temperature and currents, while they appear to preferably feed along frontal features (Druon et al., 2011;Royer et al., 2004;Schick et al., 2004). An attempt to associate presence data with several plausible factors affecting bluefin tuna distribution and abundance was recently undertaken by Druon et al. (2011) providing potential feeding and spawning habitats in the Mediterranean Sea.
In this paper, we link the ecological traits of small and large ABFTs to environmental variables (Ecological Niche Model approach, hereafter ENM) and investigate the respective feeding and spawning requirements. We used a large dataset of presence data and the literature for deriving the appropriate environmental envelope at the scale of the species distribution range (Gulf of Mexico, Mediterranean Sea and North Atlantic [defined as from 8.5°N to 74.0°N and from 82°W to 20°E outside of the other two areas]) and by size class (5-25 kg and above 25 kg, hereafter separating the small and large fish). The overlay of electronic tagging experiments and the potential habitats helped validate the model results and provide insights on the migration patterns important for understanding stocks dynamics. The seasonal and decadal habitat variability and spatial extent were discussed with respect to their potential impact on east and west ABFT stocks dynamics, as well as the utility on assessment and management.

Description of the ecological niche modelling
The methodological approach used in our ENM is essentially composed of four main steps (Fig. 1), namely: (1) identify the main behaviours and ecological traits of ABFT based on literature; (2) collect and process the ABFT presence data and environmental covariates by geographical area; (3) derive a cluster analysis to identify a suite of relevant thresholds of environmental variables related to the ABFT ecology that describe the feeding and spawning habitat characteristics and finally (4) develop a habitat model to classify on a daily basis the degree to which each portion of the study area (model grid cell) is either suitable or unsuitable for each habitat (environmental envelope). All variables were projected on the finest horizontal grid of the satellite ocean colour data which was used (NASA MODIS-Aqua sensor), i.e. at the resolution of 1/24°.

2.2.
Step 1 -specifying the ABFT habitats This first step of ENM consists of identifying the relevant ecological traits of ABFT that link behaviours to its environment. We conducted a literature review and assembled an important and widely distributed dataset of ABFT presence data across the geographical areas (see Section 2.3). Royer et al. (2004) and Druon et al. (2011) already found that feeding ABFT is preferably located in the vicinity of chlorophyll-a frontal features so that the horizontal gradient of chlorophyll-a (hereafter gradCHL) is used as a proxy for food availability. A specific range of chlorophyll-a concentration (CHL) is also associated with that proxy. ABFT is one of the rare fish species to be thermo-regulated at about 20°C (e.g. Carey et al., 1971;Carey and Lawson, 1973) and is able to occasionally dive to depths of >1000 m (Block et al., 2001) but spends most of its time in surface waters (79 ± 8% in the first 50 m from tagging studies, Walli et al., 2009). Large ABFTs have thus a rather large tolerance for temperature although it appears to be an important constraint for juvenile fish (Galuardi and Lutcavage, 2012). Therefore a specific range of sea surface temperatures (SST) was introduced to account for the thermal tolerance of both size classes while they feed. Sea surface height anomaly (hereafter SSHa) was tested as a variable potentially impacting the distribution of feeding habitat of both size classes. SSHa is mainly influenced by seasonal changes in temperature and currents that create eddies and gyres, i.e. divergent and convergent areas, potentially responsible for enhanced primary productivity and tuna prey aggregation (see examples in Hobday and Hartog (2014)). Arrizabalaga et al. (2015) and Teo and Block (2010) notably found that the temperate ABFT grows in colder and more productive environments with negative SSHa compared to tropical tuna species (with near null or positive SSHa).
ABFT spawning is known to occur in rather warm and mostly oligotrophic surface waters during spring. Parts of the Gulf of Mexico and Mediterranean Sea present such characteristics (Fromentin and Powers, 2005;Rooker et al., 2007). A specific CHL and SST range, as well as the monthly increase of SST (DSST 30 ), which simulates the spring stratification build-up, were selected to represent likely suitable spawning grounds (Druon et al., 2011). Additionally, a preferred range of SSHa was introduced for its potential role in food and larvae retention (see e.g. Bakun, 2013) as well as in transport from oligotrophic to mesotrophic areas. Intermediate levels of Eddy Kinetic Energy (EKE), which is derived from SSHa, were shown to strongly favour pelagic fish spawning (ABFT - Bakun, 2013;Teo et al., 2007;Teo and Block, 2010;Tuna species -Reglero et al., 2014; Small pelagic species -Asch and Checkley, 2013).
In order to reflect feeding opportunities within the mesotrophic (e.g. Mediterranean Sea) to eutrophic (e.g. northeast Atlantic) environments in which ABFT feed, we added an intermediate level of productivity, so that we considered three levels of feeding habitat: highly, moderately and poorly productive with daily values of 1, 0.3 and 0 respectively. The definition and parameterization of highly and moderately productive habitats and the model equations for the feeding and spawning habitats are detailed in the Supplementary Information (Section 8.4).
The boundary limits for suitable feeding habitats were defined for each ABFT size class by taking the widest favourable conditions observed spatially (typically the percentile range of 1.5-5th and 95-98.5th) in cases where the presence data was believed to appropriately cover the maximum range of an environmental variable (e.g. CHL and gradCHL). When extrema of a given variable as derived from the presence data were absent, expert knowledge was used, albeit a broader plausible range was specified than when using the literature (for e.g. the SST minimum). Spawning habitats were distinctively defined for the Gulf of Mexico and Mediterranean Sea (see Section 4). We used a cluster analysis on the presence data to objectively segregate the feeding from the spawning behaviours and to enhance habitats that are potentially underrepresented in the data. Then we set the boundary values of the spawning habitat using low and high values of the environmental variables (see details below in Step 3). In all cases, the habitat thresholds were chosen to be consistent with the suitable range across areas, with the range obtained by the cluster analysis considering the presence data distribution and with literature. For instance, the 97th percentile value for the maximum SST of the feeding habitat of small ABFT was preferred over the 85th percentile value (lower SST value) of the warmest cluster in the Mediterranean Sea because the warmer eastern basin was not well represented in the ABFT presence data. In a few cases, extreme value gathered in the literature has driven the parameter selection due to insufficient (seasonal or geographical) coverage of presence data.

Step 2 -data
The second step of our framework focuses on the collection and suitable preparation of input data for the model.

ABFT presence data
A total of 31,362 non-redundant presence data of ABFT with precise geolocation were collected between 1997 and 2014 in the studied areas; 23,644 in the North Atlantic, 5324 in the Mediterranean Sea and 2394 in the Gulf of Mexico, respectively. The presence data originates mostly from commercial fisheries (drifting surface longliners, purse seiners, rod and reel, bait boats and aerial surveys) although a few areas were observed by scientific surveys (e.g. Gulf of Lions with aerial surveys and tagging, Bay of Biscay with acoustic surveys and tagging). The conventional tagging data and the start and end locations of electronic tagging data collected by ICCAT were also included in the analysis. The location for presence data corresponds to exact GPS coordinates except for surface longliners (13% of data, two thirds of which was in the Gulf of Mexico and one third in the eastern Mediterranean) and baitboats (6% of data in the Bay of Biscay) whose position precision was estimated to be below 10-15 NM. Redundancy filtering ensured that observations on the same day were separated by more than 2.3 km (about half of the model cell). We used ABFT presence data to define the environmental envelope of each size class. We extracted the environmental information in different years depending on ocean model and satellite data availability (see next section), i.e. from 1997 to 2012 for the physical variables, from mid-2002 to 2014 for chlorophyll-a products (see Table 1 for details). Weight or length information was available for 64% of the data with 17,455 classified as large ABFT (hereafter above 25 kg) and 2589 as small ABFT (hereafter between 5 and 25 kg). When only length was available, we applied the length/weight relationship used by the ICCAT working group to estimate the corresponding weight (see ICCAT, 2013). We selected 25 kg because it corresponds to the mean weight at first maturity in the Mediterranean Sea and East Atlantic. Below this weight ABFTs are expected to only exhibit feeding behaviour. The preferred habitat in each area corresponds though to the largest individuals present in the  database (giant bluefin, hereafter above 200 kg) since they have a larger environmental envelope than that of young adults (e.g. tolerance to low temperature). As for complementary information to Table 1, the percentage of presence data among which weight was available is 13% from 5 to 25 kg, 15% from 25 to 100 kg, 26% from 100 to 200 kg and 46% above 200 kg. The presence data were not uniformly distributed across space and seasons. Consequently, we performed the environmental analysis of large known feeding and spawning areas by employing a statistical method (cluster analysis) that allowed the enhancement of underrepresented habitats. Table 1 presents, by area, the amount of presence data used to extract environmental information. The amount of data in the cluster analysis was reduced since all environmental variables must be available, including the cloudfree CHL data. Among the data with available weight (Table 1), about half of the sightings in the Mediterranean Sea were small ABFTs while most of the data in the North Atlantic were large fish and 99% in the Gulf of Mexico were giant bluefin (5th percentile value is 116 kg for the Mexican data, n = 483). The mean fork length (FL) of ABFT caught by longliners in the northern Gulf of Mexico is 226 cm FL ± 41 (ca. 238 kg, ICCAT data for area BF60, n = 2354), and therefore assumed that when weight was missing it was likely above 25 kg.
Overall, the ABFT data is mainly concentrated relatively close to land while the central and extreme latitudes of the North Atlantic show low levels of sightings (Fig. 2, see the Supplementary Information in Section 8.1 for distribution by size class and month).

The electronic tagging data
Internal archival tags record water temperature, depth and light intensity many times daily. This data was used to estimate (Lam et al., 2008) and interpolate (if no geolocation was available, Akima, 1978) the daily location of the fish. Because they are surgically implanted in the peritoneal cavity, the fish has to be recaptured to recover the detailed information recorded each minute by the tag. Pop-up archival tags (PSAT) record the same information as the archival tags, except the internal temperature. They are fixed by a tether to the dorsal musculature near the second dorsal fin and detach after a pre-set time interval. Upon release, the PSAT transmits a summary of the recorded data to the closest Argos satellites.
Information on the electronic tags used in this paper was summarized in Table 2. Several electronic tags were overlaid as independent data with the potential habitat maps from the monthly to the seasonal climatology (2003-2012) time scales in order to enhance the potential use in model validation and behavioural analysis.

Chlorophyll
Surface chlorophyll-a concentrations and fronts were used at a daily time scale from the MODIS-Aqua ocean colour sensor (http://modis.gsfc.nasa.gov) for the period from mid-2002 to 2014. The MODIS spatial resolution of 1/24°was used to identify meso-scale CHL fronts and it defines the reference grid of the habitat model. It is worth noting that any other CHL data on a daily time scale with a similar resolution may be included to increase the coverage and the length of the time series (e.g. SeaWiFS from September 1997 to December 2010), provided a specific calibration is made. We chose MODIS-Aqua as it provided the most recent time series of CHL available. Daily CHL data were pre-processed using iterations of a median filter in order to recover missing data on the edge of valid data. The median filter and Gaussian smoothing procedure (see Druon et al. (2012) for details) allowed for the recovery of ca. 8% of the CHL data. The relative gain in coverage was however much higher after the gradient calculation of CHL (+38%) and the habitat computation (+57%). The front enhancement of daily CHL data was calculated with an edge-detection algorithm which was shown to perform better than the histogram methods in detecting horizontal gradients given clear viewing conditions (Ullman and Cornillon, 2000). Note that the daily time scale was required here to allow the identification of CHL fronts which would be blurred or would disappear if using time-integrated data. We therefore only used the daily data for computing CHL fronts and the habitats. We extracted a 3-day mean CHL value however for the statistical analysis in order to substantially increase the number of presence data and analysis robustness. The 24 h variability of CHL level was thus assumed to be low.

Surface height anomaly, current velocity and temperature
SSHa, SSCV and SST were extracted from ocean model of the MyOcean Consortium (http://www.myocean.eu), a core marine service within the European Global Monitoring for Environment and Security (GMES) Program whose objective is to develop an integrated capacity for ocean monitoring and forecasting. Monthly mean data of sea surface height anomaly, temperature and current velocities for the period 1997-2012 were extracted from a Mediterranean Sea hydrodynamic model (NEMO-OPA v3.2) at the horizontal resolution of 1/16°and 72 unevenly spaced vertical levels. For the North Atlantic and the Gulf of Mexico areas, the same variables were extracted from a global model (Glorys2V3) at 1/4°horizontal resolution and 75 unevenly spaced vertical levels. Both models include a variational data assimilation scheme for temperature and salinity vertical profiles and satellite sea level anomaly (Oddo et al., 2009). Original data were interpolated on the MODIS-Aqua grid. The monthly data were linearly interpolated to daily values. Such monthly to daily interpolation is believed to produce suitable estimates of the seasonal changes that define ABFT habitats. SST was taken from the upper model layer (ca. 3 m) while SSCV was taken as the mean of the upper layers of the MyOcean models (ca. 13.5 m) in order to capture the transport of the mixed layer. The current intensity was included in the habitat model as a directionless quantity. The sea surface height anomaly from the Mediterranean model was corrected by the mean basin-wide difference found with the SSHa of the global model (À0.2842 m) to ensure a consistent calibration of the habitat model. The distribution of SSHa (all-year-round in relation to feeding) and SST and monthly difference of SST (prior and during the spawning season) are shown in Supplementary Information (Fig. SI-2).

Step 3 -ABFT environmental analysis
The third step of our ENM involved (a) identifying a set of environmental variables linked with ABFT behaviours and (b) exploring their variability to identify relevant environmental threshold values that separate favourable from unfavourable habitat. This analysis was made from 1997 to 2012 for the physical variables of MyOcean and from 2003 to 2014 for the CHL data of the MODIS-Aqua sensor.
The link of each selected environmental variable with ABFT presence was analysed with a cluster analysis following the procedures reported in (Berthold et al., 2010) and (Hartigan, 1975). The cluster analysis is particularly well suited for separating behaviours that occur in distinct environments (presently feeding and spawning) and for highlighting habitats that are marginally represented in the data and that may be interpreted as outliers with other statistical methods. The analysis was derived by area (Gulf of Mexico, North Atlantic and Mediterranean Sea) and by size class (small and large fish) except for small ABFT in the Gulf of Mexico for which there was no presence data (see Supplementary Information for more details (Section 8.3)). The selected variables were the 3-day mean CHL (log transformed), the 3-day mean horizontal gradient of CHL (gradCHL, log transformed), SST, the monthly increase of SST (DSST 30days ), SSCV, SSHa, month and weight. We tested from two to five clusters and retained the number that allowed the clearer and simpler interpretation of the ABFT seasonal behaviours in each area and size class. The selection of relevant thresholds for the habitat model was either driven by the cluster analysis (mostly spawning) or by the overall environmental envelope in situations where the ABFT presence data was suitably distributed over the seasons (partly feeding). The parameterization of the feeding habitat was completed using the literature and expert knowledge (e.g. minimum temperature). The characteristics defining the feeding habitat were consistent across the studied areas so that productive frontal features with the same CHL and gradCHL values would correspond to the same feeding capacity from one area to another. Conversely, area specific characteristics were defined for the spawning habitat because important hydrological differences exist between the Gulf of Mexico and the Mediterranean Sea, especially with respect to temperature and currents. In situations where the cluster analysis was used to define the relevant thresholds, the 15th and 85th percentile values were used as they represent relatively extreme environmental boundaries while rejecting the potentially misclassified distribution tails. Because the relative size of distribution tails is larger for the cluster analysis (lower number of presence data) than for the global dataset analysis, we chose a narrower range of percentile values for the cluster-derived thresholds to best define the core habitat (e.g. 85th compared to 97th). Care was taken to respect the consistency of selected thresholds values between the regional analysis and the literature, recognizing that our presence data was relatively unevenly distributed across seasons and geographical areas.

Step 4 -formulation of the ecological niche model
Once the environmental variables were selected and the threshold values set, the last step consisted of defining the specific ecological niche of the small and large ABFT, using the areas of favourable biotic conditions (represented by CHL and gradCHL) and abiotic preferences (SST, SSHa and SSCV). The favourable environmental envelopes predicted the daily suitability of cells within the habitat for ABFT feeding (small and large fish) and spawning (large fish) on a scale of 0-1. In practice, the favourable habitat for each behaviour are cells that meet all the suitable ranges of selected variables as defined by the parameterization (see below) and equations (see Section 8.4). The areas meeting the daily biotic and abiotic requirements of the habitat model were then integrated over time to yield seasonal suitability maps of the relative frequency of occurrence (i.e. the sum of daily habitat values [from 0 to 1] over the number of days for which the habitat was effectively estimated). The model performance was estimated by computing the distance between the presence data and the closest favourable habitat (3-day composite) for the period from mid-2002 to 2012, first by size class (6852 presence data for the large fish, 952 for the small fish) and then combined (11,312), including ABFT data for which no weight was available (+45%). For the latter exercise, the shortest distance to either small or large fish habitat was selected.

Distribution of presence data versus habitats
The presence data mainly covered the period from March to October (90%) with a peak from July to October (70%). The period from November to February was poorly represented (10%). The observations were relatively well distributed by area and season allowing the modelling of both the feeding and spawning habitats with more than 90% of data related to large ABFT from January to June in the Gulf of Mexico and from May to September in the Mediterranean Sea. Indeed, the temporal coverage included both spawning periods (from end of March to May in the Gulf of Mexico, (Rooker et al., 2007), and from mid-May to mid-July in the Mediterranean Sea) and a relatively large portion of the seasonal variability in the foraging habitat. The presence data in the North Atlantic was mostly (ca. 90%) from August to October for large fish and from July to October for small fish. The main data gaps were the winter months (all areas) and extreme latitudes in the North Atlantic so that we also used the literature to define these extreme environmental boundaries (see Section 2).

Environmental conditions of ABFT habitats by area and size class
The cluster analysis resulted in a descriptive classification of small and large ABFT habitats at the scale of each studied area. The analysis was performed by area and size class in order to group observations by main seasonal behaviours and to identify regional peculiarities. Three clusters proved to be the minimum number favouring interpretation for all considered areas and size classes. Clusters were ranked by decreasing group size so that cluster 1 was always the largest.
Small ABFTs were mainly in areas of productive fronts (high gradCHL) and intermediate CHL levels with substantially higher values in the North Atlantic compared to the Mediterranean Sea (see Supplementary Information, Fig. SI-3). The SST range was generally from 15 to 25°C depending on the season while SSHa values were always negative (mostly below À0.3 m). The Mediterranean Sea clusters represent the seasonal feeding conditions in spring and summer/autumn of two size classes of small ABFT, with median weights of 15 and 8 kg respectively. ABFT with median weights of 15 kg tend to show a wider geographical distribution in summer/autumn with lower CHL levels than smaller fish. The clusters in the Atlantic only characterize summer and autumn feeding conditions of small and intermediate juveniles (median weights of 6 and 12 kg) in the north-east basin and of intermediate size (median weight of 12 kg) in the Bay of Biscay. Note that the highest CHL and lowest temperature levels correspond to the largest juveniles in both areas.
The clusters of large ABFTs show a broader range of environmental conditions than clusters of small fish (see Supplementary  Information, Fig. SI-4). CHL and gradCHL levels were similar or substantially higher for large fish except for cluster 1 in the Mediterranean Sea and cluster 2 and 3 in the Gulf of Mexico which corresponds to particularly low levels in spring. These clusters were also distinct from the others by at least one of the following: high SST levels, a high monthly increase of SST and a narrow range of SSHa. The smaller cluster 3 (n = 82) in the Gulf of Mexico, which spans from May to November, is heterogeneous for most variables.
It notably has the highest SSCV levels, thus presumably in the Gulf Stream, and mostly occurs during the seasonal migration. Note that, similar to the small fish, the highest absolute levels of CHL and gradCHL characterize the cluster of the largest ABFTs in the Atlantic (cluster 3, median weight of 310 kg).
Taking into account the elements of ABFT ecological traits described in the methods, the clusters characterized by relatively high gradCHL levels and CHL above 0.1 mg m À3 were associated with feeding behaviour, while those with low gradCHL levels, CHL below 0.15 mg m À3 , high levels of SST and DSST 30days and a narrow range of SSHa were associated, for the large ABFTs, with spawning (cluster 1 in the Mediterranean and cluster 1 and about half of cluster 3 in the Gulf of Mexico for the large ABFTs, Fig. SI-4a and c).

Habitat modelling and parameterization
The cluster analysis described a wide range of trophic conditions in which small and large ABFTs feed, from oligotrophic in the Gulf of Mexico and the Mediterranean Sea to eutrophic in the Atlantic. Table 3 presents the habitat parameterization by size class and behaviour resulting from the cluster analysis, literature and expert knowledge (see Section 2).
We used the 15th and 85th percentiles of clusters related to spawning as boundary values for this habitat except for the maximum SST and SSHa thresholds in the Gulf of Mexico since cluster 3 in Fig. SI-4c is believed to characterize first the spawning in May and second the post-spawning migration. We chose the median values for these two cases taking into account the seasonal increase of SST and SSHa from spring to summer (Fig. SI-4c). The intermediate thresholds for CHL and gradCHL defining the levels of productive habitat were chosen using the cluster analysis and the differences between the oligo-and eutrophic environments (Gulf of Mexico and Mediterranean Sea vs North Atlantic, see Fig. SI-4). The minimum SST for the feeding habitat of large ABFTs was defined by literature (7.5°C in relation with SST fronts and fish dives, see Table 3) and by regional expertise for the small ABFT (13.0°C, which was close to the winter minimum in the western Mediterranean Sea). The maximum SST threshold for large ABFT feeding of 24.0°C was set as an intermediate level between the 85th percentile value of the winter habitat in the Gulf of Mexico (23.83°C, n = 141, red cluster) and the same value for the summer-autumn habitat in the Mediterranean Sea (24.37°C, n = 149, red cluster). The maximum SST threshold for small ABFT feeding of 26.1°C correspond to the 97th percentile value (n = 2105) while it was close to the 85th percentile value of the summer cluster in the Mediterranean Sea (25.72°C). The minimum SSHa level for the potential feeding habitat of both size classes of À0.10 m correspond to the 98.75th, 99.5th and 97th percentile values of the small, large and undifferentiated-weight fish respectively.

Outputs of the habitat model
We present in this section the spatio-temporal distribution of model results and evaluation. The 75th percentile distance from presence data to the closest favourable habitat was generally below 6 km (see Supplementary Information Section 8.5 and Table SI-1) except when large ABFTs were potentially migrating to and from the spawning grounds (Mediterranean Sea and Gulf of Mexico). In this case, the 75th percentile distances were between 16 and 133 km and the median distances were all zero, i.e. at least half of presence data were in the preferred habitat. The seasonal distribution of observations and the corresponding distances to habitat (see Supplementary Information Fig. SI-6) show that most of the large distances occurred when the amount of presence data was low, i.e. presence was possibly occasional.
When the presence data was abundant, the monthly mean distance to the closest habitat was generally below 15 km, except in July in the Mediterranean Sea (35 km) and in May-June (55-190 km) in the Gulf of Mexico when a fraction of fish were likely migrating. The model fit was best for May to November in the Mediterranean Sea and North Atlantic and for January to May in the Gulf of Mexico (Fig. SI-6). These periods indeed corresponded to relatively low monthly mean size of potential habitat with values from 1% to 12% (5.1 ± 3.5%) of the considered surface area for potential feeding and below 7% (2.4 ± 2.4%) for potential spawning (see Section 8.6).
The small ABFT potential habitat ( Fig. 3a and b) was mostly located in temperate waters and off northwest Africa in the upwelling area (down to 13°N) with a seasonal latitudinal oscillation up to the North Sea on the east Atlantic and to the Gulf of St. Lawrence in the western basin. Small ABFT feeding was predicted to be more important in summer than winter in the North Atlantic while in the Mediterranean Sea (western Alboran Sea, Gulf of Lions, western Adriatic Sea and northern Aegean Sea) potential feeding grounds were restricted to hot spots during the summer and were absent in the Gulf of Mexico. Their potential feeding habitat was widespread during winter in the western and northern Mediterranean Sea and near the shelf of the Gulf of Mexico. A notable spatial discontinuity between the winter potential habitat for small fish in the Gulf of Mexico and the closest summer feeding grounds in western Atlantic was observed. Spring and autumn habitats (not shown) had intermediate distributions between winter and summer.
The potential feeding grounds of large fish ( Fig. 3c and d) were more extended towards the north in the Atlantic compared to small ABFT, especially during summer with an extension up to the Labrador Sea in the western basin and to the Norwegian Sea in the eastern basin (up to 73°N). The seasonal distribution of feeding habitats was similar for both size classes in the Gulf of Mexico, the Mediterranean Sea and off north-west Africa.
The potential spawning habitat mainly occurred from mid-May to mid-July with a peak in June in the southern and western Mediterranean Sea and from April to May in the northern and southwest Gulf of Mexico (Fig. 4). The potential spawning grounds in the North Atlantic (off north-west Africa and near the Azores Islands) in July and August were the product of the model calibrated with the presence data of the Gulf of Mexico (no substantial habitat was obtained using the parameterization of the Mediterranean Sea). The existence of these potential spawning grounds was not supported by presence data. The monthly size of potential habitats for large ABFTs in each of the three main study areas (Gulf of Mexico, North Atlantic and Mediterranean Sea) showed high levels of variability at the seasonal and inter-annual scales (see details in Fig. SI-7). The maximum extent of potential habitat occurred in summer in the North Atlantic while the maximum size of potential feeding and spawning habitats occurred in winter and spring respectively in the Gulf of Mexico and Mediterranean Sea.
Habitat maps in Figs. 5 and 6 are similar to Figs. 3 and 4 but focus on areas where presence data and electronic tags were available. Most of the available presence data for small ABFT were available in summer (Fig. 5a, squares in magenta colour) and match the predicted habitat well, mostly in the Mid-Atlantic Bight and south of Gulf of Maine, in the Bay of Biscay and the Gulf of Lions. Panels a and b of Fig. 5 also present the overlay of archival tag most probable tracks (interpolated position for missing geolocation days) of a 5-26 kg fish initially tagged in the Bay of Biscay with the corresponding mean summer and winter potential habitat. This young fish showed a strong fidelity to the Bay of Biscay from June to October despite a priori favourable grounds further north, while the winter feeding grounds were located further south (36-43°N) and west (51-9°W) off the Azores and Portugal at the eastern area influenced by the Gulf Stream.
The overlay of presence data of large ABFT with their potential feeding (Fig. 5 panels c -summer and d -winter) and spawning habitat (mostly spring, Fig. 6) emphasizes the important north (summer) -south (winter) oscillation of populations. Note that observations were not a priori linked to a specific behaviour so that each position of large ABFT may correspond with feeding, spawning or migration. The presence data of large fish that agree with the predicted habitat are north of the Gulf Stream for latitudes above 30°N, in the Bay of Biscay and the Gulf of Lions during summer while they are located south of 36°N near Cape Hatteras and in the Gulf of Mexico during winter for feeding and spring for spawning. The overlay of two large tagged fish matches this oscillation from the Gulf of Maine area to Cape Hatteras and from the Gulf of Lions to the Balearic Islands area. The presence data that does Table 3 Model parameters defining the small and large ABFT habitats in the Mediterranean Sea (Med), North Atlantic (nAtl) and Gulf of Mexico (GoMex). The feeding habitat was identically parameterized in the three areas using the overall environmental preference while the spawning habitats in the Gulf of Mexico and the Mediterranean Sea have distinct settings (mostly 15th and 85th percentile values of the cluster analysis, see text for details). Model parameters are overlaid with the variable distributions in the Supplementary Information (Fig. SI-2 not correspond to any prediction of habitat are located in oligotrophic environments off the Bahamas and Florida in winter (and spring for which feeding habitat is not shown) and in the eastern Mediterranean Sea off Greece and Turkey in summer. In July 2011, a tagged fish in the North East Atlantic spent about two weeks in a potential spawning area identified by the model to be south of the Canary Islands (Fig. 6c) a priori unfavourable to feeding (Fig. 3d). The inter-annual variability of habitat size for both feeding and spawning is contrasted in Fig. 7. A substantial decrease of potential feeding habitat occurred mainly from 2010 to 2012 compared to previous years in the Mediterranean Sea (30-55%) and to a lesser degree in the North Atlantic (

Modelling methods
We applied important methodological changes compared to the ABFT potential habitat model of Druon et al. (2011). Notably, we substituted the satellite-derived sea surface temperature (SST) with SST from ocean circulation models because satellite SST was found to be highly variable from day to day and may not represent the mixed layer temperature well. We thus preferred the more stable ocean model SST for representing seasonal changes. For the same stability reasons, we only used the horizontal gradient of surface chlorophyll-a content (CHL) to identify productive frontal features instead of both gradients of SST and CHL simultane-ously. The statistical analysis used to identify the relevant environmental thresholds was also modified. Here we preferred a cluster analysis over an iterative optimization process (minimization of distances between presence data and habitat) as the method based on clusters was less dependent of the presence data distribution and it emphasized underrepresented groups (such as behaviours or seasons). Finally the model was considerably extended to the North Atlantic and the Gulf of Mexico using a much larger presence dataset globally and regionally (four times more in the Mediterranean Sea), including two ABFT size-classes instead of one and using more environmental variables (namely sea surface currents and sea surface height anomaly).
The parameterization of the potential feeding habitat is identical between areas, i.e. using the high and moderate productive habitats, making the reasonable hypothesis that chlorophyll-a fronts of same size (and thus of comparable duration) are similarly productive for the food chain independently of their location. In contrast, the potential spawning habitat was specifically parameterized by area since seasonal solar heating is dominant in the (a) Small ABFT potential habitat for feeding in SUMMER (July-September) (b) Small ABFT potenƟal habitat for feeding in WINTER (January-March)  Table 2 for e-tag details. Note that presence data and e-tags are not a priori associated to a specific behaviour and that presence data in the Gulf of Mexico and off Florida were associated for these maps to large fish (see text for details). The 200 m depth contour is shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Mediterranean Sea (enclosed basin) while the horizontal transport of heat by the Gulf Stream also occurs in the Gulf of Mexico (see Fig. SI-2).
The limitations of the approach mostly reside in cloud cover and distribution of presence data. Cloud coverage, that hampers the detection of CHL, is variable by region (e.g. Atlantic vs Mediterranean) and seasons (winter vs summer). Consequently, for time composite products, we used the frequency of habitat occurrence corrected by the number of possible identifications. The heterogeneous distribution of presence data by region and season may, if not properly taken into account, lead to a biased parameterization. To circumvent that issue, we first used a cluster analysis to highlight underrepresented habitats (mostly for the spawning habitat) and, in cases where presence data was lacking in marginal feeding environments (e.g. low SST), we sought ABFT observations in the literature to ensure the environmental envelope was representative of the core habitat at the population level. The addition of independent presence data (i.e. without information on fish weight) did not increase the 50th, 75th and 90th distance to the closest habitat (see Table SI-1). Finally, electronic tagging data were overlaid on the habitat maps to interpret ABFT behaviour and to show potential use in model validation, albeit with consideration given to location estimation uncertainty. The robustness of this ENM approach relies first on the primary hypothesis that feeding predators are preferably located in the vicinity of CHL fronts (see same approach used for tropical tuna (c) Large ABFT potential habitatfor feedingin SUMMER (July-September) (d) Large ABFT potential habitat for feeding in WINTER (January-March)  Druon et al., 2012, andhake recruits in Druon et al., 2014) and, second, on the cluster analysis that allowed objective segregation of presence data by habitat type as well as subsequent thresholding.
Only a few similar studies on temperate tuna habitat that used fronts are reported in the literature. Nieblas et al. (2014) used SST and CHL range and fronts and EKE to highlight potential spawning grounds of the Southern bluefin tuna in the southeast Indian Ocean. Their results are consistent with the present study notably regarding favourable conditions for SST and CHL (high and low levels respectively). Their findings on low EKE levels and low front intensity as suitable for reproduction are as well indirectly consistent with our results that enhance intermediate levels of current velocity and SSHa. The simple temperature-driven model of Hartog et al. (2011) used to model the habitats of the Southern bluefin and yellowfin tunas in the eastern Australia does not describe the core habitat of these two species but is efficient to identify the habitat overlap.

ABFT habitats and behaviour
In comparison with the previous modelling work in the Mediterranean Sea (Druon et al., 2011), the results of the present study show an overall spatio-temporal consistency, but with a substantially higher continuity and stability of the potential spawning habitat that can be linked to the substitution of satellite SST by oceanic modelled SST. The potential feeding habitat is also globally consistent although wider in winter due to a higher value for the CHL maximum which was previously set using a limited number of presence data.
The environmental envelope of the small ABFT habitat agrees with the estimated positions of the tagged fish in the eastern North Atlantic (Fig. 5), especially with respect to the minimum temperature limitation that appears to drive the winter migration off Portugal and the Azores. Note however that the entire potential habitat does not seem to be exploited by small ABFTs in summer. We suggest that presence of preferred preys (e.g. anchovy and sardine) and limited energetic cost prior to the southward/westward winter migration may favour the fidelity of small ABFT to the Bay of Biscay. In addition to that area, a large presence of small ABFT during summer in the northern area of the western Mediterranean Sea (Liguro-Provençal basin) and of the temperate U.S. shelf and Gulf Stream margin (37-43°N, Mid-Atlantic Bight and southern Gulf of Maine) highlighted the importance of these grounds for feeding in our study and is consistent with Galuardi and Lutcavage (2012). The main winter grounds for potential feeding spread further south, e.g. off Florida as attested by a recreational fishery and in the whole western Mediterranean Sea, in the north of the eastern Mediterranean (Levantine Sea), off Portugal and the Azores and off Morocco down to Senegal in agreement with reported observations (Cermeño et al., 2015;Fromentin and Powers, 2005;Galuardi and Lutcavage, 2012;Mather et al., 1995). The winter potential habitat of small ABFT in the Gulf of Mexico is yet unlikely to correspond with an effective habitat since the unsuitable summer conditions (SST levels [not shown] substantially above the model threshold of 26.1°C, in agreement with electronic tagging observations of Galuardi and Lutcavage (2012), and positive SSHa levels) would constrain them to seasonal migration into the Atlantic. We suggest that the large distance from summer grounds in the mid-Atlantic Bight, the high current velocities of the Gulf Stream and resultant energetic cost may limit the migration of small ABFT to the Gulf of Mexico in winter in preference of the suitable habitat of the more easily accessed shelf of the South Atlantic Bight. The lack of targeted fishing effort on small ABFT in the Gulf of Mexico however prevents us from asserting that this area is not an effective habitat in winter.
The wide habitat preference of large bluefin tunas, their complex behaviour and capacity to travel further than small fish was well predicted by the model although some ABFT presence in oligotrophic environments remained unexplained. In particular, known feeding patterns of large and giant ABFTs were well represented by the habitat model, such as the summer feeding grounds in the central North Atlantic, the northeast Atlantic continental shelf, the Gulf of Maine or the Gulf of St Lawrence. An exception is the central and eastern part of the Mediterranean Sea in summer where ABFT bycatches have been reported in various large pelagic fisheries operating in the area (Damalas and Megalofonou, 2012;Di Natale and Mangano, 2008;Fenech-Farrugia et al., 2004;Peristeraki et al., 2008) despite the apparent absence of favourable feeding habitat. We suggest that this false negative of the model for feeding in surface waters corresponds to a population that, as an alternative to the common epipelagic presence and foraging, is notably feeding in mesopelagic waters as attested by the stomach content study of Karakulak et al. (2009). The locations of the main bluefin tuna spawning areas in the Mediterranean Sea, the waters surrounding the Balearic Islands and Malta, the South Tyrrhenian Sea and in the northern Levantine Sea (Corriero et al., 2003;Garcia et al., 2003;Heinisch et al., 2008;Nishida et al., 1997;Oray and Karakulak, 2005;Piccinetti et al., 1997) and in the northern Gulf of Mexico (e.g. Muhling et al., 2011;Rooker et al., 2007) were correctly represented by the model both spatially and temporally (from mid-May to mid-July and from end of March to end of May respectively). In particular, the ABFT migration from the northwest Atlantic to the Gulf of Mexico and back as described by electronic tag experiments (e.g. Wilson et al., 2015) is fully consistent with the dynamics of potential feeding and spawning habitats (see Fig. SI-7 and related section).
The model also predicted uncharacteristic or secondary potential habitats that have been recently documented. These potential feeding areas were detected along the northern boundary in the Norwegian Sea in the 50s and the 60s (Fromentin, 2009), near Iceland in summer (MacKenzie et al., 2014) and during winter in the Azores area, off Morocco to Senegal (Fromentin and Powers, 2005;Suzuki and Kai, 2012). Conversely, the continental shelf of the northern Gulf of Mexico was shown by the model to be potentially favourable for feeding in winter while no such evidence was found in the literature. Feeding was however reported off the shelf in spring (Butler et al., 2015) and we suggest that if large ABFTs migrate to the Gulf of Mexico during winter for spawning in spring (Wilson et al., 2015), they do not typically feed along the shelf area where low SST levels (below 20°C) are likely to hamper the gonads' maturation (Medina et al., 2002;Schaefer, 2001). This probable false positive of the winter feeding habitat on the shelf of the northern Gulf of Mexico may thus result from a complex interaction with reproduction. The use of the Gulf of Mexico parameterization to predict potential spawning habitat in the Atlantic yielded potential spawning grounds in the Azores area, off Morocco and south of the Canary Islands which persisted throughout July and August (Figs. 4 and 6). These secondary potential spawning grounds are supported by limited historical reports (Di Natale et al., 2013) and electronic tag data from a large ABFT in July 2011 (Fig. 6). Similarly, the Central Ionian could be a potential spawning area and recent electronic tagging experiments (Cermeño et al., 2015) have indicated this possibility. More sampling in these areas is necessary to confirm these potential reproduction grounds.
The overlay of a few electronic tags notably showed the feasibility of validating the environmental envelopes estimated by the model. Additional overlay of tracks was beyond the scope of this study and may the subject of future analysis as consistent processing of geolocation and error ellipse is required. This exercise also yielded valuable information on ABFT behaviour. In particular, the small ABFT at liberty for two years in the north-east Atlantic (tag #390540) allowed one to visualize the progress in habitat extent for the same fish between its first winter migration in the Azores area ($25°W) and its successive migration further east ($50°W) using a route at the northern boundary of preferred habitat ($48.5°N). This overlay demonstrate the high capacity of small ABFTs to increase the exploration of their environment in winter while maintaining a strong fidelity to their summer feeding ground (presently the Bay of Biscay). The same type of fidelity to the summer ground in the Gulf of Lions was also shown by the tag of a young adult (80 kg) that spent most of the winter and early spring in the southern part of the western Mediterranean Sea (Fromentin and Lopuszanski, 2014).
The presence data off Florida and the Bahamas from January to June in Fig. 2 (n = 379) for which no suitable habitat was predicted and no individual weight was available corresponded with large ABFTs (213 cm FL ± 43, ca. 199 kg, ICCAT data for area BF61, n = 2309) mostly caught by longliners. Electronic tags in that area and period (result not shown) revealed that fish perform daily dives in the mesopelagic environment to about 500 m suggesting a deep feeding activity. The environmental conditions in April and May in the area off Florida and the Bahamas were favourable for spawning with respect to SST, DSST 30days and CHL, however SSCV was generally above the preferred range and, above all, SSHa levels were unquestionably higher (in the range from +0.1 to +0.4 m, see Fig. SI-2) than the preferred range (percentile values 3-97th are À0.58 and À0.10 m, n = 27,253). Positive and near null SSHa levels, which are mostly driven by high SSTs and linked with relatively low productive environments, effectively segregate worldwide the habitat of tropical from temperate tuna species (Arrizabalaga et al., 2015;Teo and Block, 2010). The area off the Bahamas and Florida was the only environment where large ABFTs were observed in positive SSHa values, and this, at times coincided with migration for spawning. The small number of larvae that was found north and east of the Bahamas in May 2013 (Lamkin et al., 2014) and the above environmental analysis suggest that it may correspond with the marginal spawning of fish transiting from the Gulf of Mexico to the northern feeding grounds or of fish that interrupted their southward migration. In addition to the documented behaviour of ABFT natal homing for spawning (Rooker et al., 2008a,b), the spatial distribution of predicted habitats delineated a continuum of favourable feeding environments suitable for the migration of the eastern stock from the North Atlantic to the Mediterranean Sea while it described a discontinuity in the area off Florida and the Bahamas (notably in relation to positive levels of SSHa; see the SSHa geographic distribution in Fig. SI-2) that may prevent part of the western large ABFT from moving towards the Gulf of Mexico. More sampling and tagging experiments off Florida and the Bahamas in spring would allow us to investigate the behaviour of large ABFT and the environmental suitability of the area for migration and spawning. We question whether a discontinuity of favourable environment may alter the reproduction success of the western stock (e.g. excessive energetic cost, unsuitable conditions for larvae survival) or if the apparent east-west stock differences for reproduction (age at first maturity, presence in the Gulf of Mexico of giant ABFT only) result from unsuitable field observations (untargeted fishing effort and lack of histological studies).

Use in management
This habitat analysis encompassed most of the known spatial domain of ABFT and the entire area currently subject to management. The ecological niche modelling approach involved a relatively fine spatial resolution of the selected oceanic processes, considered the main size classes of the species and both spawning and feeding behaviours. As such, the results help to explain most of the variability in the temporal and spatial distribution of the population, provided that migration patterns are considered. The ENM could therefore be used to develop a standardized CPUE index that would account for changes in habitat (e.g. Bigelow and Maunder, 2007). Second, the ENM could be used within a spatially explicit stock assessment model, e.g. to help quantify the parameterization of movement between different spatial units (Kerr et al., 2013;Taylor et al., 2011). The improved understanding of ABFT behaviour and movements should also help in identifying key habitat for protection. Such a model may be used as a tool to monitor the variability of ABFT habitats, notably after climate change, which is key information for adaptive management. Two fishing closure areas (see Fig. 6a) were implemented in 2015 in the northern Gulf of Mexico in April-May to limit ABFT bycatch (NOAA, 2014). These two areas are however small compared to the size of the potential spawning grounds, they should also be applied during winter months and do not seem to accommodate the inter-annual spatial dynamics of ABFT presence (Wilson et al., 2015). A dynamic, real-time approach is probably needed in which the ENM could help to define both habitats of ABFT and targeted species (e.g. yellowfin tuna) and inform the longline fleet, similar to what was implemented off Eastern Australia (Hartog et al., 2011). More generally, this habitat approach should allow for input controls of fisheries, such as dynamic spatial management (Eveson et al., 2015;Lewison et al., 2015;Maxwell et al., 2015), especially in a framework where stakeholders, fishers and scientists actively cooperate to preserve the resource and reduce costs. The impetus for detailed spatial information will nevertheless increase with the integration of fisheries in spatial planning principles and ecosystem based approaches to management that are being implemented in national and international policies (European Union, 2014;Link and Browman, 2014;Soma et al., 2015).