Introduction

Age determination is a key imput for stock assessment of aquatic living resources, specifically in populations undergoing fishery exploitation (Siskey et al. 2016; Quist & Isermann 2017). By knowing fish age from wild individuals, it is possible construct the age-structure of population, a crucial part for stock assessment models(Botsford et al. 2011; Fowler & Ling 2010).

The Chilean Jack Macherel(CJM) is a pelagic fish the widest distribution area of this genus and it is considered as a species with basin-scale distribution range, being observed mainly in the whole sub-tropical area of the South Pacific Ocean, also it is one of the biggest purse-seine fishery for small pelagic fishes in the world (Arcos et al 2001, Cardenas et al. 2005, Cardenas et al.2009, Gerlotto et al, 2012). The demographic strategy and the impotance to the fishery of this species increase the risk of overfishing, makes an international regulation agreement to be a vital importance. for this reason, the last decade the CJM has therefore become a concern for the South Pacific Region Fisheries Management Organization (SPRFMO).

For managemment purposes, the individual age was estimated since 1975 using standard sclerocronologic techniques of whole and cross-sectioned otolith in which an opaque and traslucient ring is identifying and counting as an annual increase or annulus (Panfili et al 2001). However, CJM have been shown to have particularly complicated ring structures because in many cases one annulus may contain double or multiple translucent areas(Cerna et al 2016). Wherefore, the potential of an inaccurate ageing is high and may result in a bias of age estimation, so the validation of the absolute annual growth is vital(Campana et al 1992).

The annual periodicity formation of hyaline and opaque growth zones was reported the first time for Aguayo et. al. (1981). This estimations appears to be consistent with the annuli structure of the other species of the same genus, however, discrimination between a real and a false annual growth is not an easy task and uncertainty remains in the interpretation of three first growth annuli.

Regardless of the difficulties in the annuli identificacion, the accuracy of the CJM otolith ageing had been reported acceptable precision (reproducibility) at least Age 10, but uncertainties in the age determination method still persist. Therefore, to keep accuracy and consistency in ageing, in the last years the age of CJM has been indirectly validated with sereval approaches, this studies demonstrated high growth in the first years of life which would suggest an overestimation of the ageing criteria using since 1975 (Cerna et al 2016, Araya et al 2020). This studies would be supporting the suspicion of age overestimation, and their most important results are synthesized in:

The above results allow us to confirm that the second growth mark finished its formation very close to the end of the first year of life and the second one very far from completing a new annual cycle. Therefore, The actual ageing criteria has been overestimated two years at least (Araya 2003; Araya 2020; Cerna et al 2016).

In accordance with the last evidence, it is necessary to evaluate the impact of the overestimation to which the age estimation has been subject. For which we compare the age estimation by conventional methods(historic criteria) with the new criteria resulting from the new results obtained in the validation(validation criteria)

Methods

The change of the age structure was performed with the historic data base wich contains individual age estimations of 133.943 individuals analysed since 1990 until 2018. For each individual contained in the database, the geographic position, date of capture, fork length, otolith size and radius in mm of each annulli recorded have been recorded.

To correct the overestimation of age, two radii will be subtracted from the number of radii of each specimen registered in the historic database. The subtraction of the radii registered will be done according to the following criteria:

These was the most accurate way to corrected the historical ages series accord the new ageing criteria, whose modification was automated with R software as follow:

Libraries

library(egg)
## Warning: package 'egg' was built under R version 4.0.2
library(gdata)
library(dplyr)
library(tidyr)
library(FSA) 
library(utils)

Functions

Marginal increment

IM <- function(x){
  
  im <- rep(NA, length(auxAge2$RADIO_TOTAL))
  rt <- auxAge2$RADIO_TOTAL 
  df <- auxAge2
  
  for(i in 1:dim(df)[1]){
    
    RT <- numeric(); RT = rt[i]
    temp <- numeric(); temp = df[i,c(9:32)] 
    
    if(all(is.na(as.numeric(temp))))
      
      im[i] <- 0.249
    
    else
      
      X1 = numeric(); X1 = temp[which(is.na(as.numeric(temp)))[1]-2]
      X2 = numeric(); X2 = temp[which(is.na(as.numeric(temp)))[1]-1]
    
    if(length(X1) > 0 & length(X2) > 0)
      
      if(as.numeric(X1) > 0 & as.numeric(X2) > 0)
        
        im[i] <- (RT-as.numeric(X2))/(as.numeric(X2)-as.numeric(X1))
    
    else
      
      im[i] <- 0.249
      
    else
      
      im[i] <- 0.249 
    
  }
  
  im[which(im > 1)] <- 0.249
  im[which(im <= 0)] <- 0.249
  
  return(im) 
  
}  

Age

EDAD <- function(x,y){
  
  z <- numeric(length = length(x))
  
  for(i in 1:length(x)){

    if(y[i] > 2) 
      
      z[i] <- x[i] + 1
    
    else 
      
      z[i] <- x[i]
  }
  
  return(z)
  
}

Age group

GRUPO_EDAD <- function(x){
  
  m <- as.numeric(Age$MES)
  b <- as.numeric(Age$BORDE)
  im <- as.numeric(Age$Im)
  ge <- rep(NA, length(Age$Age))
  
  for(i in 1:length(x)){
    
    if(m[i] <= 5 & b[i] <= 2) 
      
      ge[i] <-  x[i] + 1
    
    else if(m[i] == 6 & b[i] == 1 & im[i] >= 0.250)
      
      ge[i] <-  x[i] + 1
    
    else if(m[i] == 6 & b[i] == 2)
      
      ge[i] <-  x[i] + 1
    
    else if(m[i] == 7 & b[i] == 1 & im[i] >= 0.350)
      
      ge[i] <-  x[i] + 1
    
    else if(m[i] == 7 & b[i] == 2) 
      
      ge[i] <-  x[i] + 1
    
    else if(m[i] == 8 & b[i] %in% c(1,2) & im[i] >= 0.500)
      
      ge[i] <- x[i] + 1 
    
    else if(m[i] == 9 & b[i] == 2 & im[i] >= 0.700)
      
      ge[i] <- x[i] + 1 
    
    else
      
      ge[i] = x[i]
    
  }
  
  return(ge)
  
}

Length weight relationship

AB <- function(x){
  ab <- x %>%
    filter(!is.na(LONGITUD_ESPECIMEN),
           !is.na(PESO_ESPECIMEN)) %>%
    group_by(ANO,MACROZONA,TRIMESTRE) %>% 
    mutate(L_LH = log(LONGITUD_ESPECIMEN),
           L_W = log(PESO_ESPECIMEN)) %>%
    summarise(Slope = lm(L_W ~ L_LH)$coefficients[2],
              Intercept = exp(lm(L_W ~ L_LH)$coefficients[1]),
              Correlation = cor(L_W,L_LH),
              N = n()) %>%
    mutate(bmax = round(Slope,1)) %>%
    left_join(PR,by = c("bmax")) %>%
    rename(tempbmax1 = bmax,
           tempA = temp1, 
           tempC = temp2 ) %>%
    mutate(bmax= round(round(Slope,1) - 0.1,1)) %>%
    left_join(PR,by = c("bmax")) %>%
    rename(tempbmax2 = bmax,
           tempB = temp1, 
           tempD = temp2) %>%
    mutate(a1 = as.numeric(tempA - ((tempbmax1 - Slope) * (tempA - tempB)) / (tempbmax1 - tempbmax2)),
           a2 = as.numeric(tempC - ((tempbmax1 - Slope) * (tempC - tempD)) / (tempbmax1 - tempbmax2))) %>%
    dplyr::select(-starts_with('temp'))
  
  return(ab)
  
}

Catch at age

CATCH_AGE <- function(x){

  auxAge4 <- x
  
# LECTURAS DE EDAD CON TALLAS SIN ESTIMACION PARA CADA MACROZONA POR TRIMESTRE
  
  CaptAgeLength <- auxAge4 %>% 
    group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE,LONGITUD_ESPECIMEN,GRUPO_EDAD) %>%
    summarise(FREC = n()) %>%
    mutate(PROP_GRUPO = FREC/sum(FREC)) %>%
    left_join(Length, by = c("CRITERIO","ANO","MACROZONA","TRIMESTRE","LONGITUD_ESPECIMEN")) %>%
    mutate(CAPT = PROP_GRUPO*CAPT) %>%
    ungroup()
  
# CAPTURAS A LA EDAD Y TALLA POR CADA MACROZONA Y TRIMESTRE

  auxCatchAge <- CaptAgeLength %>% 
    dplyr::select(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD,CAPT) %>%
    group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD) %>%
    summarise(CAPT = sum(CAPT,na.rm=T)) %>%
    ungroup() %>%
    group_by(CRITERIO,ANO,MACROZONA, TRIMESTRE) %>%
    mutate(PORC =  CAPT/sum(CAPT,na.rm=T) *100)
  
  auxCatchAge1 <- CaptAgeLength %>%
    mutate(tmed = LONGITUD_ESPECIMEN * CAPT, 
           varl = (LONGITUD_ESPECIMEN^2) * CAPT) %>%
    filter(CAPT != 0) %>%
    group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD) %>%
    summarise(VARL = sum(varl, na.rm = T),
              TMED = sum(tmed, na.rm = T))
  
  auxCatchAge2 <- CaptAgeLength %>%
    filter(CAPT != 0) %>%
    group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD) %>%
    summarise(L_PR = weighted.mean(LONGITUD_ESPECIMEN, CAPT))

# CAPTURAS A LA EDAD
  
  CatchAge <- auxCatchAge %>%  
    left_join(auxCatchAge2,  by = c("CRITERIO","ANO","MACROZONA","TRIMESTRE","GRUPO_EDAD"))%>%
    left_join(ab, by = c("CRITERIO","ANO","MACROZONA","TRIMESTRE"))%>%
    dplyr::select(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD,CAPT,PORC,L_PR,Intercept,
                  Slope, Correlation, N,a1,a2)%>%
    left_join(auxCatchAge1,by = c("CRITERIO","ANO","MACROZONA","TRIMESTRE","GRUPO_EDAD"))%>%
    ungroup()%>%
    mutate(VAR =    abs((VARL-(L_PR*TMED))/(CAPT-1)),
           P_PR = (Intercept * ((L_PR^Slope) + 
                                  a1 * (L_PR^(Slope-2)) * VAR +
                                  a2 * (L_PR^(Slope-4)) * VAR^2))/1000,
           CAPT_TON = (CAPT* P_PR) / 1000) %>%
    group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE)%>%
    mutate(PORC_TON = (CAPT_TON / sum(CAPT_TON)) * 100)%>%
    ungroup()%>%
    dplyr::select(CRITERIO,ANO,MACROZONA,TRIMESTRE,GRUPO_EDAD,
                  CAPT,PORC,L_PR,VAR,P_PR,CAPT_TON, 
                  PORC_TON)%>% 
    ungroup() %>% 
    filter(CAPT != 0) %>%
    mutate(GE = "GE")%>% 
    unite(GRUPO_EDAD, GE,GRUPO_EDAD) %>%
    mutate(ANO = as.numeric(ANO)) %>%
    separate(GRUPO_EDAD, into = c("G","GRUPO_EDAD"), sep ="_") %>%
    mutate(GRUPO_EDAD = as.numeric(GRUPO_EDAD)) %>%
    dplyr::select(-G)

  return(list(CatchAge, CaptAgeLength))

} 

Inputs

# AJUSTE DE PENDIENTE
PR = tbl_df(read.csv2("PIENNAR_RICKER.csv", sep = ";"))

# BASE DE DATOS HISTORICA DE JUREL
RefHist <- tbl_df(read.csv2("REFERENCIA_HISTORICA_GRUPO_EDAD_JUREL.csv", sep =";")) 

# CAPTURAS EN NUMERO
CAPT_TALLA <- tbl_df(read.csv2("TALLAS_TRIMESTRALES_1990_2018.csv",sep = ";")) 

# MUESTREO BIOLOGICO
BIO_HIST <- tbl_df(read.csv2("BIOLOGICO_JUREL_1997_2018.csv",sep = ";")) %>%
  dplyr::select(PUERTO_RECALADA,ANO,MES,LONGITUD_ESPECIMEN,PESO_ESPECIMEN)

# PARAMETROS a y b DE RELACION LONGITUD PESO
ab_1990_2000 <- tbl_df(read.csv2("AB_1990_1996.csv", sep = ";"))

# ESTIMACION DE EDAD(LECTURAS DE OTOLITOS)
auxAge0 <- tbl_df(read.csv2("EDAD_1975_2018_SEGUIMIENTO.csv", sep = ";"))  %>% 
  filter(!is.na(RADIO_TOTAL))

Catch at length

# CAPTURA EN NUMERO POR ESTRATO DE TALLA
auxLength <- CAPT_TALLA %>%
  group_by(ANO,ZONA,LONGITUD_ESPECIMEN) %>%
  gather(TRIMESTRE,CAPT, starts_with("T")) %>%
  separate(TRIMESTRE, into = c("X", "TRIMESTRE"), sep = "(?<=[A-Za-z])(?=[0-9])") %>%
  mutate(TRIMESTRE = as.numeric(TRIMESTRE),
         CRITERIO = 1,
         MACROZONA = case_when(ZONA == 1 ~ 1,
                               ZONA >= 2 ~ 2,
                               TRUE ~ NA_real_)) %>% # CRITERIO HISTORICO
  ungroup() %>%
  dplyr::select(-ZONA,-X)

# CREANDO UN VECTOR IDENTICO A "Length" PARA SER UTILIZADO CRITERIO = 2

Length <- auxLength %>%
  mutate(CRITERIO = 2) %>%
  bind_rows(auxLength)

Age length key

# LECTURA OTOLITOS CRITERIO HISTORICO = 1

auxAge1 <- auxAge0 %>% 
  filter(!is.na(RADIO_TOTAL),
         !is.na(LONGITUD_ESPECIMEN),
         !is.na(BORDE), 
          ANO %in% c(1990:2018)) %>%
  group_by(ANO, LONGITUD_ESPECIMEN) %>%    
  
##########################################################
#  mutate(FREC=n()) %>%                                  # TOMA MUESTRA PROPORCIONAL
#  sample_frac(0.01, weight=FREC) %>%                    # PARA PODER PROBAR EL CODIGO
#  ungroup() %>%                                         # CON UNA FRACCION DE LA BASE DE DATOS
##########################################################
 
  mutate(ID = rep(1:length(RADIO_TOTAL)),
        CRITERIO = 1) %>% # CRITERIO HISTORICO DE LECTURA
  dplyr::select(ID,CRITERIO,ZONA,MES,ANO,LONGITUD_ESPECIMEN,BORDE,RADIO_TOTAL,
                R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15,R16,R17,R18,
                R19,R20,R21,R22,R23,R24)

# RESTANDO RADIO 1 y 3 = 2

auxAge2 <-  auxAge1 %>%
  group_by(ID,ZONA,MES,ANO) %>%
  gather(R,RADIO,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15,R16,R17,R18,
         R19,R20,R21,R22,R23,R24) %>%
  separate(R, into = c("R", "N"), sep = "(?<=[A-Za-z])(?=[0-9])") %>%
  filter(N != 1,N !=3) %>%
  mutate(N = case_when(N == 2 ~ 1,   N == 4 ~ 2,   N == 5 ~ 3, 
                       N == 6 ~ 4,   N == 7 ~ 5,   N == 8 ~ 6, 
                       N == 9 ~ 7,   N == 10 ~ 8,  N == 11 ~ 9,
                       N == 12 ~ 10, N == 13 ~ 11, N == 14 ~ 12, 
                       N == 15 ~ 13, N == 16 ~ 14, N == 17 ~ 15,
                       N == 18 ~ 16, N == 19 ~ 17, N == 20 ~ 18, 
                       N == 21 ~ 19, N == 22 ~ 20, N == 23 ~ 21, 
                       N == 24 ~ 22, TRUE ~ NA_real_)) %>%
  unite(GE,R,N) %>%
  unite(COR, ID,ZONA,MES,ANO,LONGITUD_ESPECIMEN,BORDE,RADIO_TOTAL) %>%
  spread(GE, RADIO) %>%
  separate(COR,into= c("ID","ZONA","MES","ANO","LONGITUD_ESPECIMEN","BORDE","RADIO_TOTAL"), 
           convert =TRUE) %>%
  mutate(R_23 = NA, R_24 = NA, CRITERIO = 2) %>% 
  arrange(ID) %>%
  rename(R1 = R_1, R2 = R_2, R3 = R_3, R4 = R_4, R5 = R_5, R6 = R_6, R7 = R_7,
         R8 = R_8, R9 = R_9 ,R10 = R_10, R11 = R_11, R12 = R_12, R13 = R_13,
         R14 = R_14, R15 = R_15, R16 = R_16, R17 = R_17, R18 = R_18, R19 = R_19,
         R20 = R_20, R21 = R_21, R22 = R_22,R23 = R_23, R24 = R_24) %>%
  bind_rows(auxAge1) %>%
  dplyr::select(ID,CRITERIO,ZONA,MES,ANO,LONGITUD_ESPECIMEN,BORDE,RADIO_TOTAL,
                R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15,R16,R17,R18,
                R19,R20,R21,R22,R23,R24) 

Age

To estimate the biological age, the number of radii recorded and the nature of the otolith edge are considered. If the otolith edge is hyaline (values 3 or 4 of the variable BORDE), one is added to the number of radii counted.

# VECTOR 1: EDAD

Age <- auxAge2 %>%
  mutate(MACROZONA = case_when(ZONA <= 3 ~ 1,
                               ZONA >= 4 ~ 2,
                               TRUE ~ NA_real_),
         TRIMESTRE = case_when(MES == 1 | MES == 2 | MES == 3 ~ 1,
                               MES == 4 | MES == 5 | MES == 6 ~ 2,
                               MES == 7 | MES == 8 | MES == 9 ~ 3,
                               MES == 10 | MES == 11 | MES == 12 ~ 4,
                               TRUE ~ NA_real_),
         NumRings =  apply(auxAge2[,9:32],1, function(x){length(x[!is.na(x)])}),
         Im = IM(RADIO_TOTAL),
         Age = EDAD(NumRings,BORDE))

Marginal increment

As a rule, if the marginal increment is translucent it is not counted because it is not fully formed. in the for cycle above, we calculate the proportion of the traslucent increment i in relation of the inncrement i-1.

Age group

The age group is the number of calendar years after the date of capture, for this purpose each fish will be assigned depends on the year in wich it was spawned and on the date of capture, simplifying the identification of annual class per fish.

This estimate is calculated with the number of rings observed in the otolith, the type of edge of the otolith and an arbitrary date of birth. In Chile, the first date of birth is taken as a convention on January 1.

Age <- Age %>% 
  mutate(GRUPO_EDAD = GRUPO_EDAD(Age))

# IDENTIFICANDO TALLAS SIN ESTIMACION DE EDAD

auxAge3 <- Age %>% # CT
  group_by(CRITERIO,ANO,MACROZONA,TRIMESTRE,LONGITUD_ESPECIMEN,GRUPO_EDAD) %>%
  summarise(FREC = n()) %>%
  ungroup()%>%
  mutate(ANO = as.integer(ANO),
         LONGITUD_ESPECIMEN = as.integer(LONGITUD_ESPECIMEN)) %>%
  right_join(Length, by = c("CRITERIO","ANO","MACROZONA","TRIMESTRE","LONGITUD_ESPECIMEN"))

# COMPLETANDO GRUPO DE EDAD EN TALLAS SIN ESTIMACION DE EDAD

auxAge4 <- auxAge3 %>% # DIF
  filter(is.na(FREC) & !is.na(CAPT)) %>%
  dplyr::select(-GRUPO_EDAD,-FREC,-CAPT) %>%
  left_join(RefHist, by = c("CRITERIO","LONGITUD_ESPECIMEN")) %>%
  dplyr::select(CRITERIO, ANO,MACROZONA,TRIMESTRE, LONGITUD_ESPECIMEN, GRUPO_EDAD)

Age <- Age %>% # auxC1
  dplyr::select(CRITERIO, ANO, MACROZONA, TRIMESTRE, LONGITUD_ESPECIMEN, GRUPO_EDAD)%>%
  bind_rows(auxAge4) 

Length weight parameters

The length-weight relationship is widely used in fisheries to estimate weight from the length of an individual and also to estimate condition indices. The most frequently used expression for this relationship corresponds to the allometric equation where the weight is expressed as a function of length, and its parameters are estimated by linear regression of log-transformed data. Since variability in weight usually increases with length, this transformation has the advantage that it tends to stabilize the weight variance, but introduces a bias factor in back-transformed predictions which requires correction.

The weight/length model is:

\[W_i = a*L_i^b\] We can derive a simple linear regression model applying a logarithmic transformation:

\[ln(W_i) = ln(a)+b*ln(L_i)\]

Where Wi represents the weight and Li the length of the i-th fish; a and b correspond to the intercept and slope, respectively.

BIO <- BIO_HIST %>%
  mutate(MACROZONA = case_when(PUERTO_RECALADA <= 6 ~ 1,
                               PUERTO_RECALADA >= 10 ~ 2,
                               TRUE ~ NA_real_),
         TRIMESTRE = case_when(MES == 1 | MES == 2 | MES == 3 ~ 1,
                               MES == 4 | MES == 5 | MES == 6 ~ 2,
                               MES == 7 | MES == 8 | MES == 9 ~ 3,
                               MES == 10 | MES == 11 | MES == 12 ~ 4,
                               TRUE ~ NA_real_))
auxab <- AB(BIO) %>%
  mutate(CRITERIO = 1)

ab <- auxab %>%
  mutate(CRITERIO = 2) %>%
  bind_rows(auxab,ab_1990_2000) %>%
  arrange(CRITERIO,ANO,MACROZONA,TRIMESTRE) %>%
  filter(ANO %in% c(1990:2000,2002,2005:2018))

# Borrar objetos auxiliares
rm(list=ls(pattern="aux"))

Catch at age

The age structure of the capture was performed with the quarterly age keys constructed for each fleet in order to extrapolate the catch by length stratum to catch by age group.

Once a key is available, samples of fish that were only measured for length can be distributed to age groups accordingly.

This is possible because it assumes that the sample of aged fish and the sample of fish measured for length are simple random samples from the same population, thus the probability that a fish is of a particular age, given its length is the same in both samples follows.

The individuals present in each length interval (Nj) are assigned to different ages according to a size-age key.

The number of individuals belonging to each age group according to size interval is estimated as:

\[N_{ij} = q_{ij}/N_j\] where:

The abundance by age group is obtained by:

\[N_i = \sum_{i=1}^{n} N_{ij} \]

where:

CatchAge <- CATCH_AGE(Age)

Mean length

The average length for each age group is expressed as

\[ \bar{p_{i}} = \frac{\sum_{i=1}^{n} l_{j}n_{ij}}{ni} \]

Donde:

Variance

The variance of the mean length is expressed as

\[ P_i = \sum_{j=1} p_jq_{ij} \]

\[\sigma_{p_j} = \sum_{j=1}^{L} (\frac{p_j^{2}q_{ij}(1-q_{ij})}{n_{j}-1}+ \frac{p_{j}(q_{ij}-P_{i})^{2}}{N}) \]

Donde:

  • \(n_{j}\) : Estimated number of individuals of length j
  • N : sample Size
  • \(p_{j}\) : Sample size for each length stratum j
  • \(q_{ij}\) : Probability of individuals of length j belonging to an age group i

Mean weight

Assuming that the length is a normal random variable the expected function value of \(W_L\), will be defined as:

\[E(W) = a(\bar{p_{i}}^b+a_{1}\bar{p_{i}}^{b-2}\sigma^2+a_{2}\bar{p_{i}}^{b-4}\sigma^4)\]

Where: