Research ArticleECOLOGY

Atlantic Multidecadal Oscillations drive the basin-scale distribution of Atlantic bluefin tuna

See allHide authors and affiliations

Science Advances  02 Jan 2019:
Vol. 5, no. 1, eaar6993
DOI: 10.1126/sciadv.aar6993


The Atlantic bluefin tuna (hereafter referred to as “bluefin tuna”), one of the world’s most valuable and exploited fish species, has been declining in abundance throughout the Atlantic from the 1960s until the mid-2000s. Following the establishment of drastic management measures, the stock has started to recover recently and, as a result, stakeholders have raised catch quotas by 50% for the period 2017–2020. However, stock assessments still omit the natural, long-term variability in the species distribution. Here, we explore the century-scale fluctuations in bluefin tuna abundance and distribution to demonstrate a prevailing influence of the Atlantic Multidecadal Oscillation (AMO) to provide new insights into both the collapse of the Nordic bluefin tuna fishery circa 1963 and the recent increase in bluefin tuna abundance in the Northeast Atlantic. Our results demonstrate how climatic variability can modulate the distribution of a large migrating species to generate rapid changes in its regional abundance, and we argue that climatic variability must not be overlooked in stock management plans for effective conservation.


The Atlantic bluefin tuna (Thunnus thynnus, Linnaeus 1758) is a long-lived, widespread migrating species with the largest thermal tolerance among tunas (1); it is also one of the world’s most commercially exploited marine fishes (2). The bluefin tuna, currently managed as two separated stocks in the Atlantic (3), has been declining in abundance until the mid-2000s (4) and is listed as endangered on the International Union for Conservation of Nature (IUCN) Red List. Although electronic tagging programs begun in the late 1990s have provided new insights on the species’ migratory behavior to show that eastern (western) bluefin tuna may migrate to the western (eastern) Atlantic and remain there for a few months to a few years (3), little is known about the long-term variability in the migratory behavior, regional abundance, and the spatial distribution of bluefin tuna (5).

Bluefin tuna has shown centennial periodic fluctuations (up to threefold) in abundance in the Mediterranean Sea (6), and the extremely rapid 1960s decline of the Nordic fishery remains one of the world’s most spectacular fisheries’ collapses (7). While the Nordic fishery collapse is attributed to overfishing without clear evidence of hydroclimatic influences so far (8), environmental change could have played a role. Recurrent observations of bluefin tuna after the late 1990s suggest that bluefin tuna is returning to Nordic feeding grounds after three decades of depleted abundances (9). Since stock assessments have not suggested any recovery before the mid-2000s (4), it is currently unclear why bluefin tuna has reappeared in the Nordic region. At a time when the recent increase in bluefin tuna catch quotas will exacerbate the pressure on an IUCN Red List species (5), we have investigated the basin-scale variability in the distribution of bluefin tuna. Our goal is to determine whether the current return of bluefin tuna in the Nordic region can be explained by large-scale, hydroclimatic variability and, if so, whether the environment may explain some of the historical changes in the bluefin tuna’s abundance and distribution.


The boosted regression tree (BRT) models used to examine the relationships between long-term variability in eastern bluefin tuna abundance and hydroclimatic variability (Materials and Methods and fig. S1) accurately reconstructed the century-scale fluctuations seen in the eastern bluefin tuna abundance index (Fig. 1). While several studies have previously explored the impact of the North Atlantic Oscillation (NAO), Northern Hemisphere temperature (NHT) anomalies, and total solar irradiance (TSI) on bluefin tuna abundance (10, 11), we found that most of the fluctuations are related to the Atlantic Multidecadal Oscillation (AMO) with a relative influence of 47.1% (Fig. 1B). Higher values of the eastern bluefin tuna abundance index occurred during positive (warm) AMO phases, whereas lower values were related to negative (cold) phases (fig. S2). The AMO still remained the most important hydroclimatic factor for predicting recruitment when we introduced a one-generation lag (i.e., 16 years; Materials and Methods), although its relative importance was slightly reduced (36%; fig. S2). We observed that the abundance of bluefin tuna in traps sometimes preceded shifts in AMO phase, e.g., during the Maunder minimum and at the end of the Dalton minimum, which suggests that the link between abundance and the environment may be more complex than simply through recruitment. For example, bluefin tuna can swim large distances quickly, and so, they are able to track environmental changes rapidly to follow the most favorable areas in the northeast (NE) Atlantic (3), where productivity is mainly driven by the AMO (12). So, while the relationship between the AMO and bluefin tuna may involve both recruitment and older stages, our results reveal that long-term changes in the AMO state can explain the periodic bluefin tuna fluctuations in the eastern Atlantic (Fig. 1 and fig. S3).

Fig. 1 Ensemble reconstructions of long-term fluctuations in Atlantic bluefin tuna abundance.

(A and C) Historical records (black line) and predicted long-term fluctuations in abundance (averaged prediction as an orange line with 5 and 95% confidence intervals as gray shading) calculated from (A) 1634–1929 (preindustrial tuna fishery period) with no lag and (C) with a one-generation lag (16 years). (B and D) Mean relative influence [and associated SD (standard deviation)] of the four hydroclimatic variables used to reconstruct bluefin tuna abundance.

We found that the NAO was the second or third most important hydroclimatic variable influencing bluefin tuna (relative influences of 19.5% without a lag and 21.5% with a one-generation lag; Fig. 1B and fig. S2). Higher values of the adult abundance index occur during negative and positive NAO phases and lower abundances at intermediate values, while recruitment tends to increase with the NAO index (fig. S2); this observation may explain the contrasting results found previously for the influence of the NAO on bluefin tuna recruitment and abundance (10, 13, 14). Such a pattern could arise from the NAO’s multivariate effects on biological communities (15), which may influence both the bluefin tuna survival rates and feeding grounds’ productivity. We detected a more negative impact of high NHT on recruitment than on adult abundance, although trends were comparable (fig. S2), and this agrees with studies showing that eastern bluefin tuna abundance is negatively correlated with both Mediterranean Sea and global sea temperatures, which could be explained by long-term changes in migratory behavior (10, 11).

We found that the TSI, which has been negatively correlated with abundance (11), had the smallest influence on both adult bluefin tuna abundance and recruitment (Fig. 1, B and D). The effect of TSI was limited to very low irradiance periods, i.e., bluefin tuna abundance was systematically low when the TSI was below 1360.3 for adult abundance and 1360.1 for recruitment; no changes were seen when the TSI exceeded these thresholds (fig. S2), which may also explain why the relationship vanished in the 20th century (11) after irradiance increased at the end of the Little Ice Age (fig. S3).

While the effect of TSI and NHT were similar with and without a one-generation lag, the AMO and NAO showed a different effect at their extreme values (fig. S2). The most negative NAO phases had a negative impact on bluefin tuna recruitment (i.e., lower abundance 16 years later), and the most positive AMO phases were associated with low bluefin tuna recruitment. Nevertheless, the period of highest AMO and low bluefin tuna abundance only occurred in the 1700s to the 1720s, at the very end of the Maunder Minimum (fig. S3), and corresponded to the combination of the coldest NHT and the lowest irradiance levels. As highlighted by the BRT models, the irradiance effect on bluefin tuna abundance displayed a marked threshold and so that the NHT and TSI may have caused this sudden drop in bluefin tuna abundance regardless of the AMO phase.

We next investigated the long-term (1891–2011) spatiotemporal changes in habitat suitability of bluefin tuna in the Atlantic using ecological niche modeling (ENM) (fig. S1). While studies have shown changes in the habitat distribution and seasonal size-dependent feeding grounds (16), none has assessed the link with hydroclimatic variability. This procedure was repeated 30 times per model to provide the mean values and standard deviations (SDs) of the continuous Boyce index (CBI). The best model (CBI = 0.84; table S2) was obtained when each occurrence of bluefin tuna was associated to a unique triplet of environmental parameters; this configuration reproduced the overall bluefin tuna distribution well and was therefore used for subsequent analyses.

We applied a principal components analysis (PCA) on the annual anomalies of probabilities of occurrence of bluefin tuna (hereafter referred to as habitat suitability) modeled from 1891 to 2011. The spatial pattern of the first normalized eigenvector, which showed the highest significance (37.5% of the total variance), revealed an opposition between the northeastern and southwestern North Atlantic (fig. S4). Associated long-term changes (1891–2011) in the first principal component (PC) of habitat suitability correlated strongly with the abundance index of eastern bluefin tuna (r = 0.72, P < 0.001; fig. S4), and the link with historical fluctuations was even higher when reconstructed with PC1, PC2, and PC4 (r = 0.82, P < 0.001; Fig. 2). Our analysis therefore suggests that bluefin tuna occurrence in the North Atlantic is controlled by a northeastern/southwestern “seesaw” of habitat suitability so that long-term local fluctuations in bluefin tuna abundance may reflect changes in spatial distribution rather than changes in the size of eastern and western populations.

Fig. 2 Habitat-based reconstruction of Atlantic bluefin tuna abundance.

Coefficient of linear correlation (r) and its associated probability (P) between historical records of Atlantic bluefin tuna abundance (from Fig. 1; black line) and long-term changes (1891–2011) in the species’ habitat suitability (orange line) reconstructed from the PCA computed on annual anomalies of probability of occurrence of bluefin tuna.

Spatial patterns of habitat suitability of bluefin tuna in the North Atlantic were assessed for positive (1929–1962 and 1995 to present) and negative (1896–1928 and 1963–1994) AMO phases (Fig. 3). This analysis showed that high records of bluefin tuna in the NE Atlantic coincided with high positive habitat suitability observed during positive AMO phases (Fig. 3A), while lower records occurred during negative AMO phases when habitat suitability became negative (Fig. 3B). Our analysis therefore highlights the AMO’s controlling role on bluefin tuna’s basin-scale distribution. When calculating the average habitat suitability of bluefin tuna around United Kingdom (Fig. 3B), we found that anomalies follow the main historical periods of bluefin tuna presence and absence: negative in the entire region until the mid-1920s, positive until the mid-1960s, negative until the mid-1990s, and returning to positive until today (Fig. 3C and fig. S5). Our analysis suggests that bluefin tuna persisted around the United Kingdom only during the two positive AMO phases from 1891 onward (fig. S3), which corresponded to periods of highest anomalies in habitat suitability (Fig. 3C). The spatial distribution of bluefin tuna modeled from the species’ ecological niche therefore captured the rise and fall of the Nordic fishery between the late 1920s and 1960s and its return after the late 1990s. The collapse of the Nordic bluefin tuna fishery circa 1963 coincided with an extremely rapid 2-year transition (1962–1963) from a highly positive AMO phase to its lowest recorded value (Figs. 3C and 4). The AMO reversed to a positive phase in 1996, and a fishery targeting mature fish was established around Iceland in 1997 when bluefin tuna returned to the region (17). Schools of bluefin tuna have also been seen around the United Kingdom since 1996 (9). Together, our observations suggest that the hydroclimatic variability influences bluefin tuna distribution strongly, with higher (warm AMO) and lower (cold AMO) utilization rates in the northern and southern regions of the North Atlantic, respectively. Although no time series is sufficiently long to validate the southwest/NE changes in distribution of the western stock (Materials and Methods), both stocks may follow the same pattern of variability as they have the same ecological niche.

Fig. 3 AMO phases and spatiotemporal variability in habitat suitability and distribution of Atlantic bluefin tuna.

(A and B) Anomalies of habitat suitability, reported occurrences (per 1° by 1° geographical cell; red dots), and mass centroids of occurrences (black dots) during (A) positive (1929–1962 and 1995 to present; nbins = 964) and (B) negative (1896–1928 and 1963–1994; nbins = 979) AMO phases. (C) Time series of mean anomalies (blue line) in the Nordic region [inset in (B)] from 1891 to 2011. Vertical dashed lines indicate abrupt changes in habitat suitability, and horizontal lines indicate the mean habitat suitability for each phase. The size of the tuna is proportional to its frequency of occurrence.

Fig. 4 Habitat suitability of Atlantic bluefin tuna when the Nordic fishery collapsed.

Mean anomalies of probability of occurrence of Atlantic bluefin tuna before the collapse (1959–1961; left) and when the fishery collapsed (1962–1963; right).

Long-term climate variability is known to influence fisheries production at various spatiotemporal scales (18). Here, we have shown that the AMO is an important determinant of the bluefin tuna’s spatial distribution and regional abundance in the North Atlantic and that the recent recovery of the eastern stock (4) may reflect the current positive AMO phase. Our analyses of environmental effects on bluefin tuna can therefore provide an insight to help interpret (i) the collapse of bluefin tuna in 1963 and (ii) the current increase of bluefin tuna in Nordic seas.

It is assumed that eastern bluefin tuna spawn exclusively in the Mediterranean Sea (2), where the influence of the AMO is lower than elsewhere (19). We tested the influence of hydroclimatic variability on recruitment and highlighted, in line with previous studies (10, 11), that high NHT negatively affects both the adult abundance in traps and the recruitment (fig. S2). However, the Mediterranean sea surface temperature (SST) [averaged from the Centennial Observation-Based Estimates (COBE) dataset; see Materials and Methods] and the recruitment year classes (4) were significantly and positively correlated over the period 1968–2014 (r = 0.59, P < 0.001), with the strongest correlations in summer (June, July, and August, i.e., immediately after the main spawning period), while SSB followed an opposite pattern (r = −0.63, P < 0.001). In response to rising temperatures, the average habitat suitability of bluefin tuna has also been decreasing in the Atlantic and Mediterranean Sea (correlated negatively with NHT; r = −0.49, P < 0.001). This suggests that, while higher local SST may enhance recruitment, it would still have a deleterious effect on adult bluefin tuna. If warming continues, then bluefin tuna may indeed become constrained by the upper thermal limit of its spawning preference in the Mediterranean Sea (fig. S6) (20) and could begin to use different regions in the Atlantic. This consequence of warming has been predicted in the Gulf of Mexico (21), and bluefin tuna have been discovered to spawn along the NE coast of the United States, which may constitute a recent expansion of spawning habitat (2224). Consequently, future warming may also alter ecological barriers (25) and counterbalance the influence of future cold AMO phases on the species’ distribution, causing bluefin tuna to persist in Nordic seas.

To manage Atlantic bluefin tuna and other highly migrating species sustainably, long-term recovery plans should therefore be addressed at a basin scale by combining the effects of hydroclimatic variability, such as the AMO, global climate change and fishing, and local increases in abundance should not be used as a reason to relax quotas for commercial or recreational fisheries without having done so. Since large-scale changes in the abundance of top predators integrate changes at lower trophic levels (26), regional fluctuations in bluefin tuna abundance may also serve to indicate even wider ecological changes, providing useful insights for further research.


Biological data

Occurrence data of the Atlantic bluefin tuna (T. thynnus, Linnaeus 1758). The information on Atlantic bluefin tuna distribution was compiled from georeferenced (i.e., longitude and latitude) and time-referenced data (year, month, and day) available in public databases provided by the Ocean Biogeographic Information System (n = 9288;, the Global Biodiversity Information Facility (n = 597;, the International Commission for the Conservation of Atlantic Tunas (ICCAT) Task I (n = 165,050), and the ICCAT conventional tagging (n = 87,466; Duplicates were discarded on the basis of year, month, day, longitude, and latitude. A checking procedure was applied on ICCAT data to remove occurrences with a spatial resolution lower than 1° longitude by 1° latitude and a temporal resolution lower than 1 month (i.e., the resolution of the gridded climatic data described in the “Environmental data” section). Occurrences located outside the Atlantic Ocean or south of 12°S (i.e., outside the natural distribution range of Atlantic bluefin tuna) (4) were also excluded. A total of 27,617 occurrences were retained to estimate the bluefin tuna’s habitat suitability. Distribution maps of occurrence records are presented in Fig. 3 (red dots).

Elaboration of the Atlantic bluefin tuna abundance index. The Atlantic bluefin tuna is currently managed by the ICCAT as western and eastern stocks (4), although these two stocks undertake annual migrations between the main spawning (i.e., the Gulf of Mexico and the Mediterranean Sea) and feeding (northwest, central, and NE Atlantic) regions (27). While current changes in abundance occurring in the eastern region may relate to ~90% of Atlantic bluefin tuna populations (4), the two stocks were assumed to be approximately the same order of magnitude in the 1970s until the western stock collapsed (4). However, the actual long-term fluctuations in the Atlantic-wide distribution of Atlantic bluefin tuna remain unknown. We focused our analyses on eastern bluefin tuna because historical data are available since the 16th century for this stock only (4, 6), assuming that it represented the baseline variability of about half of bluefin tuna populations over the past four centuries.

To assess the relative fluctuations in eastern Atlantic bluefin tuna from the mid-17th century onward, we used an index of abundance based on two databases: (i) the historical catches in traps located in and around the Mediterranean Sea during the period 1634–1959 that targeted large, mature individuals during their spawning migrations (6) and (ii) the last ICCAT spawning stock biomass (SSB) estimation for the period 1950–2011 (17). While the consideration of young stages in the SSB might have slightly diminished the explanatory power of analyses after the 1960s, both datasets remained the most relevant to assess the major long-term changes in eastern bluefin tuna abundance (6, 28). We scaled the SSB index to the historical index over the common period 1950–1959 and evaluated their adequacy by applying a Wilcoxon test, which revealed no significant difference (W = 23, P = 0.14). Both indices were then combined in a single abundance index spanning from 1634 to 2011, with relative abundance being averaged for the common period.

Environmental data

Environmental data used to perform the ENM. Monthly mean SSTs from 1891 onward were obtained from the COBE of SST version 2 (COBE-SST2; ( (29). Sea surface salinity (SSS) was obtained from the World Ocean Database 2013 ( We considered that salinity was constant over time in the model because the spatial variance of SSS greatly exceeds its temporal variance (30). The bottom topography was extracted from the General Bathymetric Chart of the Ocean ( All databases are available with a resolution of 1° by 1°.

Large-scale hydroclimatic indices. We investigated the relative influence of four large-scale hydroclimatic indices used in the reconstruction of the Atlantic bluefin tuna abundance index for the period 1634–2011.

1) NHT anomalies were used as a proxy of the potential effect of climate change in the Northern Hemisphere. We used the dataset developed by the Tree-ring Network for the period 918–2011 (31) (

2) The TSI index (the total spectrally integrated energy input to the top of the Earth’s atmosphere at a standard distance of one astronomical unit from the Sun) is assumed to be a good proxy of particularly cold periods (32) around the end of the Little Ice Age (14th to mid-19th century; fig. S3). Data for the period 1610–2011 were provided by the National Oceanic and Atmospheric Administration (

3) The winter NAO index describes the basin-scale gradient of atmospheric pressures over the North Atlantic between the high pressures centered on the subtropical Atlantic and the low pressures over Iceland in winter (33). We used (i) a millennium-long (1049–1969) tree-ring–based reconstruction of the NAO ( and (ii) the NAO index based on a rotated PCA performed on monthly standardized 500-mbar height anomalies over the North Atlantic sector for the period 1950–2011 (

4) The AMO index characterizes the multidecadal ocean/atmosphere natural variability in temperatures, in a range of 0.4°C, in many oceanic regions of the North Atlantic, with a periodicity ranging from 60 to 100 years (34). Here, we combined (i) a tree-ring–based reconstruction of the AMO for the period 1567–1990 ( and (ii) the index constructed from extended reconstructed SST data for the period 1856 onward and averaged in the area of 25° to 60°N and 7° to 75°W, minus regression on global mean temperature (

We combined two hydroclimatic versions of the AMO and NAO indices to cover the whole time period 1634–2011. For each version, the oldest was rescaled (range and offset) to be comparable with the recent one. Variances were then homogenized by applying a moving mean of the larger window applied in the index reconstruction [i.e., 10 years for the AMO (34) and 3 years for both NAO indices. For each index, we calculated and tested the correlation between the reconstructions based on tree-ring data and contemporary indices, which proved to be highly significant in each case (for AMO indices: F1,104 = 268.6, R2adj = 0.72, P < 0.001; for NAO indices: F1,15 = 70.0, R2adj = 0.81, P < 0.001); old and recent indices were therefore averaged per year over common periods and combined into a single index.

Time series analyses

BRT and selection of critical parameters. Long-term fluctuations in eastern bluefin tuna abundance and in hydroclimatic variability were analyzed with BRT (fig. S1, analysis 1) (35). BRT models are an ensemble method for fitting statistical models that combine the strengths of regression trees and boosting, which enables the identification of nonlinear complex relationships and interactions between variables without the need for previous data transformation or any underlying assumptions (35). The analyses were computed using R software and the gbm and dismo packages. Critical parameter values were first tested in an exploratory phase. Large numbers of trees (>1000) are required to build reliable models (35); therefore, the learning rate (lr; the contribution of each tree to the growing model) was set at a low value (lr = 0.001) to generate systematically more than 3000 trees. Using the gbm.step function, the exact number of trees was determined with a 10-fold cross-validation while increasing the model complexity. A total of 24 sets of models were built to test the two other critical parameters, tree complexity (tr; the number of nodes in a tree) and bagging fraction (bag; the proportion of data to be selected at each step). Tested values ranged from 1 to 12 for tr and from 0.4 to 0.95 for bag (randomness is removed when bag is set to 1). Model performance was assessed with the training data correlation (equivalent to R2adj in linear regression models) and the 10-fold cross-validation correlation mean deviance and coefficient (35). The sets of models with tr = 9 and bag = 0.8 systematically led to both the highest R2adj and cross-validation correlation coefficient and the lowest cross-validation deviance (table S1), and so, this configuration was used to build the two final models.

BRT models construction and performance. One hundred bootstraps, each with 10-fold cross-validation, were applied to compute the mean, SD, and confidence intervals (36) of (i) the mean reconstructed abundance of eastern bluefin tuna and (ii) the relative influence of each of the four hydroclimatic indices. Two BRT models were used to reconstruct bluefin tuna abundance for years 1634–1929, i.e., before the first industrial fishery and rapid north Atlantic warming (37), to assess the patterns of abundance while omitting the main components of human disturbances (i.e., the overfishing mortality and the rapid elevation of temperature with the increasing emissions of greenhouse gases). First, we considered the index of bluefin tuna abundance without a lag, and then, we incorporated a one-generation lag (i.e., 16-year lag selected by cross-correlation and taking into account the temporal autocorrelation of the time series). This enabled us to assess both the direct effects of hydroclimatic variability on adult distribution and the indirect effects on recruitment occurring one generation later. The performance of the two models in reconstructing eastern bluefin tuna abundances over 1634–2011 was assessed with the training data correlation (R2adj) (fig. S2).

Our BRT models reproduced well the temporal autocorrelation that was detected in the observed bluefin tuna abundance index. It slowly decreases until a 27-year lag in the historical time series and to a 29-year lag in the predicted time series, which is in line with the maximum life-span of bluefin tuna.

Species distribution analyses

Ecological niche model: NPPEN. The nonparametric probabilistic ecological niche (NPPEN) model is an ENM specifically designed to be applied on presence-only data (38). This technique allows the modeling of the ecological niche of a species and the mapping of its spatial distribution by calculating its probabilities of occurrence. The first step consisted in constructing a reference matrix with environmental data corresponding to the presence records. The reference matrix was homogenized to remove as much as possible the inaccurate reporting of occurrence records, and this resulted in a multidimensional matrix, with each dimension reflecting an environmental factor (38). In the second step, the Mahalanobis generalized distance was calculated between the observations and the homogenized reference matrix. In the third step, the model calculated the probability of each grid point to belong to the reference matrix (38). Last, the probabilities of species occurrence were mapped in the geographical space. Here, the NPPEN model was applied (i) to calculate the ecological niche (sensu Hutchinson) of Atlantic bluefin tuna and (ii) to project its spatial distribution in the North Atlantic (north of 12°S, i.e., its distribution range) over the period 1891–2011 (fig. S1, analysis 2).

Model validation. The adequacy between observed and modeled spatial distributions was assessed with the CBI (39). The CBI, developed for presence-only data and insensitive to pseudo-absences, was recommended for accurate evaluation of presence-only models such as ours (39). A model was considered “wrong” when CBI values were below −0.5, “average to random” for values ranging from −0.5 to 0.5, and “good” for values above 0.5 (i.e., occurrences in regions with high habitat suitability). To investigate whether the random selection of data influenced the modeling of the ecological niche, we randomly selected 70% of the observed data to build the ENM, and the other 30% was used to evaluate its accuracy. This procedure was repeated 30 times for each simulation to provide the mean and SD of metric values. For each month and each run, the CBI was computed and averaged per simulation.

Selection of environmental variables. Seven models were performed to determine the best combination of environmental variables (i.e., the three variables that contributed most to the spatial distribution of Atlantic bluefin tuna; table S2). The model that considered each occurrence of bluefin tuna associated to a unique triplet of environmental parameters (SST, SSS, and bathymetry) displayed the highest adequacy with the observed spatial distribution (CBI = 0.84; table S2); we therefore retained this optimal configuration to assess the spatiotemporal variability in the habitat suitability of Atlantic bluefin tuna over the period 1891–2011.

Estimation of the potential influence of localized intensive sampling or duplicates. Fisheries data are typically biased toward undersampling in remote areas and oversampling in known fishing grounds (40). To consider this potential bias, we tested 5 levels of homogenization of the reference matrix (see the “Ecological niche model: NPPEN” section) by using an increasing number of occurrences (from 1 to 3, 5, and 10) to evaluate the representativeness of each set of environmental conditions (i.e., SST, SSS, and/or bathymetry) on the model accuracy. The model that considered one species occurrence per triplet of environmental conditions outperformed other models (highest mean values and lowest SDs of the evaluation metric; table S2) and was selected for subsequent analyses.

Long-term changes in habitat suitability of Atlantic bluefin tuna. We then investigated long-term changes in anomalies of probability of occurrence of bluefin tuna (also termed habitat suitability) obtained from the NPPEN model (fig. S1). We first calculated monthly anomalies for each geographical cell of the spatial domain, i.e., the difference between each probability of occurrence of bluefin tuna from 1891 to 2011 and the climatology of probabilities of occurrence of bluefin tuna calculated for the period 1891–2011. We then performed a standardized PCA to examine long-term changes in habitat suitability in the distribution range of bluefin tuna (i.e., the Atlantic Ocean and Mediterranean Sea from 100°W to 37°E and from 12°S to 75°N). The first four PCs (55% of the total variance) reflected the main long-term changes in habitat suitability, and the normalized eigenvectors represented the correlation between changes in habitat suitability in each geographical cell and the corresponding PCs (fig. S4).

We also evaluated long-term changes in habitat suitability by averaging monthly anomalies of probability of occurrence of bluefin tuna for different time periods: (i) negative (climatology from 1896 to 1928 and from 1963 to 1994) and positive (climatology from 1929 to 1962 and from 1995 to 2011) AMO phases (Fig. 3), (ii) before and after the Nordic fisheries collapse circa 1963, (iii) per decades from 1896 to 2011 (fig. S5), and (iv) at an annual scale from 1896 to 2011 in the Nordic region (from 30°W to 12°E and 48°N to 70°N; Fig. 3C). Abrupt changes in the Nordic region were then investigated by applying the method of optimal detection of change points with a linear computational cost (41). Mean change in habitat suitability was also assessed at an annual scale in the whole species’ distribution range (Atlantic Ocean and Mediterranean Sea from 100°W to 37°E and from 12°S to 75°N), and the relationship with NHT was investigated with linear correlation.

Relationships between long-term changes in habitat suitability and abundance

We calculated the partial correlations between the observed Atlantic bluefin tuna abundance index and the first four PCs obtained from the PCA applied on anomalies of probability of occurrence of bluefin tuna (fig. S4). After application of a forward/backward procedure, the first, second, and fourth PCs were used to elaborate a habitat suitability index of the Atlantic bluefin tuna (Fig. 2). The relationship between changes in our habitat suitability index and the observed Atlantic bluefin tuna abundance index was then investigated by a linear correlation analysis (Fig. 2).


Supplementary material for this article is available at

Fig. S1. Summary of the two-step analysis conducted in this study.

Fig. S2. Relationships between each hydroclimatic index and the abundance of Atlantic bluefin tuna from the BRT analysis.

Fig. S3. Long-term changes in the Atlantic Multidecadal Oscillation and in the abundance of eastern Atlantic bluefin tuna.

Fig. S4. Spatiotemporal variability in the habitat suitability of Atlantic bluefin tuna.

Fig. S5. Decadal monthly anomalies of habitat suitability of Atlantic bluefin tuna in the NE Atlantic.

Fig. S6. Thermal niche of Atlantic bluefin tuna.

Table S1. Influence of tree complexity and bagging fraction on the performance of the BRT models.

Table S2. Effects of different combinations of environmental variables and of different thresholds of occurrence on the performance of the NPPEN model applied on Atlantic bluefin tuna.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank the reviewers for their valuable comments. Funding: The authors acknowledge the Contrat de plan État-région (CPER) CLIMIBIO program for funding R.F. and the regional program INDICOP for funding E.G. Author contributions: R.F., G.B., and R.R.K. designed the study. R.F. and E.G. performed statistical analyses. R.F. wrote the manuscript and produced display items. G.B., E.G., and R.R.K. discussed the results and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article