Data collection

Otoliths and length data were collected throughout the Regional Tuna Tagging Project (RTTP), a large scale mark-recapture program, and the West Sumatra Tuna Tagging Project (WSTTP), a simple capture program, that were carried out by the IOTC. The tagging operations of RTTP were conducted between 2005–2007 on two pole-and-line vessels that were chartered to operate in the Western Indian Ocean. Field operations consisted in catching tunas, tagging them on a vinyl-covered cradle, measuring their fork length (fish length from the front to the fork in the center of the tail; ) through marks printed directly on the cradle and finally releasing them at sea [47]. The tagging operations of RTTP were conducted between 2005–2007 on two pole-and-line vessels that were chartered to operate in the Western Indian Ocean. Field operations consisted in catching tunas, tagging them on a vinyl-covered cradle, measuring their fork length (fish length from the front to the fork in the center of the tail; ) through marks printed directly on the cradle and finally releasing them at sea [47]. The fish were tagged with dart tags inserted into the musculature, below the second dorsal fin. The date and geographic location were recorded for each tag event. Some fish received a OTC injection, an antibiotic that is rapidly incorporated into calcified parts such as bones, scales, and otoliths and leaves a permanent fluorescent mark in the growth increment being formed at the time of tagging. Depending on fish size, between 1.5–3 mL of OTC were injected with a syringe in the intramuscular region of the back [47].

Recovery operations took place across the entire IO basin between 2005–2012. The majority of reported recoveries came from fish caught by the European purse seiners (IOTC 2011). of recovered fish was measured with calipers or a tape measure to the nearest 0.5 cm. The accuracy in the date and location of recaptures is dependent on the place and the process in which the tag is recovered. About 20% of the recoveries were made during purse seine fishing operations, which resulted in the recovered fish being associated with only one position and date. In contrast, tunas recovered during the unloading of purse seiners could be associated with several dates and catch locations due to the process of storing tunas in refrigerated wells that contain up to 5 sets collected over an entire fishing trip. The recovery can also occur downstream of the unloading process or in the canneries. The range of dates associated with each recapture was derived from logbook data and well maps and conducted in close collaboration with the IOTC and the purse seine fishing industry.

The WSTTP program was conducted in August 2007 on a pole-and-line vessel in the off western Indonesia. The fork length of fish measured either in a marked cradle or using a calipers, and the precisely date and geographic location were recorded for each fish caught.

Otolith sampling, preparation and reading

Sagittal otoliths were collected from 128 yellowfin recovered in the Western Indian Ocean through the RTTP, including 124 OTC-tagged fish, and measuring between 43 to 85 cm at tagging and 47.9 to 146.5 cm at recapture. Otoliths were also collected from 38 yellowfin captured during the WSTTP measuring between 19 to 46.6 cm (Table S1, Fig. S1). Otoliths were extracted, rinsed in water to remove tissue, and stored dry.

Yellowfin have fragile, thin and elliptical otoliths that require particular care during their preparation and subsequent interpretation of microstructural features [48]. All the otoliths collected were analyzed at the "Laboratoire de Sclérochronologie des Animaux Aquatiques" (LASAA) in Brest, France. Otoliths were prepared for age analysis using the following method [49]; [20]; [14]. They were cleaned in sodium hypochlorite and rinsed with distilled water before being embedded in resin blocks and transversally cut on both sides of the nucleus. The sections containing the nucleus were then fixed to a glass slide using thermoplastic glue and sanded to the level of the nucleus using different alumina grains (0.3 to 3 µm). The operation was performed on each side of the section until a slice of about 100 µm thickness was decalcified with EDTA (tri-sodium-ethylene-diamine-tetra-acetic acid) to increase the contrast between increments. The thin slides were examined under a microscope (1000×magnification) to count increments along the counting path on the sagitta, i.e. from the primordium, original point of growth, to the last increment deposited on the maximal growth axis.

Otoliths collected from OTC-tagged fishes, an increment count was made for different otolith sections: (i) between the nucleus and the OTC mark (), (ii) between the OTC mark and the edge () and (iii) between the nucleus and the edge () (Fig. 1). For fish that were not chemically tagged were read in full (). All otolith readings were performed by the same reader. Each otolith was read two to five times without prior knowledge on length or time-at-liberty of the individuals sampled so as to maintain certain independence between the multiple readings.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Otoliths of yellowfin tuna (external right and internal left; a) and the different sections used for reading the number of increments.

OTC: Oxytetracycline; : section from the nucleus to the OTC mark; : section from the OTC mark to the edge; : section from the nucleus to the edge; : Time-at-Liberty.

https://doi.org/10.1371/journal.pone.0060886.g001

Ageing error model

In this section, a hierarchical model was developed to estimate the age of each fish. The stochastic processes associated with otolith preparation and reading were modeled with the choice of an error structure and informative prior density functions based on expert judgment (Tables 1 and 2). In a first step, the hypothesis of daily increment deposition in otoliths, that has been observed in eastern Pacific yellowfin tuna [48]; [50], was tested using a subset of OTC-tagged otolith data. Information on the increment deposition process was subsequently used for estimating yellowfin age based on counts of otolith increments.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Parameters and variables used in the ageing error and somatic growth models.

https://doi.org/10.1371/journal.pone.0060886.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Deterministic and stochastic processes used in the ageing error and growth models. All variables are defined in table*#146;1.

https://doi.org/10.1371/journal.pone.0060886.t002

Evaluating model performance through simulations

Different simulations were performed to test and validate the ageing error model (Fig. S2). The first simulation aimed to assess the accuracy of the model in estimating . For this simulation, two alternative runs with and without individual variability in increment formation were considered, i.e. fixed to 0.95 and varying according to a normal distribution around 0.95. This simulation was repeated three times using datasets of 25 individuals whose time-at-liberty were generated according to a uniform distribution between 30 and 970 days. Four noisy readings of the number of increments were then generated according to Eq. D1 and S1 (2).

A second set of simulations was performed using 2 to 5 readings of the same otolith to evaluate the model's ability to accurately estimate fish age and its relevance relative to the traditional method which is based on averaging individual-specific increment counts. An intermediate method was also considered that based on averaging individual increment counts, but also account for potential bias in the increment deposition. For these simulations, was fixed to 0.95. The approach consisted of simulating realistic ages-at-recapture from which increment counts were derived with respect to , , and (2 Eq. D4). Noise was then added to the increment counts by randomly generating repeated readings (Eq. S3). Five age classes were considered, from 6 months to 5 years. Three simulations were performed using datasets composed of 500 individuals, i.e. 100 individuals for each age class whose ages were generated according to a uniform distribution between the minimum and the maximum values of the class. The accuracy of age estimates was assessed with the relative root mean square error (RMSE), a normalized indicator that measures the discrepancy between simulated and estimated ages. For each individual, the RMSE was calculated as follows:where is the estimated age and is the simulated age.

Coupling the ageing error and the growth models

In this section, we implemented a hierarchical model that allows for growth variations according to an individual's specific stochastic process. For IO yellowfin, modal progression analysis [30]; [32]; [33] and the preliminary analysis of the RTTP data [51] indicated a succession of phases of growth deceleration and acceleration, supporting the use of a two-stanza growth model. Accordingly, we used the Von Bertalanffy logK model (VB-logK) developed for southern bluefin tuna (Thunnus maccoyii), which allows for a smooth transition between two different growth rate coefficients ( and ) by modeling changes in growth using a logistic function [52]; [53]. The fork length of fish at the opportunity of capture , at tagging and at recapture, was then modeled as:where errors in length measurements, , were assumed to be independent and normally distributed around zero with a common variance .

The somatic VB-logK growth model was coupled to the ageing error model and fitted to the data using Bayesian inference (Table 2, Eq. D6 and D7). Consider the vector of growth parameters,  = {, , , , , , }, the vector of ageing parameters  = {, , }, and the relative percentage of misread otolith increments, . , , denote the posterior distributions of these parameters and , , denote their prior distributions. Here, and are not directly observable latent variables. The full model corresponds to the joint distribution of parameters and latent variables. As and are independent, this joint posterior distribution can be written as:(L1)

where represents the conditional Gaussian likelihood of the observed lengths. Thus, the length values were predicted from the marginal posterior distribution:(L2)

The growth rate coefficients and are partly model-specific, thus weakly informative priors were assigned to them. was assumed to vary according to a gamma prior distribution with a mean and coefficient of variation determined from the literature. was set equal to with following a uniform distribution (Table 2 Eq. P6 and P7). The transition rate between and , , which is specific to the VB-logK model, and the theoretical age of zero length, , which depends on the data, were assigned uniform distributions (Eq.P9 and P10). The parameter is the mean age relative to at which the change in growth occurs. This was assigned a weakly informative prior gamma distribution with a mean defined from the literature (Eq. P8; [36]; [37]; [32]; [33]). The standard deviation of size measurement errors, , was determined from the differences in the fork length of RTTP-IO fish released and recaptured several times with time-at-liberty less than or equal to seven days. These individuals were not included in subsequent analyses and therefore constitute an independent data set (Eq. P11).

The asymptotic length, , is a particularly important parameter because it determines the shape of the second part of the growth curve. As the dataset included limited information on the asymptotic part of the growth curve, auxiliary information was provided for this parameter consistently with the available knowledge on the biology of the species. An informative prior distribution was defined for using a generalized extreme value distribution (GEV). This allows for the extrapolation of the behavior of distribution tails from the greatest values of a sample and thus estimates the probability of the occurrence of extreme events [54]. The choice of this distribution is motivated by the fact that yellowfin grow throughout their life, such that the largest observed sizes should correspond to the oldest fish. The distribution was fit based on size measurement data from fresh fish that were collected during 1952–2011 from the European and Seychelles purse seine fisheries, Maldivian pole and line vessels, and Taiwanese and Japanese longliners. The observed maximum fork length from each measurement platform, i.e. either on vessels board or in cannery, fishery, and year was considered to represent independent random variables (,...,) with common continuous distribution function . Asymptotic length was then estimated from the approximation of the upper tail of by using the distribution:where is a location parameter, is the scale parameter (0) and is a tail index (shape parameter). These parameters were estimated using the maximum likelihood method [55].

The estimates of ages and growth parameters were evaluated based on 350,000 samples, thinned to one draw every 1000 sample, from Markov Chain Monte Carlo (MCMC) simulations of the joint posterior distribution. A burn-in period of 5,000 iterations was rejected. Three MCMCs were produced using a Gibbs sampler as implemented in OpenBugs version 3.2.1 [56]. The convergence of the MCMC to a stationary posterior distribution was visually evaluated from the Gelman-Rubin diagnostic, which was based on the ratio of inter-chain variance on intra-chain variance, i.e. it must be close to one in order to converge [57].

Testing the influence of the ageing technique on growth modeling

Simulations were performed to evaluate the contribution of the error model on growth parameter estimates and resulting age classification of commercial catches. To illustrate, a simple tagging dataset was simulated and the relationship specified from posterior marginal mode that was obtained with the growth model coupled with the ageing error model (section and Table 3) was used as reference growth curve. In this approach, four fork lengths for each 2 cm class, from 20 to 146 cm, were simulated to obtain a data set of 252 fishes. As the VB-logK curve is irreversible, corresponding ages were deduced by minimizing the difference between the simulated and expected length (E1). These ages were converted into numbers of increments from which noisy repeated readings were generated as per the approach previously developed (section and Fig. S2). The growth parameters were then re-estimated using the coupled VB-logK model and a classical VB-logK model, which estimates age using the traditional method. As the parameter is difficult to estimate from a tagging dataset, its value was fixed using the value obtained by fitting the otolith data to the coupled growth model.(E1)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Attributes of marginal posterior distributions from VB log K growth model coupled with ageing error model fit to yellowfin otolith data.

https://doi.org/10.1371/journal.pone.0060886.t003

The mean squared error (MSE) was used to compare the marginal posterior distributions of each parameter with its actual value and evaluate the accuracy of models fits. Accuracy is related to the similarity between the marginal modes and the true values, while precision refers to the uncertainty around the modal value, i.e. the discrepancies between the samples generated by the MCMC simulation. The MSE is a measure of the average of the square of errors and was calculated as follows:where is the size of the MCMC sample, is the parameter estimate at iteration and is the actual value of the parameter.

In addition, the mean relative error (MRE) was used to give information on the bias in parameter estimates:

Finally, an age-length key was derived from both the coupled model and the classical model and was used to convert the catch-at-size into catch-at-age. The two conversions were then compared. For this example, we used 2010 IOTC catch data from French purse seiners in the Western Indian Ocean.

Simulation evaluation

The fit of the ageing model to the simulated data showed good agreement between the estimated and simulated age, as well as a strong ability to estimate the increment deposition rate . From the simulations conducted in the first model step, the marginal posterior distribution of , derived from MCMC outputs, followed a normal distribution with a mode and variance very close to the simulated distribution. The fixed simulated values were well within the 95% credibility interval of the posterior distributions (Fig. S3). The model showed difficulties in estimating the edge bias by systematically underestimating its variability.

Simulations that focus on the number of readings showed that the ageing error model provided accurate age estimates, with accuracy increasing with the number of readings (Fig. S4). The RMSE values were significantly lower when there were more readings of the same otolith (Wilcoxon: p-value0.05; Table S2). From five to three readings, the loss of accuracy was generally low, i.e. an average loss of 1%. However, the loss of accuracy became more important between three to two readings, suggesting that three readings was a good compromise between accuracy of estimates and reading cost.

Regardless of the number of readings, the ageing error model estimated age more accurately than the traditional method. The RMSE values obtained with the ageing error model were significantly lower than those obtained with the traditional method (Wilcoxon: p-value0.01; Fig. S5 and Table S3). The ageing error model and the intermediate ageing method provided age estimates as accurate. The RMSE values obtained with the ageing error model and the intermediate method did not differ significantly (Table S3), with the exception of second simulation dataset where the RMSE values significantly lower with the ageing error model (Wilcoxon: p-value0.01).

Testing the hypothesis of daily increment deposition

The ratio and the mean edge bias were precisely estimated in the regression model (Fig. S3). was normally distributed around 0.94 with little variability (Table 3). The value of one was not included in the Bayesian 95% credibility interval which resulted in a significant failure of the hypothesis of daily increment deposition. The reader tended to underestimate the fish ages. The edge bias was close to zero, with a lower variability than expected from the prior distribution, resulting in error values ranging from 1.6% for a six month-old fish to 0.16% for a five year-old fish (Table 3).

Yellowfin tuna growth

The model supported a two-stanza growth for the IO yellowfin with two distinct phases over the fish lifespan of the fish (Figure 2). The first stanza was characterized by relatively slow growth, which decreased gradually to a minimum of 1.75 cm.month for fish up to 1.85 years of age (63 cm ). In the second stanza, growth accelerated up to a maximum of 3.54 cm.month for fish up to 2.38 years of age (79 cm ) and then progressively decreased with size, becoming very slow when the size approached the asymptotic length.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Yellowfin growth curve as estimated from the VB log K model coupled with the ageing error model.

The solid line correspond to the mean growth curve, the dashe to the uncertainty around the mean curve and the points to observation data.

https://doi.org/10.1371/journal.pone.0060886.g002

The fit of the coupled growth model revealed strong correlations among some parameters. A positive correlation was found between and (0.8) and a negative correlation was found between these growth rate coefficients and the asymptotic length (−0.93 and −0.8 with and respectively; Table S4). These correlations resulted in some model instability and difficulty fitting the posterior distributions of the parameters, particularly in the absence of information from the prior distribution (Fig. S6). The expected asymptotic length was estimated at 151.8 cm , which was close to the maximum fork length of the data, 146.5 cm , but comparatively low to the mean asymptotic fork length estimated at about 173 cm from the catch of the purse seiners and longliners and to the maximum lengths of 200 cm that has been observed for yellowfin in the Indian Ocean. In contrast, the data provided more information for the first part of the growth curve (3 years), where the model fit was good. The uncertainty associated with the mean growth curve that was derived from marginal posterior modes increased as the fish aged.

Sensitivity of growth parameters to biased age estimates

Except for , the MRE values derived from the simulations indicated that estimates were more accurate for the coupled model than for the classical model (Table 4). For both models, was estimated with a negative bias of an average value of 19% and 10% for the coupled and classical models respectively. was estimated with a positive bias for the classical model. For the other parameters, the mean biases were less than 10%. The MSE values were high for the asymptotic length . The influence of the GEV prior distribution was greater for the asymptotic part of the curve than for the fit to the real data as this part had very few data points. This resulted in the asymptotic lengths being overestimated. For the other parameters, the MSE values were low (Table 4). The MSE values for the estimates of , and were higher for the coupled model than for the classical model indicating a greater variability in the estimates of the coupled model. These results suggested that the coupled model provided more accurate, but less precise, estimates than the classical model because it accounted for uncertainties associated with age estimates while the ages were fixed in the classical model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Features of both growth models fit to simulated data.

https://doi.org/10.1371/journal.pone.0060886.t004

The growth curve obtained from the coupled model fit coincided with the simulated data in the first stanza but tended to underestimate the growth between 3.5 and 6.5 years. The growth curve of the classical model consistently overestimated the simulated growth (Figure 3). Such deviations in growth patterns led to a significant divergence in the age-length key that was derived from both models. This divergence resulted in potentially major differences when converting catch-at-size into catch-at-age. As illustrated by the 2010 yellowfin catch data, the proportions-at-age by 10 cm class lengths were significantly different according to the age-length key, particularly for lengths of less than 100 cm (Figure 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Growth curves obtained from the fit to simulated data with the coupled growth model (A) and the classical growth model (B; up) and their use to convert the size frequencies from fishing catches into age frequencies (down).

Gray levels correspond to the age classes (quarter).

https://doi.org/10.1371/journal.pone.0060886.g003

Age estimation method

The traditional ageing method, which converts an average count of increments into age, depends on the assumption that each increment corresponds to one day [60]. Otherwise, the use of this method will lead to systematically biased age estimates. In the case of the European hake (Merlucccius merluccius), it has recently been demonstrated that the assumption of a daily deposition rate is not met. As such, although the use of the traditional method for this species is internationally accepted, its use has resulted in a significant overestimation of fish age with important implications for estimates of key biological parameters [61]. This example highlights that validation of the increment deposition rate for the entire age range of interest should be major prerequisite for any ageing study. Various age validation techniques exist including the release of tagged fish of a known age, nuclear or radiochemical dating, modal progression of catch-at-size data over time, and captive rearing [25]. For yellowfin, the assumption of daily deposition was tested through a mark-recapture experiment of wild fish with oxytetracycline (OTC) in the length interval from 49.7–131 cm . The result of this experiment suggested that the deposition rate estimated by the reader within this length interval is significantly different than the assumed single increment per day, with an increment-over-day ratio of 1∶0.939. The marginal posterior distribution revealed little variation in this ratio, suggesting that the reading method tended to systematically underestimate the actual fish ages. However, additionnal readings performed by new readers for some of the otolith data used in the present study however indicate a 1∶0.997 ratio consistent with a daily deposition rate as shown for the Pacific Ocean yellowfin [48]. As tuna otoliths are structurally complex, a high degree of interpretation is required reading when them, which can lead to a reading bias and may account for the differences in estimated deposition rates [62]; [63]. This highlights the need for inter-comparison exercises between research institutes, using reference otolith collections, so as to harmonise reading methods [23]. This first step of the model can be adapted to other species or to other age validation techniques.

It is also necessary to assess the precision of the ageing technique [25]. The main uncertainty in age estimates arises from discrepancies between repeated readings. The magnitude of this uncertainty is partly dependent on the proficiency of the reader [16]. Thus, when a calcified structure is read by various readers, our model must be adapted to take account of the probability of misinterpretation of each reader. In addition, for some otoliths, the interpretations of the first and marginal growth increments can generate further uncertainty, particularly as these rings may obliterate during a poor otolith preparation. These additional uncertainties were relfected in the model through the choice of stochastic error structures and informative prior distributions based on the expertise to compensate for the lack of information in the data. In this context, Bayesian analysis appears particularly adapted to include such information into the statistical inference procedure.

The model evaluation by simulation suggested that it was capable of estimating both the ageing biases and an accurate age regardless of the age of the fish. The mean and variance of the deposit bias were consistently well-estimated. The model performed less well in estimating the bias linked to the misinterpretation of the marginal rings. However, edge level errors are partly estimated by the overall reading error. For yellowfin, expert judgement suggests that the nucleus and edge biases could be very low, i.e. less than five increments, resulting in error ranging from 4% for a six month-old fish to 0.4% for a five year-old fish. As such, these biases would become negligible when the age in days values are converted to age in year values. Thus, the intermediate method in which the average number of increments is corrected for the actual deposit frequency provided age estimates as accurately as the ageing error model. However, the range of uncertainty is fully ignored which may have serious implications on the final outcome of the analysis [43].

Growth of Indian Ocean yellowfin

The method developed in this study appears suitable for modeling the growth of Indian Ocean yellowfin. The fit of the growth model appeared to be adequate in the first stanza. In this first phase, the expected first growth rate, 2.11 cm.month for fish from 22 to 63 cm , and length from which growth accelerates, 63 cm , were consistent with several growth studies for the same population. Growth rates of small yellowfin have been estimated to be between 1.3 and 2.9 cm.month [30]; [63]; [31]; [65]; [32] with growth acceleration occurring around 62 to 66 cm [31], 70 cm [65], 56 to 66 cm [32] and from 60 cm [33]. However, in the absence of older fish, the model seemed to underestimate the expected asymptotic length leading to an apparent underestimation of yellowfin growth over 2.7 years of age (95.72 cm ). Increasing the uncertainty in the second stanza reflected greater individual variability in growth. A primary source of variability in growth rates may arise from sexual dimorphism, characterized by faster growth and a larger asymptotic length for males than females [66]; [36]. According to Wild [60], in the eastern Pacific, yellowfin females grow faster than males until 94.9 cm , i.e. at about 2 years old, and then the trend reverses. Several authors have also shown that in IO yellowfin, males are largely dominant above 145 cm [67]; [68] which is close to the modal value of the asymptotic length estimated in this study at 146.5 cm .

Estimating the growth parameters with the growth model coupled with the ageing error model was more accurate, but less precise, than with the classical growth model. This latter used the traditional ageing method which underestimated the “actual” age of fish and resulted in a consistent overestimation in the growth curve. In contrast, taking into account uncertainties associated with age estimates led to a growth curve that fitted very closely to the simulated data in the first stanza. Such differences greatly influenced the age distributions of the fisheries catches. In addition, by disregarding the ageing errors, the classical model was overconfident in determining uncertainties of growth estimates. It is critical to represent the full range of uncertainty to adequately evaluate management alternatives and the expected consequences of decision-making [69]; [70]. Developed within a hierarchical Bayesian framework, our method appears particularly suited to quantifying the uncertainties resulting from otolith reading or observational errors and their propagating in the estimates of biological parameters.

However, there are some technical difficulties associated with the use of Bayesian inferences. One of the main criticisms is the need to specify prior distributions, which can have some influence on the parameter posterior distributions [71; 38]. Commonly, sensitivity analyses are performed to evaluate the effect of prior change on the posterior. As increasing data quality and quantity reduces the prior influence. Therefore, these sensitivity analyses are above all an assessment tool for appreciate the amount of information in the data [40]. In this study, extensive work was undertaken on prior specification and weakly informative distributions were assigned to the parameters for which there was no established scientific understanding. This was done to obtain estimates that were more consistent with the biology of yellowfin.

The general approach developed in this study can be applied to other species and can be adapted to any structure that produces periodic growth increments in response to various ecological issues.