Modeling and analysis of the recruitment of Peruvian anchovy (Engraulis ringens) between 1961 and 2009

R Oliveros–Ramos^1*, C Peña²

¹ Centro de Investigaciones en Modelado Oceanográfico y Biológico Pesquero, Instituto del Mar del Perú (IMARPE), Esquina Gamarra y General Valle s/n, Chucuito, Callao, Perú.* Corresponding author. E–mail: roliveros@imarpe.pe

² Unidad de Investigaciones en Dinámica Poblacional, Instituto del Mar del Perú (IMARPE), Esquina Gamarra y General Valle s/n, Chucuito, Callao, Perú.

Received September 2010 ]]> Accepted May 2011

RESUMEN

El objetivo principal de este trabajo es presentar un nuevo modelo para estimar el reclutamiento del stock norte–centro de la anchoveta peruana (Engraulis ringens) y usarlo para estimar series mensuales de reclutas y prereclutas entre 1961 y 2009. El modelo está basado en longitudes y tiene una resolución mensual, permitiendo modelar la variabilidad en la abundancia y la estacionalidad del reclutamiento a partir de información disponible de la pesquería y cruceros de evaluación científica. La población de anchoveta se modela con dos reclutamientos por año, parametrizando el reclutamiento por su abundancia en número, longitud media y dispersión de la longitud media al reclutamiento. Esta forma de modelar el reclutamiento permite representar mejor la continuidad temporal del proceso de reclutamiento que un enfoque basado sólo en el número de individuos reclutados. Al analizar las series mensuales de prereclutas (4.0–7.5 cm) y reclutas (8.0–11.5 cm) obtenidas se identificaron dos cambios de fase en la dinámica del reclutamiento (1971 y 1991), hallándose además que los patrones estacionales del reclutamiento obtenidos por el modelo son coherentes con los patrones descritos por otros estudios biológicos. Se discuten las ventajas y limitaciones del enfoque utilizado en el marco del estudio del proceso de reclutamiento de la anchoveta peruana y la integración de información reproductiva en los modelos de evaluación de stocks.

Palabras clave: reclutamiento, anchoveta peruana, series de tiempo, algoritmos evolutivos.

ABSTRACT

The principal objective of this work is to present a new model for the assessment of recruitment of the northern–central stock of Peruvian anchovy (Engraulis ringens), and use it to estímate monthly time series of recruits and pre–recruits between 1961 and 2009. The model is length–based and has monthly temporal resolution, allowing the variability in abundance and seasonality of recruitment to be modeled based on fishery and scientific survey information. The anchovy population is modeled with two recruitments per year, parameterizing recruitment using abundance, mean length, and dispersion of the mean length at recruitment. This way of modeling recruitment enables a better representation of the temporal continuity of the recruitment process in comparison to an approach that only considers the abundance of recruits. The monthly time series of pre–recruits (4.0–7.5 cm) and recruits (8.0–11.5 cm) obtained were analyzed, and two regime shifts in recruitment dynamics (1971 and 1991) were identified. Seasonal patterns given by the model are consistent with those obtained by independent direct studies. We discuss the advantages and limitations of this approach within the framework of the study of the recruitment process and the integration of reproductive data in stock assessment models.

Key words: Peruvian anchovy, recruitment time series, evolutionary algorithms.

La principal pesquería del ecosistema de la Corriente de Humboldt es la de la anchoveta peruana Engraulis ringens, con desembarques que superan los 5 millones de toneladas por año y que alcanzaron hasta 2.5 millones de toneladas en noviembre de 2005 (Ñiquen y Freón 2006). Esta pesquería es altamente variable a diversas escalas temporales y sus fluctuaciones tienen grandes impactos sociales y económicos (Agüero 1987, Aranda 2009). Por lo tanto, en este sentido es particularmente importante poder entender el efecto de la variabilidad temporal de la población de esta especie. En particular en Perú, el manejo de la anchoveta incluye vedas reproductivas (Perea et al. este número) por lo que es especialmente importante poder identificar cambios en la estacionalidad de los momentos de máximo desove y reclutamiento, a fin de poder proteger a los reproductores en sus momentos de mayor actividad reproductiva, así como a los organismos juveniles.

La pesca intensiva y la variabilidad ambiental inducen cambios en la estructura y el funcionamiento de los ecosistemas marinos, y es difícil predecir la respuesta de los ecosistemas a dichas perturbaciones debido a que tal respuesta resulta de complejas interacciones entre los diversos componentes del ecosistema. El clima y la pesca afectan el ecosistema en diferentes formas, como por ejemplo la luz sobre el crecimiento del fitoplancton, la temperatura en la distribución espacial del zooplancton (Beaugrand et al. 2002), la pesca en la distribución de tallas de las poblaciones (Shin et al. 2005), y la concentración de oxígeno (Bertrand et al. 2010) y las corrientes (Mullon et al. 2003) en la distribución y el éxito del reclutamiento de los pelágicos pequeños, respectivamente. La anchoveta peruana es especialmente sensible a estos cambios (Pauly y Tsukayama 1987) que influyen en el proceso de reclutamiento, que a su vez involucra una cadena de procesos y estadios naturales que incluyen desove, fertilización, desarrollo de las fases de huevo, larvas, crecimiento y supervivencia de prereclutas y reclutas (Bakun y Broad 2001).

Los reclutas son los individuos que ingresan por primera vez a la fase explotable del stock y comienzan a tener un efecto decisivo en la dinámica de su población (Csirke 1989); por lo tanto, la estimación del reclutamiento es muy importante en la evaluación de los stocks ya que influye significativamente en las decisiones de manejo. Así, Gulland (1968) manifestó que el reclutamiento modula cambios en la actividad pesquera, por lo que desde los inicios de la pesquería de la anchoveta se han realizado de forma permanente estudios sobre su reclutamiento. Estimaciones indirectas a través de índices de concentración del reclutamiento indicaron que la magnitud de éste dependía tanto de la biomasa total de la población desovante como del grado de concentración que alcanza la población ante la ampliación o reducción de su area normal de distribución (Csirke 1979), que a su vez varía por influencia del medio ambiente. Otros estudios han analizado el efecto de factores ambientales en la magnitud del reclutamiento en el norte del sistema de la Corriente de Humboldt (Cahuín et al. 2009), aunque estimando el reclutamiento a partir de un análisis de población virtual (APV) que presenta la limitación de trabajar con estimaciones anuales que enmascaran los efectos de los dos periodos de reclutamiento anuales de la anchoveta peruana (Jordán 1976), y sólo explora la variabilidad del reclutamiento en términos de su abundancia.

Durante los inicios de la pesquería de la anchoveta a inicios de la década de los 60, la temporada de reclutamiento a la pesquería ocurría en el otoño austral, especialmente de abril a junio (Satersdal y Valdivia 1964). Luego de El Niño 1972–1973 este patrón cambió de manera que el reclutamiento iniciaba a partir de noviembre/diciembre y alcanzaba sus valores más altos en febrero/marzo (Jordán 1976). Así es como la respuesta a los cambios ambientales puede modificar también la estacionalidad del reclutamiento, ya sea por el desplazamiento de los momentos de máximo desove o la ampliación de éstos ante cambios repentinos, como los que forman parte del fenómeno El Niño, o variaciones interdecadales en la productividad del ecosistema (Ñiquen et al. 1999, Alheit y Ñiquen 2004, Ñiquen y Bouchon 2004). En general, la variabilidad del reclutamiento sucede como una respuesta más o menos inmediata a las condiciones imperantes durante o inmediatamente después del desove (Csirke 1979), por lo que estos cambios pueden observarse en la distribución mensual de tallas de la población, pero se pierden al agregar los datos sobre una base anual y luego de aplicarles las claves talla–edad, lo cual es común en modelos basados en la edad como el APV. Considerando que la previsión del reclutamiento es esencial para el manejo sostenible de la pesquería de anchoveta, es importante contar con métodos alternativos que permitan estimar tanto la magnitud como la estacionalidad del reclutamiento, algo especialmente importante para un desovador parcial de vida corta como la anchoveta peruana.

Así, el objetivo principal de este trabajo es presentar un nuevo modelo para estimar el reclutamiento del stock norte–centro de la anchoveta peruana E. ringens y usarlo para estimar series mensuales de reclutas y prereclutas entre 1961 y 2009. El modelo está basado en longitudes, tiene una resolución mensual, y permite modelar la variabilidad en la abundancia y la estacionalidad del reclutamiento a partir de información disponible de la pesquería y cruceros de evaluación científica. Los resultados obtenidos se comparan con la información biológica obtenida por métodos directos y se discute el alcance y limitaciones del método propuesto, el cual se considera puede contribuir en un futuro a una integración de la información reproductiva en los modelos de evaluación de stocks.

MATERIALES Y MÉTODOS

Datos

Se utilizaron datos de capturas por tallas mensuales (clases de 0.5 cm, con un intervalo de tallas entre 4.0 y 20.0 cm) del stock norte–centro (4–16° S) de la anchoveta peruana entre enero de 1961 y diciembre de 2009 (IMARPE 2010). Se utilizó también información procedente de 47 cruceros de evaluación acústica de recursos pelágicos realizados anualmente entre 1963 y 2009, entre enero y abril por corresponder a los meses cuando la distribución del stock hace más confiables las estimaciones acústicas (Gutiérrez 2000). De éstos se obtuvieron estimaciones de biomasa total por métodos acústicos del stock norte–centro (Pauly y Tsukuyama 1987, Gutiérrez 2000, IMARPE 2010) y la estructura de tallas de la población estimada a partir de estudios biométricos.

La dinámica del stock de anchoveta peruana fue modelada considerando como variable de estado la estructura de cada cohorte, definida a partir del número de individuos (N) que la conforman, su longitud media (L) y su dispersión (S') medida por la desviación estándar de las longitudes. El paso de tiempo del modelo fue de un mes y se consideró una vida de tres años, con dos reclutamientos al año, para la anchoveta (Jordán 1976), por lo que en el modelo coexistían simultáneamente seis cohortes. Los procesos considerados fueron el reclutamiento, el crecimiento en longitud, la mortalidad natural y la mortalidad por pesca (bajo la forma de capturas por longitud). Estos procesos fueron descritos en dos submodelos: (1) reclutamiento y (2) dinámica de cohortes.

Reclutamiento

El modelo considera que el reclutamiento a la edad 0 (R₀) sucede dos veces al año, con puntos máximos que ocurren teóricamente el 1 de febrero y el 1 de agosto. La fecha teórica de los máximos de desove es usada internamente como punto de referencia para calcular la fecha real de nacimiento. Cada R₀ es modelado con una función tipo gaussiana (Schnute y Fournier, 1980) y se caracteriza por tres números: N₀, L₀ y N₀, donde N0 representa la magnitud del reclutamiento (en número de individuos), L₀ representa la longitud media de la cohorte al momento del desove y S₀ está relacionada con la dispersión de la cohorte en términos de su longitud, la cual toma en cuenta el desplazamiento temporal del máximo de desove con respecto a la fecha teórica supuesta. El número de reclutas de longitudes entre l1 y l2 será denotado por R₀([l1, l2]) y es calculado por

Para el modelo, se consideraron intervalos fijos de longitud de 0.5 cm, los cuales se describen a partir de sus marcas de clase. Así, considerando l₁ = l – 0.25 y l₂ = l + 0.25, denotaremos R₀ ([l₁, l₂]) por R¹₀. En principio, las tallas medias al momento del máximo teórico de desove pueden ser negativas, lo que indica que la fecha real del punto de máximo desove es posterior a la fecha teórica supuesta. De manera similar, tallas positivas al momento del máximo teórico indican que el máximo real ocurre antes de la fecha teórica supuesta. Los individuos con tallas negativas crecen según la ecuación de crecimiento, pero no son afectados por la mortalidad natural. A fin de considerar en la dinámica a los individuos que todavía no nacen, se consideró el intervalo de tallas efectivo entre –10.0 y 22.0 cm.

Dinámica de cohortes

Para la dinámica de las cohortes, se consideró el incremento en la longitud por efecto del crecimiento y la disminución del número de individuos por efecto de la mortalidad natural y de la pesca. Cada cohorte fue separada en clases de longitud de 0.5 cm utilizando la ecuación 1, y la mortalidad natural y el crecimiento se aplicaron a los individuos de cada intervalo de clase por la mitad del paso de tiempo de un mes (Δt = 15 días). El efecto de la mortalidad natural fue modelado por decaimiento exponencial (ecuación 2) y el crecimiento mediante el modelo de von Bertalanffy parametrizado según la ecuación de Brody (Schnute y Richards 2002) (ecuación 3) sobre cada una de las clases de longitud en las cuales se supuso que los individuos crecen y mueren de manera homogénea, escribiendo así:

en donde N^l_y y L^l_y (t) representan la abundancia y la longitud medias de los individuos en el intervalo de longitud con marca de clase l en el tiempo t, M es la mortalidad natural, L_∞ es la longitud asintótica de crecimiento y k es la tasa a la que se alcanza la longitud asintótica. Para la anchoveta peruana se supuso un valor de M = 0.7 año^–1, dentro del intervalo de 0.6 a 0.8 año^–1 usado por IMARPE (1974, 2010). Para la ecuación de crecimiento, se supusieron k = 0.07 mes^–1 y L_∞ = 21.5 cm (IMARPE 2010).

Hasta aquí se ha simulado el efecto del crecimiento y la mortalidad natural sobre la población por un periodo de 15 días, tras el cual se aplica el efecto de la pesca instantáneamente suponiendo que toda la captura es eliminada de la población en el día 15, es decir, a la mitad del paso de tiempo de la simulación. Para ello, se aplicó la ecuación 1 nuevamente a cada cohorte para agruparla en clases de longitud de 0.5 cm y se descontó el efecto de la pesca restando la captura por longitud mensual observada, que también se cuantifica en intervalos de clase de longitud de 0.5 cm con las mismas marcas de clase ponderando el efecto sobre cada cohorte según su abundancia relativa en cada clase de longitud. Cuando la captura observada en una clase de longitud era mayor al número de individuos presentes, se utilizó una captura equivalente al 99% de los individuos disponibles. Así, las capturas simuladas por el modelo son equivalentes a las capturas observadas excepto en los casos en que la población predicha era menor a la necesaria para soportar dicha captura.

Finalmente, se aplicó a cada cohorte el efecto de la mortalidad natural y el crecimiento por 15 días más, calculando nuevamente el número de individuos, longitud media y dispersión de la cohorte al final del mes mediante las ecuaciones 4, 5 y 6, concluyendo con la simulación de la dinámica poblacional para un paso de tiempo (30 días). El proceso se repitió hasta que una cohorte alcanzaba la edad de 36 meses, luego de lo cual se supuso que todos los peces de la cohorte morían de senectud y eran eliminados del modelo.

Para cada mes se calculó la abundancia de reclutas como el número de individuos comprendidos en los intervalos con marcas de clase de entre 8.0 y 11.5 cm, y la abundancia de prereclutas como el número de individuos comprendidos en los intervalos con marcas de clase entre 4.0 y 7.5 cm (Jordán 1976, IMARPE 1965).

Estimación de parámetros y calibración

El modelo se inicializa con la estructura poblacional dada por seis cohortes semestrales de entre 0 y 3 años de edad, cada una caracterizada por su abundancia, longitud media y dispersión. Así, para cada uno de los 98 reclutamientos semestrales y las seis cohortes iniciales se estimaron tres parámetros: N₀, L₀ y S₀, que resultaron en un total de 312 parámetros. Con la finalidad de prevenir el sobreajuste, se introdujeron restricciones al espacio de los parámetros estableciéndoles cotas superiores e inferiores. La abundancia de reclutas a la edad 0 (N₀) de cada nueva cohorte se acotó entre 1.0 × 10⁵ y 5.0 × 10¹², siendo este máximo equivalente a 1.5 veces el reclutamiento a la edad 0 estimado mediante la separación de cohortes de las estructuras de tallas de los cruceros y corrigiendo el efecto de la mortalidad natural hasta la fecha teórica de nacimiento. La longitud media a la edad teórica de desove (L₀) se acotó entre –4.5 y 4.0 cm, que corresponde a una diferencia de entre –3 y 3 meses con respecto a la fecha teórica de nacimiento, de manera que las fechas reales de nacimiento pudieran ser efectivamente cualquier día del año tomando en cuenta las dos cohortes. La dispersión (S₀) se acotó entre 0.25 y 3.0 cm, que implican que en un radio de una desviación estándar los individuos estén en un mismo intervalo de clase y dispersos sobre un intervalo de seis meses de crecimiento, respectivamente. Se estimaron dos funciones de error para biomasas y capturas, calculando (1) los errores entre las biomasas observadas durante los cruceros de evaluación acústica y las biomasas mensuales calculadas por el modelo y (2) los errores entre las capturas mensuales y las capturas estimadas por el modelo. Así, la biomasa total estimada por métodos acústicos es el principal calibrador del modelo, considerando que las capturas observadas son introducidas directamente como un forzante, por lo que la calibración sólo tiene el efecto de reproducir estructuras de tallas poblacionales que permitan reproducir las capturas observadas.

El problema de optimización no lineal asociado a minimizar ambas funciones de error simultáneamente fue solucionado utilizando el algoritmo evolutivo PBIL (population–based incremental learning) con una versión modificada del operador genético del individuo óptimo (Duboz 2010). Se utilizó un algoritmo evolutivo debido a la no linealidad del modelo y a la no derivabilidad introducida al truncar las capturas cuando la población no es capaz de soportarlas, lo que impide utilizar métodos de optimización convencionales basados en derivadas y no garantiza la convergencia a óptimos globales, condiciones en las que los algoritmos evolutivos se desempeñan mucho mejor (Andrzej y Stanislaw 2006). En el algoritmo, los "individuos" corresponden a las diferentes posibles soluciones (i.e., combinaciones de parámetros del modelo) para el problema de minimizar los errores entre los datos observados y los datos simulados. En cada iteración, el algoritmo calcula un "individuo óptimo" que reune el conjunto de parámetros que resulta de la combinación de los "individuos" de cada generación, y que resuelven mejor la minimización de los errores en la biomasa o en las capturas. Dicho "individuo óptimo" se calcula ponderando cada parámetro de acuerdo a su variabilidad durante el transcurso de la calibración. Con la finalidad de considerar la diferencia de escala entre los distintos parámetros (e.g., abundancia y longitud), los pesos se calcularon utilizando el coeficiente de variación en vez de la desviación estándar de los parámetros. El algoritmo se ejecutó durante 6000 generaciones, considerando para el análisis las últimas 4000, en las cuales se observó una estabilización de los errores cuadráticos y la correlación entre los valores observados y predichos. Se utilizaron como criterio de decisión para seleccionar las posibles soluciones aquellos conjuntos de parámetros que produjeran biomasas y capturas mensuales simuladas con una correlación mayor a 0.85 y 0.9, respectivamente, con los datos observados. Adicionalmente, debido a que la captura no fue modelada explícitamente, es posible observar conjuntos de parámetros que simulen mortalidades por pesca muy altas en algunos años, indicando que casi toda la cohorte desapareció luego de la explotación. Un estudio independiente utilizando un modelo estadístico de captura estructurado por edad para el periodo 1962–2008 (Díaz et al. 2010) estimó que las mortalidades por pesca fueron siempre inferiores a F = 1.8 año^–1. Así, se incluyó un criterio adicional para seleccionar las posibles soluciones, eliminando las que produjeran mortalidades por pesca superiores a F = 3 año^–1, lo que equivale a 2.5% de sobrevivientes a partir de la edad de un año.

Análisis de las series de reclutamiento

Se realizó el test STARS para cambios de régimen (Rodionov 2004, Rodionov y Overland 2005) para determinar la presencia de estos cambios en las series de parámetros del reclutamiento (N₀, L₀ y S₀), utilizando una longitud del periodo de búsqueda de 10 años y un nivel de significancia del 5% (α = 0.05). Dentro de cada régimen se calculó el valor medio ponderado para dichos parámetros, disminuyendo el efecto de los valores atípicos al utilizar los pesos dados por la fórmula de Huber para cinco desviaciones estándar (los pesos disminuyen para valores alejados más de cinco desviaciones estándar de la media aritmética) (Rodionov 2004, Rodionov y Overland 2005). Se analizó también la existencia de diferencias estacionales en los parámetros del reclutamiento, utilizando como variables explicativas la estación (verano y primavera) y los regímenes identificados. Se calcularon los patrones estacionales de las series de tiempo mensuales, para el periodo total y por regímenes. Para los regímenes detectados, se analizaron las diferencias en las varianzas en los distintos regímenes utilizando el test no paramétrico de Fligner–Killeen implementado en R y las diferencias entre las medias por estación (otoño y primavera) dentro de cada régimen utilizando el test no paramétrico de Wilcoxon tambien implementado en R (R Development Core Team 2010). Finalmente, se calculó la abundancia mensual de pre–reclutas y reclutas para todo el periodo estudiado utilizando la ecuación 1, y los patrones estacionales (promedios por mes) de las series de prereclutas y reclutas para todo el periodo estudiado y según los regímenes identificados.

RESULTADOS

Se encontraron 26 soluciones que cumplen con los criterios exigidos para un candidato a solución, con los cuales se estimó la distribución de los parámetros estimados (N₀, L₀ y S₀ de cada reclutamiento) y se construyeron intervalos de confianza. La mediana de los parámetros estimados y los intervalos de confianza al 95% se muestran en la figura 1.

La aplicación del algoritmo STARS para detección de cambios de régimen sugiere la existencia de dos cambios de régimen, uno en 1971 y otro en 1991. Otro cambio de régimen es potencialmente detectado en 2009 pero el test no es concluyente por estar al extremo más reciente de la serie. Los cambios de régimen sólo se presentaron en la serie de abundancia de reclutas, no observándose diferencias en la longitud ni la dispersión del reclutamiento a edad 0. En la figura 2 se muestran la serie de abundancias de reclutas a la edad 0 y los promedios ponderados para cada régimen.

En base a estos resultados se consideraron para el análisis tres regímenes: 1961–1971, 1971–1991 y 1991–2009. Se encontraron diferencias significativas en la varianza de la abundancia del reclutamiento entre los regímenes, con diferencias significativas entre las varianzas del régimen 2 y las de los demás (P < 10^–6), mas no entre los regímenes 1 y 3 (P = 0.2271). Considerando cada régimen por separado, el test de Fligner–Killeen no encontró diferencias significativas entre estaciones (otoño o primavera). En la tabla 1 se presentan los resultados para las diferencias entre las medias por estación (otoño y primavera) dentro de cada régimen. En la figura 3 se comparan las distribuciones de las abundancias del reclutamiento a la edad 0 entre regímenes y estaciones. A partir de la comparación de medias entre estaciones en el primer régimen (P = 0.0753) y la asimetría opuesta de la distribución del reclutamiento para ambas estaciones (fig. 3), y considerando además la alta variabilidad del reclutamiento, se considera que existen diferencias significativas entre estaciones (otoño y primavera) sólo en el primer régimen.

Las series de tiempo de prereclutas y reclutas se muestran en la figura 4, donde se pueden identificar los tres periodos cualitativamente distintos, con abundancias altas, bajas e intermedias, respectivamente. Las figuras 5 y 6 muestran diagramas de cajas de los patrones estacionales por mes de las series de prereclutas y reclutas, respectivamente, para todo el periodo estudiado y para cada uno de los regímenes detectados, en los que se nota que en la década de los 60 las diferencias estacionales fueron más marcadas que en décadas posteriores. Además de los cambios en la abundancia entre los diferentes regímenes, también se observan cambios en la estacionalidad de los patrones de reclutamiento y prereclutamiento.

DISCUSIÓN

Los parámetros N₀, L₀ y S₀ permiten modelar la variabilidad temporal observada en el reclutamiento manteniendo una estructura discreta del modelo. El parámetro S₀, la dispersión en las tallas de la cohorte, está directamente relacionado con la duración de los máximos de desove y permite representar de manera más apropiada el comportamiento de un desovador parcial como la anchoveta peruana. El parámetro L₀ representa la talla media de la cohorte en el momento teórico de desove. Bajo esta interpretación, los valores positivos y negativos de L₀ corresponden a una cohorte que tuvo un máximo de desove anterior o posterior a la fecha teórica supuesta, respectivamente. Conocidas esta longitud L₀ y las fechas teóricas de desove (fijadas arbitrariamente), es posible calcular el punto máximo real del desove a partir de la relación de crecimiento utilizada. En conjunto con S0, esto hace posible estimar la distribución de las fechas de nacimiento de cada una de las cohortes. Así, a partir de un modelo de dinámica poblacional que usa información de la pesquería y estimaciones acústicas de biomasa, es posible predecir atributos de la población que pueden ser contrastados con información reproductiva, la cual normalmente no se incluye de manera directa en los modelos de evaluación de stocks. Más aún, la forma de parametrizar el reclutamiento utilizada posibilita que la información reproductiva (producción de huevos, fracción desovante) pueda ser usada directamente en la calibración de modelos de dinámica poblacional, lo cual será explorado en futuros estudios. Esto es especialmente importante para la evaluación de un stock como el de la anchoveta peruana, en el que el efecto del ambiente no es sólo sobre el reclutamiento a nivel de su abundancia, sino que también puede influir en los patrones de desove y la duración de éste. Asimismo, parametrizar al reclutamiento en función de N₀, L₀ y S₀ nos permite darle una resolución de un mes a las estimaciones de prereclutas y reclutas. Una limitación del presente trabajo ha sido utilizar parámetros biológicos constantes (M, k, L_∞) para todo el periodo estudiado. Sin embargo, la variabilidad temporal de estos parámetros puede ser fácilmente incluida en el modelo y permitiría reproducir más adecuadamente la variabilidad del reclutamiento. Otra limitación ha sido usar la misma relación de crecimiento para los estadios larvarios, ya que el modelo considera tallas más pequeñas que las usadas para deducir las relaciones de crecimiento utilizadas. Ecuaciones de crecimiento obtenidas experimentalmente para las larvas de anchoveta podrían ayudar a reducir este sesgo.

El número de parámetros incluidos en el modelo se encuentra dentro del promedio utilizado en los modelos de dinámica poblacional, que pueden llegar a incluir cientos de parámetros (Quinn 2003). Una diferencia importante está en la elección de parametrizar el reclutamiento de la manera propuesta en vez del efecto de la pesca (selectividad variable en el tiempo, mortalidad por pesca). Por otro lado, las cotas impuestas a los parámetros disminuyen el riesgo de sobre–parametrización, al limitar las posibilidades de cada parámetro a un rango razonable según la información biológica de la especie. En la figura 1 se observa que las series de parámetros estimados no se acercan, en general, a los límites impuestos, lo que sugiere que la elección de dichas cotas fue eficiente. Sin embargo, definir cotas variables en el tiempo que sean estimadas a partir de la información directamente observada (duración de los máximos de desove para S₀, máximos de desove para L₀) sería una forma de mejorar las estimaciones de los parámetros de reclutamiento.

En los patrones estacionales del primer régimen (1961–1971) se observó un ingreso constante de reclutas desde octubre, con un promedio sostenido en 0.4 × 10¹² individuos, que alcanzó su punto máximo entre los meses de marzo y abril. Estos resultados son similares a los descritos por Jordán (1976), quien manifestó que el "patrón general de reclutamiento" se observa a partir de noviembre a diciembre y alcanza sus valores más altos de febrero a marzo como consecuencia de los desoves de invierno–primavera y un reclutamiento secundario en invierno proveniente de los desoves de verano. La regularidad del patrón estacional en la década de los 60 se puede atribuir a la marcada relación entre la densidad de la población adulta y el reclutamiento subsiguiente documentada por Csirke (1979). Esta regularidad se vio alterada durante 1972–1973 debido a la drástica disminución de la población de anchoveta debido al incremento de la presión pesquera que condujo al colapso del reclutamiento a comienzos de 1971, y que fue intensificada por la presencia de El Niño 1972–1973 (Pauly y Palomares, 1989). Durante el segundo régimen (1971–1991) la anchoveta peruana se vio fuertemente afectada por continuos eventos cálidos de El Niño durante los años 1976, 1982–1983 y 1987–1988 (Lagos 1984, Zuta et al. 1984, Macharé y Ortlieb 1993, Morón 2000, Trasmonte y Silva 2008) que, junto a los bajos niveles de abundancia del stock adulto no posibilitaron el éxito reproductivo y, por ende, el éxito del reclutamiento (Jordán 1976), observándose un patrón estacional sin máximos bien definidos y con reclutamientos medios similares durante la mayor parte del año en los resultados.

La recuperación de la anchoveta después de este periodo cálido fue lenta debido a la presencia de sardina que se constituyó en la especie dominante del sistema (Ñiquen et al. 1988). Esta recuperación marca el inicio de un nuevo régimen para la anchoveta, el tercero observado (desde 1991) en los resultados obtenidos en este trabajo y en el que la abundancia del reclutamiento aumenta con respecto al régimen anterior, pero sin alcanzar los niveles de la década de los 60. La estrategia reproductiva durante este último régimen mostró alternancias en los desoves, con periodos marcados de uno o dos puntos máximos de desove (Mori et al. este número, IMARPE 2010), que también son consistentes con los resultados en este trabajo donde los patrones estacionales del reclutamiento muestran dos ligeros máximos en noviembre y mayo (fig. 6), pero sin que ninguno predomine sobre el otro. Esta modificación en el proceso reproductivo de la anchoveta y la presencia de un ambiente favorable para su desarrollo (ENFEN 2009) marcan el ingreso de un mayor contingente de individuos a la pesquería en forma continua, con mayores valores de octubre a febrero y con un ligero repunte entre abril y mayo como se puede observar también en la figura 6.

Esta discusión nos motiva a pensar que los primeros resultados obtenidos son razonables a pesar de las limitaciones antes descritas. Es importante destacar que los patrones de variación mensual de prereclutas y reclutas no podrían ser descritos mediante modelos clásicos estructurados por edades. Adicionalmente, la reconstrucción de la estructura de la población con una resolución mensual abre muchas posibilidades para un estudio más detallado del proceso de reclutamiento y los factores que influyen en el mismo. Por otro lado, al parametrizar el reclutamiento dos veces al año, el modelo tiene un horizonte predictivo de hasta 6 meses sin necesidad de acoplarlo con modelos stock–reclutamiento. Esto hace que el modelo sea apropiado para la predicción de la estructura de la población sin ninguna suposición paramétrica sobre el proceso de reclutamiento, y para la reconstrucción de la evolución de los individuos juveniles que no son bien representados en las evaluaciones directas. Por lo tanto, este modelo se presenta como una alternativa interesante para el manejo de la anchoveta, complementaria a la evaluación directa que se realiza regularmente durante los periodos reproductivos más importantes (verano y primavera).

CONCLUSIONES

La forma en que el reclutamiento es parametrizado en el modelo presentado permite representar más adecuadamente la continuidad temporal del proceso de reclutamiento a pesar de modelar la dinámica poblacional en tiempo discreto. La parametrización usada, al modelar la variabilidad del reclutamiento no solamente en términos de su abundancia, sino de su longitud y dispersión, abre la posibilidad de relacionar directamente los parámetros del modelo con información reproductiva.

Los cambios de régimen y los patrones estacionales del reclutamiento obtenidos por el modelo son coherentes con los resultados descritos por estudios biológicos basados en observaciones directas.

El modelo permite la reconstrucción mensual de la estructura de la población, lo cual abre muchas posibilidades para un estudio más detallado del proceso de reclutamiento y los factores que influyen en el mismo; asimismo, es apropiado para la predicción de la estructura de la población, en especial de los individuos juveniles que no son bien representados en las evaluaciones directas.

Al comité organizador del I Simposio Iberoamericano de Ecología Reproductiva, Reclutamiento y Pesquerías por el financiamiento para la presentación de este trabajo en el citado evento. Se agradece al personal del Instituto del Mar del Perú (IMARPE) el apoyo brindado durante la realización de este trabajo. Agradecemos también a dos revisores anónimos sus valiosos comentarios. Este trabajo es una contribución al LMI DISCOH (Dinámicas del Sistema de la Corriente de Humboldt) y al proyecto ANR–VULNS–007 PEPS (Peru Ecosystem Projection Scenarios), en el marco de cooperación entre el IMARPE y el Institut de Recherche pour le Développement (IRD) de Francia.

REFERENCIAS

Agüero M. 1987. A biooeconomic model of the Peruvian pelagic fishery. In: Pauly D, Tsukayama I (eds.), The Peruvian Anchoveta and its Upwelling Ecosystem: Three Decades of Change. ICLARM Studies and Reviews 15, pp. 307–324.         [ Links ]

Alheit J, Ñiquen M. 2004. Regime shifts in the Humboldt Current ecosystem. Prog. Oceanogr. 60: 201–222.         [ Links ]

Andrzej O, Stanislaw K. 2006. Evolutionary algorithms for global optimization. In: Pintér J (ed.), Global Optimization: Scientific and Engineering Case Studies. Springer, pp. 267–300.         [ Links ]

Aranda M. 2009. Developments on fisheries management in Peru: The new individual vessel quota system for the anchoveta fishery. Fish. Res. 96: 308–312.         [ Links ]

Ayon P, Swartzman G, Espinoza P, Bertran A. 2011 Long–term changes in zooplankton size distribution in the Peruvian Humbodt Current System: conditions favouring sardine or anchovy. Mar. Ecol. Progr. Ser. 422: 211–222.         [ Links ]

Bakun A, Broad K. 2001. Climate and Fisheries: Interacting paradigms, scales, and policy approaches. The IRI–IPRC Pacific Climate–Fisheries Workshop, Honolulu, Hawaii, 14–17 November 2001.         [ Links ]

Beaugrand G, Reid PC, Ibáñez F, Lindley JA, Edwards M. 2002. Reorganization of North Atlantic marine copepod biodiversity and climate. Science 296: 1692–1694.         [ Links ]

Bertrand A, Ballón M, Chaigneau A. 2010. Acoustic observation of living organisms reveals the oxygen minimum zone. PLoS ONE 5: e10330; doi:10.1371/journal.pone.0010330.         [ Links ]

Cahuín S, Cubillos L, Ñiquen M, Escríbano R. 2009. Climatic regimes and the recruitment rate of anchoveta, Engraulis ringens, off Peru. Estuar. Coast. Shelf Sci. 84: 591–597.         [ Links ]

Csirke J. 1979. El reclutamiento en la anchoveta peruana y su dependencia de la población adulta en los años 1961–1976. Bol. Inst. Mar Perú 4(1): 29 pp.         [ Links ]

Csirke J. 1989. Introducción a la dinámica de poblaciones de peces. FAO Doc. Téc. Pesca 192: 82 pp.         [ Links ]

Díaz E, García C, Espinoza D, Guevara–Carrasco R, Csirke J, Ñiquen M, Vargas N, Argüelles J. 2010. Evaluación del stock norte–centro de la anchoveta peruana (Engraulis ringens) por un modelo estadístico estructurado por edades. Bol. Inst. Mar Perú 25: 57–62.         [ Links ]

Duboz R, Versmisse D, Travers M, Ramat E, Shin Y. 2010. Application of an evolutionary algorithm to the inverse parameter estimation of an individual based model. Ecol. Model., doi:10.1016/j.ecolmodel.2009.11.023.         [ Links ]

ENFEN. 2009. Comité multisectorial encargado del estudio nacional del fenómeno El Niño. Comunicado oficial No. 09/ 2009, http://www.imarpe.pe.         [ Links ]

Gulland J. 1968. Informe sobre la dinámica de la población de la anchoveta peruana. Bol. Inst. Mar Perú 1(6): 305–346.         [ Links ]

Gutiérrez M. 2000. Estimados de biomasa hidroacústica de los cuatro principales recursos pelágicos en el mar peruano durante 1983–2000. Bol. Inst. Mar Perú 19 (1–2): 136–156.         [ Links ]

IMARPE. 1965. Efectos de la pesca en el stock de anchoveta. Informe del Instituto del Mar de Perú 7, 16 pp.         [ Links ]

IMARPE. 1974. Informe de la cuarta sesión del panel de expertos de la evaluación del stock de anchoveta peruana. Bol. Inst. Mar Perú 2(10): 605–719.         [ Links ]

IMARPE. 2010. Informe sobre la quinta reunión de expertos en dinámica de población de la anchoveta peruana. Bol. Inst. Mar Perú 25 (1–2): 1–84.         [ Links ]

Jordán R. 1976. Biología de la anchoveta. Parte I: Resumen del conocimiento actual. FAO Tech. Rep. 185: 359–399.         [ Links ]

Lagos P. 1984. El Niño 1982–1983 y las alteraciones climáticas en el Perú. Rev. Com. Perm. Pacífico Sur 15: 95–119.         [ Links ]

Macharé J, Ortlieb L. 1993. Registro del fenómeno El Niño en el Perú. Bull. Inst. Francais d'Etudes Andines (Lima, Peru) 22(1): 35–52.         [ Links ]

Mori J, Buitrón B, Perea Á, Peña C, Espinoza C. 2009. Variabilidad interanual en la estrategia reproductiva de la anchoveta peruana en la región norte–centro del litoral del Perú. Cienc. Mar. 37: 513–525.         [ Links ]

Morón O. 2000. Características del ambiente marino frente a la costa peruana. Bol. Inst. Mar Perú 19(1–2): 179–204.         [ Links ]

Mullon C, Fréon P, Parada C, van der Lingen CD, Huggett JA. 2003. From particles to individuals: Modelling the early stages of anchovy in the southern Benguela. Fish. Oceanogr. 12: 396–406.         [ Links ]

Ñiquen M, Bouchon M. 2004. Impact of El Niño event on pelagic fisheries in Peruvian waters. Deep–Sea Res. II 51: 563–574.         [ Links ]

Ñiquen M, Freón P. 2006. A new record set by the Peruvian fishery. Globec International Newsletter, October 2006.         [ Links ]

Ñiquen M, Echevarría A, Mariategui L. 1988. Variaciones de desembarques de anchoveta (Engraulis ringens) y sardina (Sardinops sagax sagax) en el Perú en enero–abril de 1987. In: Salzwedel H, Landa A (eds.), Recursos y Dinámica del Ecosistema de Afloramiento Peruano. Bol. Inst. Mar Perú (special issue), pp. 265–271.         [ Links ]

Ñiquen M, Bouchon M, Cahuin S, Valdez J. 1999. Efectos de El Niño 1997–1998 sobre los principales recursos pelágicos de la costa peruana. Rev. Peruana Biol. (special issue): 85–96.         [ Links ]

Pauly D, Tsukayama I (eds.). 1987. The Peruvian Anchoveta and its Upwelling Ecosystem: Three Decades of Change. ICLARM Studies and Review 15, 351 pp.         [ Links ]

Pauly D, Palomares ML. 1989. New estimates of monthly biomass, recruitment and related statistics of anchoveta (Engraulis ringens) off Peru (4–14°S) 1953–1985. In: Pauly D, Muck P, Mendo J, Tsukayama I (eds.), The Peruvian Upwelling Ecosystem: Dynamics and Interactions. ICLARM Conference Proceedings 18, pp. 189–206.         [ Links ]

Perea A, Peña C, Oliveros–Ramos R, Buitrón B, Mori J. 2011. Potential egg production, recruitment, and closed fishing season of the Peruvian anchovy (Engraulis ringens): Implication for fisheries management. Cienc. Mar. 37: 585–601.         [ Links ]

Quinn II TJ. 2003. Ruminations on the development and future of population dynamics models in fisheries. Nat. Resour. Model. 16: 341–392.         [ Links ]

R Development Core Team. 2010. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3–900051–07–0, URL http://www.R–project.org.         [ Links ]

Rodionov SN. 2004. A sequential algorithm for testing climate regime shifts. Geophys. Res. Lett. 31: L09204; doi:10.1029/2004GL019448.         [ Links ]

Rodionov SN, Overland JE. 2005. Application of a sequential regime shift detection method to the Bering Sea ecosystem. ICES J. Mar. Sci. 62: 328–332.         [ Links ]

Saetersdal G, Valdivia E. 1964. Un estudio del crecimiento, tamaño y reclutamiento de la anchoveta (Engraulis ringens) basado en datos de frecuencia de longitud. Bol. Inst. Invest. Recurs. Mar. (Callao, Peru) I(4): 85–136.         [ Links ]

Schnute J, Fournier D. 1980. A new approach to length–frequency analysis: Growth structure. Can. J. Fish. Aquat. Sci. 37: 1337–1351.         [ Links ]

Schnute J, Richards L. 2002. Surplus production models. In: Hart P, Reynolds J (eds.), Handbook of Fish Biology and Fisheries. Vol. 2: Fisheries. Blackwell Publishing, 426 pp.         [ Links ]

Shin YJ, Rochet MJ, Jennings S, Field J, Gislason H. 2005. Using size–based indicators to evaluate the ecosystem effects of fishing. ICES J. Mar. Sci. 62: 394–396.         [ Links ]

Trasmonte G, Silva Y. 2008. Evento La Niña: Propuesta de definición y clasificación según las anomalías de la temperatura de la superficie del mar en el area Niño 1+2. Inf. Inst. Mar Perú 35(3): 199–208.         [ Links ]

Zuta S, Farfán M, Morón O. 1984. Características principales del mar frente al Perú durante la evolución del evento El Niño 1982–1983. Rev. Com. Perm. Pacífico Sur 15: 141–178.         [ Links ]

NOTA

English translation by Christine Harris.