Comparing the distribution of tropical tuna associated with drifting fish aggregating devices (DFADs) resulting from catch dependent and independent data

Journal Pre-proof Comparing the distribution of tropical tuna associated with drifting fish aggregating devices (DFADs) resulting from catch dependent and independent data Blanca Orue, Jon Lopez, Maria Grazia Pennino, Gala Moreno, Josu Santiago, Hilario Murua PII:

S0967-0645(19)30069-4

DOI:

https://doi.org/10.1016/j.dsr2.2020.104747

Reference:

DSRII 104747

To appear in:

Deep-Sea Research Part II

Received Date: 21 February 2019 Revised Date:

31 January 2020

Accepted Date: 31 January 2020

Please cite this article as: Orue, B., Lopez, J., Pennino, M.G., Moreno, G., Santiago, J., Murua, H., Comparing the distribution of tropical tuna associated with drifting fish aggregating devices (DFADs) resulting from catch dependent and independent data, Deep-Sea Research Part II (2020), doi: https:// doi.org/10.1016/j.dsr2.2020.104747. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Blanca Orue: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Roles/Writing – original draft, Writing – review & editing. Jon Lopez: Conceptualization, Investigation, Methodology, Validation, Writing – review & editing. Maria Grazia Pennino: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Writing – review & editing. Gala Moreno: Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing – review & editing. Josu Santiago: Funding acquisition, Investigation, Project administration, Validation, Writing – review & editing. Hilario Murua: Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing – review & editing

1

Comparing the distribution of tropical tuna associated with drifting fish

2

aggregating devices (DFADs) resulting from catch dependent and independent

3

data

4

Blanca Oruea,*, Jon Lopez b, Maria Grazia Penninoc, Gala Morenod, Josu Santiagoa, Hilario Muruaa-d

5 6 7 8 9 10

a

11

*corresponding author.

12

E-mail address: [email protected] (B.Orue)

AZTI-Tecnalia, Herrera kaia portualdea z/g 20110 Pasaia (Gipuzkoa), Spain

b

Inter-American Tropical Tuna Commission 8901 La Jolla Shores Drive, La Jolla CA 92037-1509

c

Instituto Español de Oceanografía, Centro Oceanográfico de Vigo. Subida a Radio Faro, 50-52 36390 Vigo, Spain

d

International Seafood Sustainability Foundation (ISSF), 1440 G Street NW Washington DC 20005

13 14 15

Keywords: DFAD, Indian Ocean, INLA, Species distribution model, Tuna

16 17 18

ABSTRACT

19 20

Species distribution models (SDMs) are used for a variety of scientific and management applications.

21

For species associated with drifting fish aggregating devices (DFADs), such as tuna, spatial models

22

can help tuna Regional Fisheries Management Organizations (t-RFMOs) understand their habitat

23

characteristics and dynamics. DFADs are monitored and tracked with satellite linked echo-sounder

24

buoys, which remotely provide fishers rough estimates of the abundance of fish underneath them.

25

Although this type of catch-independent data has been recently used in scientific studies, SDMs

26

using these data have never been compared with models using catch-dependent data (i.e. nominal

27

catch data). This study investigates the results obtained with both data sources using Bayesian

28

Hierarchical spatio-temporal models, allowing to analyze their advantages and disadvantages, as well

29

as compare the predicted distributions. Although the two model outputs show, in general, similar

30

areas of tuna presence under the DFADs, the most remarkable result of the comparison between the

31

models derived from the two different data sources is the precision of the hotspots identified in the

32

prediction maps. The maps obtained with acoustic data allow identifying areas of high probability of

33

tuna presence under the DFADs with greater precision, whereas the maps derived from catch data do

34

not allow observing any variation on a finer scale. The application of spatio-temporal models of tuna

35

associated with DFADs using acoustic data provided by fishers´ echo-sounder buoys appears

36

promising to identify the distribution dynamics of the species in a cost-effective way and may help

37

designing integrated spatial programs for more efficient fishery management.

38 39 40

1. Introduction

41 42

Tropical tuna, such as skipjack (Katsuwonus pelamis), yellowfin (Thunnus albacares) and

43

bigeye (Thunnus obesus), aggregate under floating objects in the Atlantic, Indian and Pacific Oceans

44

(Castro et al., 2002). Man-made floating objects, also called drifting fish aggregating devices

45

(DFADs), have been deployed and used by the tropical tuna purse-seine fishery since the 1990s to

46

increase their catches and fishing efficiency (Fonteneau et al., 2013; Davies et al., 2014). These

47

devices have been extensively used since then, and especially in recent years, it has been shown a

48

shift in strategy towards the use of more DFADs (Lennert-Cody et al., 2018). For example, around

49

80% of the Spanish tuna purse-seine sets in the Indian Ocean in 2017 were made on DFADs (Báez et

50

al., 2018), while the rest of the catch comes from sets on unassociated tuna schools, also called free-

51

swimming schools (FSC). Currently, DFADs are monitored and tracked with satellite linked echo-

2

52

sounder buoys (Lopez et al., 2014), which remotely inform fishers in near real-time about the

53

geolocation of the object and provide rough estimates of the abundance of fish underneath them.

54

Because of the enormous importance and value of this fishery, increasing the knowledge of

55

the spatio-temporal distribution of tuna associated with DFADs may help to develop adequate

56

science-based advice for tuna Regional Fisheries Management Organizations (t-RFMOs). Species

57

distribution models (hereinafter SDMs) are an essential tool for science and management, as they

58

provides a clear picture of the distribution dynamics and extent of marine resources (Pennino et al.,

59

2014; Robinson et al., 2017; Martínez-Minaya et al., 2018). In general, information to build marine

60

SDMs can be obtained from two main sources: catch-independent and catch-dependent data.

61

To date, most studies of the distribution of tuna species are conducted using catch-dependent

62

data, such as catch data coming from fishers logbook (Chen et al., 2005; Lee et al., 2005; Rajapaksha

63

et al., 2013) or observers´ data onboard commercial vessels (Sequeira et al., 2012; Lezama-Ochoa et

64

al., 2016; Coelho et al., 2017). Catch-dependent data, due to its mandatory collection nature, are

65

generally easier to access, less expensive to collect and usually involve a relatively large spatially and

66

temporally sample size as it covers the whole distribution of the fishery. Yet, there are some biases

67

associated with such data. One of the main concerns associated with these data is that target-species´

68

samples are collected by preferential sampling as the fishing fleets are commercially driven (Conn et

69

al., 2017). In addition, they suffer from quality control and assurance issues due to legibility,

70

problems with standardization in data collection, aggregation of species into generic groups

71

transcription errors and misreporting of catches (Lowell et al., 2015; Pennino et al., 2016; Bradley et

72

al., 2019).

73

On the other hand, catch-independent data usually have fewer biases, but one of the major

74

associated problems is that they are often collected at high economic and human cost, resulting in

75

smaller sample sizes covering generally specific areas. In addition, for large migratory pelagic

76

species like tuna, catch-independent surveys covering their whole spatio-temporal distribution are 3

77

very difficult and expensive to carry out due to their enormous geographic range. Recent works have

78

highlighted the importance that fishers’ echo-sounder buoys may have to investigate several issues of

79

scientific relevance (Lopez et al., 2016; Moreno et al., 2016). Although these data are recorded by

80

tools that are not originally developed for scientific analyses, they present many notable advantages.

81

They are continuously recording and collecting geolocated acoustic data and, hence, have the

82

potential to stream huge amounts of information in a very cost-effective manner. Indeed, information

83

from these devices have already been used in recent years to investigate tuna dynamics at different

84

temporal and spatial scales (Lopez et al., 2010; Robert et al., 2013; Lopez et al., 2017a; Lopez et al.,

85

2017b; Orue et al., 2019b).

86

This study compares the resulting spatio-temporal distribution of tuna in the Western Indian

87

Ocean by implementing Bayesian Hierarchical spatial models (hereinafter B-HSM) using both catch-

88

dependent (i.e. nominal catch data) and independent (i.e. acoustic data from fishers´ echo-sounder

89

buoys) data. The results obtained will allow analyzing the advantages and disadvantages of both data

90

sources, as well as determining the suitability of using fishers´ echo-sounders buoys acoustic data for

91

species distribution models, and more generally, the validity of the same data for scientific purposes.

92

Moreover, results of this exercise will contribute to increase the knowledge of the spatial dynamics

93

and distribution of tuna in the Western Indian Ocean and assist the scientific-advice to t-RFMOs.

94 95 96

2. Materials and methods

97 98

2.1. Study area

99

4

100

Our study area is bounded by longitude 30ºE to 80ºE and latitude 15ºN to 30ºS in the Western

101

Indian Ocean, where tropical tuna purse seine fishery operates throughout the year using DFADs

102

(Davies et al., 2014) (Fig. 1). DFADs are not distributed evenly in space and their trajectories are

103

influenced by surface currents and winds (Dagorn et al., 2013; Davies et al., 2014). The ocean

104

circulation in the Indian Ocean is related to the marked monsoon system that greatly affects the

105

oceanography and production of the area (Schott and McCreary, 2001; Wiggert et al., 2006; Schott et

106

al., 2009). The drastic changes in circulation of the surface currents affect biophysical factors (e.g.

107

chlorophyll, temperature, salinity, dissolved oxygen) (Tomczak and Godfrey, 2013) and could affect

108

the presence and relative species composition in an area (Jury et al., 2010). In addition, it is known

109

that the aggregation dynamics of tuna associated with DFADs differed between monsoon periods

110

(Orue et al., 2019b).

111

112

2.2. Catch-independent data

113 114

The acoustic data were collected by Satlink buoys (SATLINK, Madrid, Spain,

115

www.satlink.es) attached to DFADs and were provided by the Spanish purse seine vessel company

116

Echebastar. The echo-sounder transmits to the user the amount of biomass aggregated under each

117

object using an Inmarsat satellite connection and a depth layer echo-integration procedure

118

(Simmonds and MacLennan, 2005). The observation depth range extends from 3 to 115 m, composed

119

of ten homogeneous layers, each with a resolution of 11.2 m (see Lopez et al. (2016) for more

120

technical details on the buoy used) (Fig. 2). A vertical depth limit of 25 m was established as a

121

potential boundary between tuna and non-tuna species, based on experimental evidences from

122

tagging and acoustic surveys around DFADs in the Indian Ocean (Moreno et al., 2007; Govinden et

123

al., 2010; Filmalter et al., 2011; Forget et al., 2015). Although temporary overlap may exist between 5

124

species, it is assumed that, in this region, most of the time tuna is below 25 m while non-tuna species

125

remain in shallower waters. Other studies using the same echo-sounder buoy, addressing similar

126

aggregation dynamics and distribution studies, used the vertical limit of 25 m to separate tuna from

127

non-tuna species (Robert et al., 2013; Lopez et al., 2017a; Lopez et al., 2017b; Orue et al., 2019b).

128

Therefore, if there was an acoustic record between 25 and 115 m, it was tuna presence and, if not,

129

absence.

130

In addition to the acoustic information, buoy data include information about the owner vessel,

131

buoy ID (unique alphanumeric code given by the manufacturer), spatial information (i.e. latitude and

132

longitude), and the date and GMT hour of the sampling. All this information was available for 962

133

buoys with 42,322 acoustic samples, including presence and absence information, which were active

134

in the Indian Ocean from January 2012 to May 2015. Although echo-sounders provide daily data and

135

precise positions, catch-independent information were aggregated at 1ºx1º degrees per month to have

136

the same spatio-temporal resolution as catch-dependent data. All the buoys were attached to newly

137

deployed DFADs (i.e. virgin DFADs; objects with no aggregation associated when first deployed),

138

with a maximum of 60 days at sea (i.e. time constrained based on a preliminary data distribution

139

analysis that showed that after 60 days only 50% of the DFADs were available). Data cleaning was

140

carried out following the protocol proposed by (Orue et al., 2019a). Buoy data and FAD and fishing

141

logbook information were crossed to assure that no fishing activity occurred in our virgin DFAD.

142

Potential fishing sets were identified on the logbooks based on the information of buoy code, position

143

and date. If a fishing set was identified to occur during the DFAD trajectory, the buoy information

144

after the set was eliminated

145 146

2.3. Catch-dependent data

147

6

148

Nominal tropical tuna catch data (i.e. any species; skipjack, yellowfin and bigeye) of the

149

purse-seine fisheries in the Indian Ocean were obtained from the Indian Ocean Tuna Commission

150

(IOTC; www.iotc.org). The dataset includes information on EU DFAD catches over the study period

151

(2012–2015) with a spatial resolution of 1° latitude-longitude grid cell and a monthly temporal

152

resolution. A total of 7357 strata (1ºx1º grid cell and month) with catch information reported by the

153

EU fleets were observed over a period of 4 years (i.e. 2012-2015).

154

Only fishery strata with positive samples (i.e. presences) are available in the original dataset.

155

To generate a presence/absence dataset, pseudo-absences were generated randomly for the study

156

area. The set of pseudo-absences was generated 100 times using the “srswor” function (simple

157

random sampling without replacement) from the “sampling” package (Tille and Matei, 2016) of the

158

R software (R Development Core Team, 2017). In each case, the number of generated pseudo-

159

absences was the same as the number of real presences (Sequeira et al., 2013). Pseudo-absences were

160

then combined with real presences into a single presence-absence dataset to be used in the binomial

161

modeling phase. Although the generation of pseudo-absence datasets is an active research topic and

162

is subject to some debate (e.g. Hastie and Fithian (2013)), we preferred to use such a binomial

163

distribution with a Bayesian spatial model instead of a less accurate model that allows the use of

164

presence-only data (e.g. BIOCLIM, DOMAIN) (Roos et al., 2015).

165 166

2.4. Environmental data

167 168

The environmental variables were obtained from the EU Copernicus Marine Environment

169

Monitoring Service (CMEMS; http://marine.copernicus.eu/). Eight abiotic and biotic variables were

170

downloaded and processed based on previous studies analyzing tropical tuna environmental

171

preferences (Fonteneau et al., 2008; Song et al., 2008; Fraile et al., 2010; Arrizabalaga et al., 2015;

7

172

Lopez et al., 2017b) (Table 1). For the acoustic data, the environmental variables were extracted for

173

each record, according to their position and date. In the case of tuna catch data, as they are monthly

174

data aggregated in 1x1°, we extract the average value of the strata for each variable. All

175

environmental variables were tested for correlation, collinearity, outliers and missing values before

176

their use in the models. Velocity was eliminated as it is highly correlated to eddy kinetic energy

177

(Pearson correlation, r > 0.9; p-value < 0.001). The rest of the variables presented low correlation and

178

collinearity values. Moreover, variables were standardized using the function “decostand” in the

179

“vegan” package (Oksanen et al., 2013) of the R software, in order to facilitate interpretation and to

180

enable comparison of relative weights between variables (Kinas and Andrade, 2017).

181 182

2.5.Species distribution model

183 184

B-HSMs were used to predict the probability of tuna presence with respect to the selected

185

environmental variables using both catch-dependent and independent data. In the case of models

186

performed using catch-independent data, the buoy ID was included as a random effect to account for

187

the potential autocorrelation of the sampling structure, as acoustic data are collected repeatedly in

188

various days by the same buoy for each DFAD. The days spent at sea by the DFAD was also

189

considered as a covariate for the catch-independent model. In both models the response variable was

190

a binary variable that represents the presence (1) or absence (0) of tuna (Yi) in each location i, and

191

then the occurrence is modelled as: ~

=

= ,…., +

+



~ !"# , $# % 8

=

0, ' (, )

2

( =

"+ , ,+

2

) =

"- , ,-

192

where πi represents the probability of the species presence for a given location (i), Xiβ represents the

193

matrix of the fixed effects of the linear predictor, Zi is the buoy random effect and Wi represents the

194

spatially structured random effect at the location i. Prior Gaussian distributions with a zero mean and

195

covariance matrix (Q) was assumed for the spatial component, which depend on the hyperparameters

196

k and τ, and determine the range of the effect and the total variance, respectively. Hyperpriors for k

197

and τ are cantered in values such that the range is about 20% of the diameter of the region and the

198

variance is equal to 1 (Lindgren et al., 2011).

199

We use the Integrated Nested Laplace Approximations (INLA) approach (Rue et al., 2009)

200

and the package INLA (http:\\www.r-inla.org) that is implemented in the R software to obtain

201

Bayesian parameter estimates and predictions. The spatial effects (W) were computed using the

202

Stochastic Partial Differential Equations (SPDE) approach implemented in INLA (Lindgren et al.,

203

2011), which ensures that the continuous spatial domain (also known as Gaussian Random Field;

204

GRF) is discretized into smaller spatial units (known as Gaussian Markov Random Field; GMRF). In

205

such a set-up, the GMRF is reduced to a multivariate Gaussian distribution that is characterized by a

206

zero mean and a spatial dependent covariance matrix that is defined by the Matérn function

207

(Lindgren et al., 2011). Default zero-mean Gaussian prior distributions with a variance of 100 were

208

used for all of the parameters involved in the fixed effects as recommended by Held et al. (2010).

209

Following the Bayesian framework, posterior distributions were obtained for each parameter.

210

This type of distribution, in contrast to the mean and confidence interval produced by classical

211

frequentist analyses, enables explicit probability statements about the parameters. Therefore, the 9

212

region bounded by the 0.025 and 0.975 quantiles of the posterior distribution indicates that the

213

unknown parameter is 95% likely to fall within this range of values (Dell'Apa et al., 2017).

214 215

2.5.1. Model selection

216

The selection of explanatory variables was conducted by comparing all possible interactions,

217

but only the best combination of variables was chosen, based on the Watanabe-Akaike (WAIC)

218

information criterion (Watanabe, 2010) and Log-Conditional Predictive Ordinations (LCPO) (Roos

219

and Held, 2011). Specifically, lower WAIC values indicate a better fit, while lower LCPO scores

220

represent better predictive quality. The best compromise between fit, parsimony and predictive

221

quality is the smaller the WAIC and LCPO values are. Thus, the best models were selected based on

222

the mentioned compromise between low WAIC and LCPO values, containing only relevant

223

predictors (i.e. those with 95% credibility intervals excluding zero (Fonseca et al., 2017)).

224 225

2.5.2. Model prediction

226

Once the inference was carried out, we predicted monthly probability tuna presence in the area

227

of interest using Bayesian kriging (Muñoz et al., 2013). The prediction in INLA was performed

228

simultaneously with the inference, considering the prediction locations as points where the response

229

is missing. The INLA SPDE module allows the construction of a Delaunay triangulation covering the

230

region of interest for the prediction. This has at least three advantages over a regular grid. First, the

231

triangulation is denser in regions where there are more observations and consequently there is more

232

information, and more detail is needed. Second, it saves computing time, because prediction

233

locations are typically much lower in number than those in a regular grid. Thirdly, it is possible to

234

take boundary effects into account by generating a triangulation with small triangles in the domain of

235

interest and using larger triangles in the extension used to avoid boundary effects. Once the

236

prediction is performed in the observed location, there are additional functions that linearly 10

237

interpolate the results within each triangle into a finer regular grid. As a result of the process, a

238

faceted surface prediction is obtained which approximates to the true predictive surface (Muñoz et

239

al., 2013). In order to have the same spatio-temporal resolution, all environmental data were

240

aggregated at 1ºx1º degrees per month, using the ‘raster’ package (Hijmans et al., 2017) of the R

241

software. The two final established models derived from the two different data sources (see section

242

3.1 and 3.2) were used to compute monthly predictions on the monthly means of the selected

243

environmental variables, to take into account possible intraannual variation based on the marked

244

oceanographic dynamics of the area.

245 246

2.5.3. Model validation

247

Model predictions were validated using the leave-group-out cross-validation procedure (Kuhn

248

and Johnson, 2013). Both datasets were randomly split into two main subsets: a training dataset

249

including 80% of the total observations, and a validation dataset containing the remaining 20% of the

250

data (Fielding and Bell, 1997). The training data were used to model the relationship between

251

observed data and the explanatory variables while the second dataset was used to assess the quality of

252

predictions. We repeated the validation procedure 10 times for the best model of each data source and

253

month. Model performance was assessed using the area under the receiver-operating characteristic

254

curve (AUC) (Fielding and Bell, 1997) and the true skill statistic (TSS) (Allouche et al., 2006). AUC

255

measures the ability of the model to discriminate species presence’s locations from those in which

256

the species is absent. AUC values range from 0 to 1, where values closer to 1 indicate better

257

predictions. TSS corrects AUC for the dependence of prevalence on specificity (i.e. ability to

258

correctly predict absences) and sensitivity (i.e. ability to correctly predict presences). TSS ranges

259

from -1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance

260

no better than random (Pennino et al., 2017). Both AUC and TSS are used in combination when

261

evaluating the predictive power of a SDM (Pearson et al., 2006; Brodie et al., 2015). Model 11

262

validation was performed using the ‘PresenceAbsence’ package (Freeman and Moisen, 2008) of the

263

R software.

264 265

2.6. Comparison of model predictions

266 267

The monthly spatial structure similarity of predictions of species' distributions derived from

268

the two different data sources were compared using the Schoener's D (Schoener, 1968) and Warren's

269

I statistics (Warren et al., 2008). Both statistics assume probability distributions defined over

270

geographic space and are calculated using the difference between species suitability score at each

271

grid cell (Pennino et al., 2018). The two metrics ranged from 0, where two distributions have no

272

overlap, to 1, where two distributions are identical. These analyses were carried out using the

273

“nicheOverlap” function of the “dismo” package (Hijmans et al., 2016) of the R software. The

274

Kruskal-Wallis H test was used to elucidate whether there were significant differences between

275

months for each statistic index.

276

Moreover, Spearman rank correlation between monthly predictions derived from the two

277

different data sources was also calculated, as Schoener's D and Warren's I statistics may tend to

278

overestimate similarity between rasters when many grid cells are of similar values (e.g. when a large

279

part of the region contains low probability values). The Spearman rank correlation coefficient ranged

280

from -1 (perfectly negatively correlated) and 1 (perfectly positively correlated).

281 282 283

3. Results

284 285

3.1. Catch-independent data

286 12

287

The final B-HSM, after testing all the possible variable combinations, included SST, SSH,

288

EKE and days at sea as relevant predictors, included as linear terms. The spatial effect (W) and buoy

289

identification were included as random effects (Table 2).

290

According to the model, higher probability of tuna presence is expected in warmer waters (posterior

291

mean = 0.12; 95% CI = [0.09, 0.16]), with higher SSH values (posterior mean = 0.17; 95% CI =

292

[0.14, 0.20]) and low eddy kinetic energy values (posterior mean = -0.17; 95% CI = [-0.21, -0.13]).

293

Moreover, days at sea also explained the probability of presence of tuna, finding greater probabilities

294

with higher days at sea values (posterior mean = 0.03; 95% CI = [0.02, 0.03]) (Table 3).

295

Maps of the predicted probability of occurrence of tuna using catches-independent data show

296

relatively small and disaggregated hotspots. During autumn and part of winter (i.e. October to

297

January), hotspots are mostly distributed in the equatorial zone. During the rest of the year, small

298

hotspots distributed from -10°S to 10°N are shown, with some variation in February, where a marked

299

hotspot in the Horn of Africa is found, and in March, where highest probabilities were found close to

300

the coast of Somalia and Kenya (Fig. 3). The standard deviation maps (Fig. S.1) show very low

301

values in the area where data were collected, while the error increases along the edges and off the

302

data collection area.

303

The analyses of the model prediction performances obtained in the cross-validation showed AUC

304

values of 0.71±0.05, indicating a reasonable degree of discrimination between locations in which

305

tuna were present and those in which were absent. Likewise, TSS values of 0.31±0.05 were found,

306

which represents a reasonable ability of the model to predict true negative and true positive

307

predictions (Brodie et al., 2015).

308 309

3.2. Catch-dependent data

13

310 311

The main predictors of tuna habitats in the Western Indian Ocean were SAL, SST, SSH and a

312

stochastic spatial component that accounted for the residual spatial autocorrelation (W) (Table 4).

313

Specifically, results showed a positive relationship between SST (posterior mean = 0.17; 95% CI =

314

[0.07, 0.26]), SAL (posterior mean = 0.26; 95% CI = [0.18, 0.34]), and the probability of tuna

315

occurrence. On the contrary, the SSH showed a negative relationship (posterior mean = -0.20; 95%

316

CI = [-0.30, -0.09]) with tuna occurrence (Table 5).

317

From January to May the tuna probability hotspot moves from the central Western Indian Ocean to

318

the south, creating a second hotspot in the Mozambique Channel. In June, the hotspot in the

319

Mozambique Channel disappears and it becomes smaller and limited to the equatorial zone. From

320

July to December the hotspot remains more or less stable in the same area, with small variations.

321

Another less marked hotspot in the east can be also observed, varying in intensity throughout the year

322

and even disappearing in certain months (i.e. April, May and November) (Fig. 4). In this case, the

323

standard deviation maps (Fig. S.2) show low values along the edges and practically zero values in the

324

rest of the area.

325

The generation of pseudo-absences is an open and well-known discussion in the scientific

326

community. To check the reliability of our results, we generated randomly for the study area three

327

other sets of pseudo-absences and we ran again the models. The results were similar in terms of final

328

covariates selected, and model validation and prediction capabilities (i.e. performance scores) (Table

329

S.1, Fig. S.3), underlining that the model has been correctly established and that it can be a

330

representative selection of the average results.

331 332

AUC and TSS values showed reasonable values, with averages of 0.74±0.08) and 0.45±0.05, respectively, indicating good performance of the models for prediction as well. 14

333 334

3.3. Model comparison statistics

335 336

On the one hand, the similarity statistics indicate a good level of overlap in the predictions of

337

the two model outputs. In fact, for all months, Schoener's D ranged from about 0.43 to 0.67 (mean

338

0.58±0.07), while Warren's I ranged from 0.74 to 0.91 (mean 0.85±0.06). On the other hand, the

339

Spearman rank correlation coefficient indicated positive but not very high scores, which ranged from

340

0.43 to 0.67 (mean 0.58±0.07) (Table 6). Although intra-annual variations were detected, the

341

differences found between months are not statistically significant for either index (Kruskal Wallis,

342

p>0.05).

343 344 345

4.

Discussion

346 347

This study represents the first investigation of spatio-temporal distribution of tropical tuna

348

aggregated under the DFADs using observations from catches-dependent and independent data

349

sources. Moreover, this study provides new information about spatio-temporal habitat preferences of

350

tuna aggregated under the DFADs in the Western Indian Ocean.

351 352

4.1. Model comparison

353 354

Although the similarity statistics computed in this study (i.e. Schoener's D and Warren's I)

355

indicated a good level of overlap in the predictions of the two model outputs, Spearman rank

356

correlation showed low values. These differences could most likely be due to the areas of low

357

probability. This study considers large-scale habitat preference predictions in which, outside the main 15

358

fishing area, very low probabilities are estimated for both models using the two data sources, which

359

may result in misleading similarity statistics. This highlights the need of working with various

360

indices of similarity, and not only Schoener's D and Warren's I or Spearman rank correlation

361

statistics individually, to better understand the spatial structure of the results and obtain non-

362

misleading results.

363

Despite the two model outputs showing, in general, similar areas of tuna presence under the

364

DFADs, some discrepancies are present. While temperature is relevant, and it has a positive effect on

365

both models, SSH has the opposite effect in each model (i.e. positive in catch-independent model and

366

negative in catch-dependent model). Tuna’s preference for warm surface waters is consistent with

367

previous works, where warmer waters have been linked to both higher catch rates of large tuna

368

(Zagaglia et al., 2004; Lan et al., 2011; Lan et al., 2017) or acoustic biomass values (Lopez et al.,

369

2017b). The opposite effect found for SSH is more complex to interpret. Tropical tuna usually live in

370

warmer environments that display near null or positive SSH values (Arrizabalaga et al., 2015; Lan et

371

al., 2017; Lopez et al., 2017b), which is consistent with catch-independent model results (positive

372

SSH). Catch-dependent data also revealed a positive relationship with salinity, which is not identified

373

by catch-independent data as one of the main covariates. Also, the buoy-based model found a

374

negative relationship with EKE whereas the catch-based model did not find any significant

375

relationship. These differences in the relevance of environmental variables between the models could

376

be due to several factors. On the one hand, it must be taken into account that the environmental

377

variables were considered to play a less important role in DFADs than in for example free school

378

communities in the Indian Ocean (Lezama-Ochoa et al., 2015), which suggests that other ecological

379

processes may also be involved at different scales on the aggregation process of the species at

380

DFADs. On the other hand, in this study we have not considered the specific composition of the

381

species aggregated under the object, and it could be possible that the behavior of the species in

382

relation to the environmental variables varies between species. Current echo-sounder buoys provide a 16

383

single biomass value without determining species or size composition of the fish underneath the

384

DFAD, although work is being conducted to enable echo-sounder buoys provide abundance of tuna

385

by species in the near future (Boyra et al., 2018). Therefore, once this becomes possible, the specific

386

composition of the species aggregated under the DFADs should be included in the models, to

387

consider possible differences in species behavior and environmental variables.

388

It is also possible that the differences found between models are related to the characteristics

389

of the data (e.g. strata differences), data processing (e.g. aggregated or fine-scale environmental data)

390

or areas where catches occur and the trajectories of the DFADs. For example, there is a significant

391

difference in sampling size between the two datasets: for the same period, catch-independent data

392

have 42,322 observations compared to the 7357 of the catch-dependent data. These may be the main

393

reasons for the difference in the precision of the hotspots identified in the prediction maps. The maps

394

obtained with acoustic data allow identifying areas of high probability of tuna presence under the

395

DFADs with greater precision, whereas the maps derived from catch data do not allow observing any

396

variation on a finer scale. Maps illustrating the distribution of tuna associated with DFADs are

397

fundamental sources to help design management and conservation measures, so it is very useful to

398

obtain maps with as much detailed as possible within the large area occupied by tropical purse-seine

399

fisheries. With the availability of high-resolution data of tuna presence derived from fishers’ echo-

400

sounder buoys and environmental variables from satellite data, broad-scale studies can be made to

401

estimate tuna distribution at fine resolution.

402

One of the most important limitations of catch data is that they are constrained to the areas

403

where fishing occurs. Therefore, there are areas where, for several reasons (e.g. exclusive economic

404

zones without fishing agreements, distance to port, lower catch rates, etc.), fishing is not occurring

405

and, hence, there are no observations to build SDMs. In contrast, echo-sounders buoys are

406

continuously recording and providing information on the DFADs trajectories and biomass of fish

407

aggregated underneath in a non-invasive manner, covering areas and seasons where catch data are not 17

408

available. Based on our results, the application of spatio-temporal models for tuna associated with

409

DFADs using acoustic data provided by fishers´ echo-sounder buoys appears promising for the

410

identification of the spatial distribution of tuna in specific areas where fishing effort is limited or does

411

not exist. Nevertheless, echo-sounder buoys are sometimes switched off when they move away from

412

the fishing zone, to avoid communication fees for DFADs outside the practical range. However, in

413

this study, the presence of DFADs outside regular fishing grounds throughout the year suggest that

414

that is not the case for all the buoys. Despite all the advantages of acoustic data derived from fishers´

415

echo-sounder buoys, there are also some limitations. For example, using data obtained in DFADs,

416

only the distribution of small fish of the population can be modeled as large fish are mainly

417

unassociated to DFADs and, presently can only be modeled using a mixture of 3 species of tunas

418

(skipjack, yellowfin and bigeye tuna) and not a single species as it cannot still differentiate using

419

buoy’s echo-sounder biomass information.

420 421

4.2. Spatio-temporal distribution of tuna

422 423

Results showed month-specific distribution patterns, suggesting that tuna species may have

424

different habitat preferences depending on the season. The Indian Ocean is characterized by strong

425

environmental fluctuations associated with monsoon regimes that affect ocean circulation and

426

biological production (Schott and McCreary, 2001), seasonal variability in fishing grounds (Davies et

427

al., 2014), and catches (Kaplan et al., 2014). Changes in biophysical factors associated with

428

seasonality (i.e. chlorophyll, temperature, salinity, dissolved oxygen) may play an important role in

429

the distribution of tuna aggregated under floating objects.

430

In the case of the maps obtained using catch data, we observe a more general hotspot that

431

varies throughout the year but does not allow the identification of more precise areas. On the other

432

hand, the maps obtained with echo-sounder buoy data identify hotspots at a finer scale and the intra18

433

annual changes in tuna distribution. For example, an intermediate zone of low probability of tuna

434

presence was found between the two hotspots identified by the acoustic data (one in horn of Africa

435

and the another one around Seychelles and west of Seychelles) during February. This intermediate

436

zone of low probability of tuna presence may be related to the relatively high speed of the currents in

437

the area. In fact, tuna could be forced to leave floating objects in this area because of the high

438

energetic cost being associated with fast-DFADs carried by the South Equatorial Countercurrent

439

flowing eastwards at speeds > 1- 1.5 knots (Schott and McCreary, 2001; Tomczak and Godfrey,

440

2013). This finding is consistent with the study of Lopez et al. (2017b), that observed greater

441

abundances of tuna aggregated to DFADs at low speeds (i.e. 0.5-1 knots or lower). From July to

442

September, tuna distribution patterns seem to be related to the upwelling system of the Somali coast.

443

Indeed, during the summer monsoon (i.e. June to September) there is an extensive coastal upwelling

444

over the margin of continental areas in the northwestern Indian Ocean (Schott et al., 2009). Cold and

445

high-salinity waters shallow as a result of the upwelling with high nutrient contents that stimulate the

446

increase of primary production in the area (Veldhuis et al., 1997; Hitchcock et al., 2000). In addition,

447

one of the main characteristic of this upwelling is the rapid export of surface waters offshore reaching

448

over 500 km from the coast (Hitchcock et al., 2000; Wiggert et al., 2006), due to the presence of a

449

large anticyclonic eddy, called Great Whirl, in the vicinity of the upwelling zone (Beal and Donohue,

450

2013). These patterns are mainly reflected in the July and August prediction maps, which show the

451

highest probability of tuna in Somali waters. Moreover, in September higher tuna probability areas

452

move to the east. This may be related to the offshore export of surface productive waters after the

453

summer upwelling, indicating that tuna may aggregate to DFADs that follow this nutritive current

454

heading eastward. Although in this study chlorophyll has not been identified as a relevant variable,

455

previous studies have found positive relationship between tuna associated with DFADs and

456

chlorophyll concentration, particularly in the Indian Ocean (Fraile et al., 2010), but also in the

457

Atlantic Ocean (Lopez et al., 2017b). 19

458 459 460 461

5. Conclusions

462

The present study compares, for the first time, species distribution models developed using

463

catch-dependent (i.e. catch) data with models using catch-independent data (i.e. acoustic data derived

464

from fishers´ echo-sounder buoys), concluding that the maps obtained with acoustic data allow

465

identifying areas of high probability of tuna presence under the DFADs with greater precision than

466

maps derived from catch data. This highlights the great potential fishers´ echo-sounder buoys have

467

for addressing key scientific questions on tuna ecology and behavior. Furthermore, these data could

468

significantly contribute to the development of possible management measures, for example time–area

469

closures, in which, by having greater precision in the spatio-temporal distribution of tunas associated

470

with DFADs, management measures could be more efficient. The integration of novel sources of

471

information in science and management, such as fishers´ echo-sounder buoys data, are a valuable

472

tool.

473 474 475

Acknowledgments

476 477

The authors would like to thank the Spanish fishing company of tuna purse seiners in the

478

Indian Ocean ECHEBASTAR who kindly provided the acoustic data from their echo-sounder buoys

479

used in the present study.

480 481 482

References 20

483 484 485 486

Allouche, O., Tsoar, A., Kadmon, R., 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J Appl Ecol. 43, 1223-1232. Arrizabalaga, H., Dufour, F., Kell, L., Merino, G., Ibaibarriaga, L., Chust, G., Irigoien, X., Santiago,

487

J., Murua, H., Fraile, I., 2015. Global habitat preferences of commercially valuable tuna. Deep

488

Sea Res. Part 2 Top. Stud. Oceanogr. 113, 102-112.

489

Báez, J.C., Fernández, F., Pascual-Alayón, P.J., Ramos, M.L., Deniz, S., Abascal, F., 2018. Updating

490

the statistics of the EU-Spain purse seine fleet in the Indian Ocean (1990-2017). IOTC-2018-

491

WPTT20-15.

492 493 494

Beal, L., Donohue, K., 2013. The Great Whirl: Observations of its seasonal development and interannual variability. J Geophys Res-Oceans 118, 1-13. Boyra, G., Moreno, G., Sobradillo, B., Pérez-Arjona, I., Sancristobal, I., Demer, D.A., 2018. Target

495

strength of skipjack tuna (Katsuwanus pelamis) associated with fish aggregating devices (FADs).

496

ICES J Mar Sci. 75, 1790-1802.

497

Bradley, D., Merrifield, M., Miller, K.M., Lomonico, S., Wilson, J.R., Gleason, M.G., 2019.

498

Opportunities to improve fisheries management through innovative technology and advanced

499

data systems. Fish Fish. 20 564-583.

500

Brodie, S., Hobday, A.J., Smith, J.A., Everett, J.D., Taylor, M.D., Gray, C.A., Suthers, I.M., 2015.

501

Modelling the oceanic habitats of two pelagic species using recreational fisheries data. Fish

502

oceanogr. 24, 463-477.

503

Castro, J., Santiago, J., Santana-Ortega, A., 2002. A general theory on fish aggregation to floating

504

objects: an alternative to the meeting point hypothesis. Rev. Fish Biol. Fish. 11, 255-277.

505

Coelho, R., Mejuto, J., Domingo, A., Yokawa, K., Liu, K.M., Cortés, E., Romanov, E.V., Silva, C.,

506

Hazin, F., Arocha, F., 2017. Distribution patterns and population structure of the blue shark

507

(Prionace glauca) in the Atlantic and Indian Oceans. Fish Fish. 19 90-106.

508

Conn, P.B., Thorson, J.T., Johnson, D.S., 2017. Confronting preferential sampling when analysing

509

population distributions: diagnosis and model‐based triage. Methods Ecol Evol. 8, 1535-1546.

510 511 512 513

Chen, I., LEE, P.F., TZENG, W.N., 2005. Distribution of albacore (Thunnus alalunga) in the Indian Ocean and its relation to environmental factors. Fish oceanogr. 14, 71-80. Dagorn, L., Bez, N., Fauvel, T., Walker, E., 2013. How much do fish aggregating devices (FADs) modify the floating object environment in the ocean? Fish oceanogr. 22, 147-153.

21

514 515 516

Davies, T.K., Mees, C.C., Milner-Gulland, E., 2014. The past, present and future use of drifting fish aggregating devices (FADs) in the Indian Ocean. Mar Policy. 45, 163-170. Dell'Apa, A., Pennino, M.G., Bonzek, C., 2017. Modeling the habitat distribution of spiny dogfish

517

(Squalus acanthias), by sex, in coastal waters of the northeastern United States. Fish. Bull. 115,

518

89-100.

519 520

Fielding, A.H., Bell, J.F., 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ conserv. 24, 38-49.

521

Filmalter, J.D., Dagorn, L., Cowley, P.D., Taquet, M., 2011. First descriptions of the behavior of

522

silky sharks, Carcharhinus falciformis, around drifting fish aggregating devices in the Indian

523

Ocean. Bull. Mar. Sci. 87, 325-337.

524 525 526

Fonseca, V.P., Pennino, M.G., de Nóbrega, M.F., Oliveira, J.E.L., de Figueiredo Mendes, L., 2017. Identifying fish diversity hot-spots in data-poor situations. Mar. Environ. Res. 129, 365-373. Fonteneau, A., Chassot, E., Bodin, N., 2013. Global spatio-temporal patterns in tropical tuna purse

527

seine fisheries on drifting fish aggregating devices (DFADs): Taking a historical perspective to

528

inform current challenges. Aquat. Living Resour. 26, 37-48.

529 530

Fonteneau, A., Lucas, V., Tewkai, E., Delgado, A., Demarcq, H., 2008. Mesoscale exploitation of a major tuna concentration in the Indian Ocean. Aquat. Living Resour. 21, 109-121.

531

Forget, F.G., Capello, M., Filmalter, J.D., Govinden, R., Soria, M., Cowley, P.D., Dagorn, L., 2015.

532

Behaviour and vulnerability of target and non-target species at drifting fish aggregating devices

533

(FADs) in the tropical tuna purse seine fishery determined by acoustic telemetry. Can. J. Fish.

534

Aquat. Sci. 72, 1398-1405.

535

Fraile, I., Murua, H., Goni, N., Caballero, A., 2010. Effects of environmental factors on catch rates of

536

FAD-associated yellowfin (Thunnus albacares) and skipjack (Katsuwonus pelamis) tunas in the

537

western Indian Ocean, IOTC-2010-WPTT-46.

538 539 540 541 542 543 544 545

Freeman, E.A., Moisen, G., 2008. PresenceAbsence: An R package for presence absence analysis. Journal of Statistical Software. 23 (11): 31 p. Govinden, R., Dagorn, L., Filmalter, J., Soria, M., 2010. Behaviour of Tuna a ssociated with Drifting Fish Aggregating Devices (FADs) in the Mozambique Channel, IOTC-2010-WPTT-25. Hastie, T., Fithian, W., 2013. Inference from presence‐only data; the ongoing controversy. Ecography 36, 864-867. Held, L., Schrödle, B., Rue, H., 2010. Posterior and cross-validatory predictive checks: a comparison of MCMC and INLA. Statistical modelling and regression structures, 91-110.

22

546 547 548 549 550 551 552 553 554

Hijmans, R., Phillips, S., Leathwick, J., Elith, J., 2016. dismo: species distribution modeling. R package ver. 1.0-15. Hijmans, R., van Etten, J., Cheng, J., Mattiuzzi, M., Sumner, M., Greenberg, J., 2017. Raster: Geographic Data Analysis and Modeling. R package version 2.3–33; 2016. Hitchcock, G.L., Key, E.L., Masters, J., 2000. The fate of upwelled waters in the Great Whirl, August 1995. Deep Sea Res. Part 2 Top. Stud. Oceanogr. 47, 1605-1621. Jury, M., McClanahan, T., Maina, J., 2010. West Indian ocean variability and east African fish catch. Mar. Environ. Res. 70, 162-170. Kaplan, D.M., Chassot, E., Amandé, J.M., Dueri, S., Demarcq, H., Dagorn, L., Fonteneau, A., 2014.

555

Spatial management of Indian Ocean tropical tuna fisheries: potential and perspectives. ICES J

556

Mar Sci. 71, 1728-1749.

557

Kinas, P.G., Andrade, H.A., 2017. Introdução à análise bayesiana (com R). Consultor Editorial.

558

Kuhn, M., Johnson, K., 2013. Applied predictive modeling. Springer.

559

Lan, K.-W., Lee, M.-A., Lu, H.-J., Shieh, W.-J., Lin, W.-K., Kao, S.-C., 2011. Ocean variations

560

associated with fishing conditions for yellowfin tuna (Thunnus albacares) in the equatorial

561

Atlantic Ocean. ICES J Mar Sci. 68, 1063-1071.

562

Lan, K.-W., Shimada, T., Lee, M.-A., Su, N.-J., Chang, Y., 2017. Using Remote-Sensing

563

Environmental and Fishery Data to Map Potential Yellowfin Tuna Habitats in the Tropical

564

Pacific Ocean. Rem Sen. 9, 444.

565 566

Lee, P.-F., Chen, I., Tzeng, W.-N., 2005. Spatial and temporal distribution patterns of bigeye tuna (Thunnus obesus) in the Indian Ocean. Zool stud. 44, 260-270.

567

Lennert-Cody, C.E., Moreno, G., Restrepo, V., Román, M.H., Maunder, M.N., Poos, H.e.J.J., 2018.

568

Recent purse-seine FAD fishing strategies in the eastern Pacific Ocean: what is the appropriate

569

number of FADs at sea? ICES J Mar Sci. 75, 1748-1757.

570

Lezama-Ochoa, N., Murua, H., Chust, G., Ruiz, J., Chavance, P., de Molina, A.D., Caballero, A.,

571

Sancristobal, I., 2015. Biodiversity in the by-catch communities of the pelagic ecosystem in the

572

Western Indian Ocean. Biodivers conserv. 24, 2647-2671.

573

Lezama-Ochoa, N., Murua, H., Chust, G., Van Loon, E., Ruiz, J., Hall, M., Chavance, P., Delgado

574

De Molina, A., Villarino, E., 2016. Present and future potential habitat distribution of

575

Carcharhinus falciformis and Canthidermis maculata by-catch species in the tropical tuna purse-

576

seine fishery under climate change. Front. Mar. Sci. 3, 34.

23

577

Lindgren, F., Rue, H., Lindström, J., 2011. An explicit link between Gaussian fields and Gaussian

578

Markov random fields: the stochastic partial differential equation approach. J. R. Stat. Soc. 73,

579

423-498.

580

Lopez, J., Moreno, G., Boyra, G., Dagorn, L., 2016. A model based on data from echosounder buoys

581

to estimate biomass of fish species associated with fish aggregating devices. Fish. Bull. 114,

582

166-178.

583

Lopez, J., Moreno, G., Ibaibarriaga, L., Dagorn, L., 2017a. Diel behaviour of tuna and non-tuna

584

species at drifting fish aggregating devices (DFADs) in the Western Indian Ocean, determined

585

by fishers’ echo-sounder buoys. Mar. Biol. 164, 44.

586

Lopez, J., Moreno, G., Lennert-Cody, C., Maunder, M., Sancristobal, I., Caballero, A., Dagorn, L.,

587

2017b. Environmental preferences of tuna and non-tuna species associated with drifting fish

588

aggregating devices (DFADs) in the Atlantic Ocean, ascertained through fishers’ echo-sounder

589

buoys. Deep Sea Res. Part 2 Top. Stud. Oceanogr. 140, 127-138.

590

Lopez, J., Moreno, G., Sancristobal, I., Murua, J., 2014. Evolution and current state of the technology

591

of echo-sounder buoys used by Spanish tropical tuna purse seiners in the Atlantic, Indian and

592

Pacific Oceans. Fish. Res. 155, 127-137.

593 594 595 596 597 598 599

Lopez, J., Moreno, G., Soria, M., Cotel, P., Dagorn, L., 2010. Remote discrimination of By-cath in purse seine fishery using fishers' echo-sounder buoys. IOTC-2010-WPEB-03. Lowell, B., Mustain, P., Ortenzi, K., Warner, K., 2015. One name, one fish: why seafood names matter. Oceana, 1-12. Martínez-Minaya, J., Cameletti, M., Conesa, D., Pennino, M.G., 2018. Species distribution modeling: a statistical review with focus in spatio-temporal issues. Stoch. Env. Res. Risk A, 1-18. Moreno, G., Dagorn, L., Capello, M., Lopez, J., Filmalter, J., Forget, F., Sancristobal, I., Holland, K.,

600

2016. Fish aggregating devices (FADs) as scientific platforms. Fish. Res. 178, 122-129.

601

Moreno, G., Josse, E., Brehmer, P., Nøttestad, L., 2007. Echotrace classification and spatial

602

distribution of pelagic fish aggregations around drifting fish aggregating devices (DFAD).

603

Aquat. Living Resour. 20, 343-356.

604

Muñoz, F., Pennino, M.G., Conesa, D., López-Quílez, A., Bellido, J.M., 2013. Estimation and

605

prediction of the spatial occurrence of fish species using Bayesian latent Gaussian models. Stoch.

606

Env. Res. Risk A. 27, 1171-1180.

607

Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R., O’hara, R., Simpson, G.L.,

608

Solymos, P., Stevens, M.H.H., Wagner, H., 2013. Package ‘vegan’. Community ecology

609

package, version 2. 24

610

Orue, B., Lopez, J., Moreno, G., Santiago, J., Boyra, G., Uranga, J., Murua, H., 2019a. From fisheries

611

to scientific data: A protocol to process information from fishers’ echo-sounder buoys. Fish. Res.

612

215, 38-43.

613

Orue, B., Lopez, J., Moreno, G., Santiago, J., Soto, M., Murua, H., 2019b. Aggregation process of

614

drifting fish aggregating devices (DFADs) in the Western Indian Ocean: Who arrives first, tuna

615

or non-tuna species? PloS one 14, e0210435.

616

Pearson, R.G., Thuiller, W., Araújo, M.B., Martinez‐Meyer, E., Brotons, L., McClean, C., Miles, L.,

617

Segurado, P., Dawson, T.P., Lees, D.C., 2006. Model‐based uncertainty in species range

618

prediction. J. Biogeogr. 33, 1704-1711.

619

Pennino, M.G., Conesa, D., López-Quílez, A., Munoz, F., Fernández, A., Bellido, J.M., 2016.

620

Fishery-dependent and-independent data lead to consistent estimations of essential habitats.

621

ICES J Mar Sci. 73, 2302-2310.

622

Pennino, M.G., Mérigot, B., Fonseca, V.P., Monni, V., Rotta, A., 2017. Habitat modeling for

623

cetacean management: Spatial distribution in the southern Pelagos Sanctuary (Mediterranean

624

Sea). Deep Sea Res. Part 2 Top. Stud. Oceanogr. 141, 203-211.

625 626

Pennino, M.G., Muñoz, F., Conesa, D., López-Quílez, A., Bellido, J.M., 2014. Bayesian spatiotemporal discard model in a demersal trawl fishery. J. Sea Res. 90, 44-53.

627

Pennino, M.G., Rufener, M.-C., Thomé-Souza, M.J., Carvalho, A.R., Lopes, P.F., Sumaila, U.R.,

628

2018. Searching for a compromise between biological and economic demands to protect

629

vulnerable habitats. Sci. Rep. 8, 7791.

630 631 632

R Development Core Team, 2017. R: A language and environment for statistical computing. R.Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Rajapaksha, J.K., Samarakoon, L., Gunathilaka, A.A.J.K., 2013. Environmental Preferences of

633

Yellowfin Tuna in the North East Indian Ocean: An Application of Satellite Data to Longline

634

Catches. Int. J. Fish. Aquat.Stud 2(4): 72-80.

635

Robert, M., Dagorn, L., Lopez, J., Moreno, G., Deneubourg, J.-L., 2013. Does social behavior

636

influence the dynamics of aggregations formed by tropical tunas around floating objects? An

637

experimental approach. J. Exp. Mar. Biol. 440, 238-243.

638

Robinson, N.M., Nelson, W.A., Costello, M.J., Sutherland, J.E., Lundquist, C.J., 2017. A Systematic

639

Review of Marine-Based Species Distribution Models (SDMs) with Recommendations for Best

640

Practice. Front. Mar. Sci. 4, 421.

641 642

Roos, M., Held, L., 2011. Sensitivity analysis in Bayesian generalized linear mixed models for binary data. Bayesian Anal. 6, 259-278. 25

643 644 645 646 647 648 649 650 651 652 653 654 655

Roos, N.C., Carvalho, A.R., Lopes, P.F., Pennino, M.G., 2015. Modeling sensitive parrotfish (Labridae: Scarini) habitats along the Brazilian coast. Mar. Environ. Res. 110, 92-100. Rue, H., Martino, S., Chopin, N., 2009. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. J. R. Stat. Soc. B 71, 319-392. Schoener, T.W., 1968. The Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49, 704-726. Schott, F.A., McCreary, J.P., 2001. The monsoon circulation of the Indian Ocean. Prog. Oceanogr. 51, 1-123. Schott, F.A., Xie, S.P., McCreary, J.P., 2009. Indian Ocean circulation and climate variability. Rev. Geophys. 47. Sequeira, A., Mellin, C., Rowat, D., Meekan, M.G., Bradshaw, C.J., 2012. Ocean‐scale prediction of whale shark distribution. Divers. Distrib. 18, 504-518. Sequeira, A.M., Mellin, C., Delean, S., Meekan, M.G., Bradshaw, C.J., 2013. Spatial and temporal

656

predictions of inter-decadal trends in Indian Ocean whale sharks. Mar. Ecol. Prog. Ser. 478, 185-

657

195.

658 659 660

Simmonds, J., MacLennan, D., 2005. Fishery acoustic theory and practice, Blackwell Scientific Publications, Oxford, UK. Song, L.M., Zhang, Y., Xu, L.X., Jiang, W.X., Wang, J.Q., 2008. Environmental preferences of

661

longlining for yellowfin tuna (Thunnus albacares) in the tropical high seas of the Indian Ocean.

662

Fish oceanogr. 17, 239-253.

663 664

Tille, Y., Matei, A., 2016. sampling: Survey Sampling. R package version 2.8. https://CRAN.Rproject.org/package=sampling.

665

Tomczak, M., Godfrey, J.S., 2013. Regional oceanography: an introduction. Elsevier.

666

Veldhuis, M.J., Kraay, G.W., Van Bleijswijk, J.D., Baars, M.A., 1997. Seasonal and spatial

667

variability in phytoplankton biomass, productivity and growth in the northwestern Indian Ocean:

668

the southwest and northeast monsoon, 1992–1993. Deep Sea Res. Part I: Oceanogr. Res. Pap. 44,

669

425-449.

670 671

Warren, D.L., Glor, R.E., Turelli, M., 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868-2883.

672

Watanabe, S., 2010. Asymptotic equivalence of Bayes cross validation and widely applicable

673

information criterion in singular learning theory. J. Mach. Learn. Res. 11, 3571-3594.

26

674

Wiggert, J., Murtugudde, R., Christian, J., 2006. Annual ecosystem variability in the tropical Indian

675

Ocean: Results of a coupled bio-physical ocean general circulation model. Deep Sea Res. Part 2

676

Top. Stud. Oceanogr. 53, 644-676.

677

Zagaglia, C.R., Lorenzzetti, J.A., Stech, J.L., 2004. Remote sensing data and longline catches of

678

yellowfin tuna (Thunnus albacares) in the equatorial Atlantic. Remote sens. environ. 93, 267-

679

281.

680 681 682 683

Table 1. Predictor variables used for modeling tuna occurrence using catches dependent and independent data in the Indian Ocean. Variables Variable name

Spatial

Temporal

resolution

resolution

Unit

acronym SST

Sea Surface Temperature

ºC

0.25°

Daily

SAL

Salinity

PSU

0.25°

Daily

SSH

Sea Surface Height

m

0.25°

Daily

VEL

Velocity

m s-1

0.25°

Daily

EKE

Eddy Kinetic Energy

m2s2

0.25°

Daily

HEADING

Direction

degrees

0.25°

Daily

CHL

Chlorophyll concentration

mg m-3

0.50°

Weekly

O2

Oxygen Concentration

mmol m-3

0.50°

Weekly

684 685 686 687 688 689

Table 2. Comparison of the most relevant models for the tuna B-HSMs selected using catchindependent data. Statistics acronyms are: Watanabe Akaike Information Criterion (WAIC) and Logarithmic Cross Validated Score (LCPO). Predictor acronyms are: W= spatial effect, SST = Sea Surface Temperature, SSH = Sea Surface Height, SAL = Salinity, O2 = Oxygen Concentration, EKE= Kinetic Energy, CHL= Chlorophyll. The best model is highlighted in bold. MODEL

WAIC

LCPO

b0+W+SST+SSH+EKE+Days at sea+ buoy ID

23415.36

0.53

b0+W+SST+SSH+EKE+ buoy ID

25489.22

0.55

27

b0+W+EKE+ buoy ID

25989.14

0.55

b0+W+SST+SSH+EKE+Days at sea

29475.19

0.57

b0+ CHL+SAL+SST+SSH+O2+EKE

30451.66

0.57

690 691 692 693 694 695

Table 3. Numerical summary of the marginal posterior distribution of the fixed effects for the best tuna B-HSMs selected using catch-independent data. For each variable, the mean, standard deviation, and a 95% credible central interval (Q0.025- Q0.975) is provided, containing 95% of the probability under the posterior distribution. Predictor acronyms are: W= spatial effect, B = buoy identification, SST = Sea Surface Temperature, SSH = Sea Surface Height, EKE = Eddy Kinetic Energy.

696 VARIABLE MEAN

MODEL

SD

Q0.025

Q0.975

SSH

0.17

0.02

0.14

0.20

SST

0.12

0.01

0.09

0.16

EKE

-0.17

0.02

-0.21

-0.13

Days at sea

0.03

0.01

0.02

0.03

b0 +W +SSH +SST +EKE +Days at sea+ B

697 698 699 700 701 702

Table 4. Comparison of the most relevant models for the tuna B-HSMs selected using catchdependent data. Statistics acronyms are: Watanabe Akaike Information Criterion (WAIC) and Logarithmic Cross Validated Score (LCPO). Predictor acronyms are: W= spatial effect, SST = Sea Surface Temperature, SSH = Sea Surface Height, SAL = Salinity, HEADING= Direction. The best model is highlighted in bold. MODEL

WAIC

LCPO

b0 +W +SST +SSH +SAL

8553.15

0.33

b0 +W +SST +SSH + HEADING+SAL

8599.92

0.33

b0 +W +SST +SSH

8600.34

0.34

b0+W

10214.16

0.42

b0

17878.28

0.69

703 28

704 705 706 707 708

Table 5. Numerical summary of the marginal posterior distribution of the fixed effects for the best tuna B-HSMs selected using catch-dependent data. For each variable, the mean, standard deviation, and a 95% credible central interval (Q0.025- Q0.975) is provided, containing 95% of the probability under the posterior distribution. Predictor acronyms are: W= spatial effect, SST = Sea Surface Temperature, SSH = Sea Surface Height, SAL = Salinity. MODEL

b0 +W +SST +SSH +SAL

VARIABLE

MEAN

SD

Q0.025 Q0.975

SST

0.17

0.05 0.07

SSH

-0.20

0.05 -0.30 -0.09

SAL

0.26

0.04 0.18

0.26

0.34

709 710

Table 6. Similarity statistics and spearman rank correlation coefficients obtained for each month Month Schoener's D Warren's I 0.647 0.907 January 0.624 0.894 February 0.523 0.792 March 0.512 0.731 April 0.586 0.871 May 0.655 0.917 June 0.432 0.743 July 0.586 0.874 August 0.864 September 0.558 0.662 0.902 October 0.874 November 0.601 0.671 0.911 December

Spearman rank correlation (p-value) 0.647 (<0.001) 0.624 (<0.001) 0.523 (<0.001) 0.512 (<0.001) 0.586 (<0.001) 0.645 (<0.001) 0.431 (<0.001) 0.586 (<0.001) 0.558 (<0.001) 0.662 (<0.001) 0.601 (<0.001) 0.671 (<0.001)

711

29

712 713

Fig. 1. Map of the study area (Western Indian Ocean) and spatial distribution of the tropical tuna

714

catch data of the purse-seine fisheries (green points) with a spatial resolution of 1° latitude/longitude

715

grid cell and a monthly temporal resolution.

30

716 717

Fig. 2. Characteristics of the Satlink echo-sounder buoy: Beam angle (a), depth range (h), and

718

diameter (d) at 115 m. An example of the echogram display for the 10 depth layers (ranging from 3

719

m to 115 m) (taken from Lopez et al. (2016)).

720

31

721 722 723

Fig. 3. Maps of the posterior mean of tuna occurrence probability for different months using catchindependent data.

724

32

725 726 727

Fig. 4. Maps of the posterior mean of tuna occurrence probability for different months using catchdependent data.

728 729 730

Supplementary material

731 732

Table S.1. Final covariates selected and model validation and prediction capabilities for each of the pseudo-absence creation tests. Model

Selected variables

WAIC

LCPO

Test 1

SST +SSH +SAL

8591.25 0.333

Test 2

SST +SSH +SAL

8549.14 0.331

33

Test 3

SST +SSH +SAL

8557.90 0.331

Model selected

SST +SSH +SAL

8553.15 0.330

733

734 735 736

Fig. S.1. Maps of the standard deviation of tuna occurrence probability for different months using catch-independent data.

34

737 738 739

Fig. S.2. Maps of the standard deviation of tuna occurrence probability for different months using catch-dependent data.

35

740 741 742

Fig. S.3. Maps of the posterior mean of tuna occurrence probability for different months for each of the pseudo-absence creation tests.

36

This manuscript describes an original work and is not under consideration by any other journal. All authors approved the manuscript and this submission. The authors have no conflicts of interest to declare.