Comparison of demographic and direct methods to calculate probabilistic maturation reaction norms for Flemish Cap cod (Gadus morhua).

Age and length at maturation have declined in many fish populations and this has been hypothesized to be a genetic change caused by high fishing mortality. Probabilistic Maturation Reaction Norms (PMRNs) have been used as a tool to gain a better understanding of the possible genetic nature of these changes. The demographic and direct methods are two ways to calculate PMRNs. The data requirements are more often met for the demographic method than for the direct method which requires the identification of recruit spawners. However, the demographic method relies on more assumptions than the direct method, typically assuming equality of growth and mortality rates for immature and mature individuals within an age class. This study provides the first direct comparison of demographic and direct methods and shows that both methods produce comparable results. Differences between methods are hypothesized to be owed to possible differences in growth rate between mature and immature individuals in Flemish Cap cod.

Introduction

The maturation process is the most important life-history transition in plants and animals which comes with large changes in energy allocation (Dieckmann and Heino 2007; Roff 1992, Stearns 1992; Bernardo 1993). While a fish is immature, energy is allocated to increase survival and growth. After maturation, energy is required for reproduction and this may decrease survival, deplete body reserves and reduce growth (Roff 1992; Berner and Blanckenhorn 2007; Saborido-Rey and Kjesbu, unpublished manuscript). For early maturing individuals, the decrease in lifetime reproductive potential and in lifespan can be even more marked than in late maturing ones. Indeed, as fecundity, egg size and egg quality can increase with fish size and age (Trippel 1995; Marshall and Frank 1999), females maturing at smaller size and younger age suffer from reduced reproductive success (Solemdal 1997; Trippel et al. 1997). At a population level, the decrease in age and length at maturation of individuals can produce a decrease in reproductive potential, (Domínguez-Petit 2007), and make the recruitment process more sensitive to environmental variability (Saborido-Rey and Kjesbu, unpublished manuscript).

In recent decades several studies describing a reduction in age and length at maturation in exploited species have appeared, including Pacific salmon Oncorhynchus spp. (Ricker 1981), European sole Solea solea (De Veen 1976), European plaice Pleuronectes platessa (Rijnsdorp 1993b), several populations of Atlantic Cod Gadus morhua (Jørgensen 1990; Trippel et al. 1997; Saborido-Rey and Junquera 1998), haddock Melanogramnus aeglefinus (Taylor and Stefánsson 1999) and Atlantic halibut Hippoglossus hippoglossus (Haug and Tjemsland 1986). In many exploited stocks this change in age and size at maturation has been attributed to high fishing pressure (Trippel 1995; Saborido-Rey and Junquera 1998; Law 2000; Grift et al. 2003; Olsen et al. 2004). The maturation process is determined by an ecological component that operates through phenotypic plasticity and an evolutionary component that works through genetic adaptation (Dieckmann and Heino 2007). Through their size selective action, fisheries can produce differential selection favouring genotypes coding for high fecundity and early maturation (Law and Grey 1989; Rijnsdorp 1993b). In addition, through the changes that fisheries produce in the ecosystem, inter and intraspecific competition can be altered, leading to an enhanced availability of food and producing earlier maturing individuals through phenotypic plasticity (Rijnsdorp 1993a; Taylor and Stefánsson 1999; Law 2000). Thus, fisheries can trigger changes in the maturation process through two pathways, through direct environmental effects manifested as phenotypic plasticity and through genetic selection.

Determining if changes in age and size at maturation have solely a phenotypic basis, have only a genetic basis or are actually a mixed effect of both phenomena is not a trivial issue. Reversal in maturation trends is expected to be much slower if genetic change is the major cause of a decrease in age and length at maturation (Law and Grey 1989; Law 2000; Olsen et al. 2005; Enberg et al. 2009). Therefore, knowing the origin of changes in reproductive parameters is very important and should be considered in management decisions. In the last two decades, four different methods have been used to separate phenotypic and genetic causes of shifts in maturation: experimental manipulation, comparative studies, countergradient or countertrend variation, and reaction norm analysis (Trippel 1995; Olsen et al. 2005; Dieckmann and Heino 2007; Sharpe and Hendry 2009; Conover and Baumann 2009).

The concept of norms of reaction was developed a century ago by Woltereck (1909). Stearns and coworkers (Stearns and Crandall 1984; Stearns and Koella 1986) developed the first bivariate reaction norm to examine the effect of growth rate on age and size at maturation. This methodology assumes that maturation is a deterministic process. However, since the microenvironment experienced for each individual can be different and because age and size are only proximate measures of maturation tendency, maturation should be considered as a variable process (Heino et al. 2002; Dieckmann and Heino 2007). To account for this individual variation, probabilistic maturation reaction norms (PMRNs) were developed (Heino et al. 2002). The first approach to this was called the direct method to calculate PMRNs (DirPMRNs) (Heino and Dieckmann 2008). DirPMRNs estimate the probability of becoming mature using the age and size-specific densities of immature and newly maturing individuals (hereafter called recruit spawners). Heino et al. (2002) applied this method to northeast Arctic cod. A variant of the direct method was applied for PMRNs estimation on herring Clupea harengus (Engelhard and Heino 2004) and smallmouth bass Micropterus dolomieu (Dunlop et al. 2005). In these works, determination of size and age when individuals mature was derived by means of the study of the distance between two consecutive annual rings in the fish otolith or scale. Morita and Fukuwaka (2006), applied the direct method to chum salmon (Oncorhynchus keta), for which maturing individuals were determined by return of salmon from the sea to their spawning territories in the rivers. However, data separating recruit spawner individuals from those that matured in earlier seasons are rarely available (Heino and Dieckmann 2008). Barot et al. (2004) overcame this difficulty by developing the demographic method (DemPMRN hereafter) (Heino and Dieckmann 2008). To deal with the lack of recruit spawner individual information DemPMRN usually makes two important assumptions: (i) mature and immature individuals have the same survival rates and (ii) all individuals within an age class have identical growth increments in the year considered (Barot et al. 2004). With these assumptions and auxiliary information on growth trajectories and maturation ogives it is possible to reconstruct the recruit spawner size and probability distribution by age class. The DemPMRN method permitted extension of the use of PMRNs to various populations of Atlantic cod Gadus morhua (Barot et al. 2004; Olsen et al. 2004, 2005), European plaice Pleuronectes platessa (Grift et al. 2003), American plaice Hippoglossoides platessoides (Barot et al. 2005), haddock Melanogrammus aeglefinus (Wright 2005) and sole Solea solea (Mollet et al. 2007).

The purpose of this study was the estimation and comparison of PMRNs obtained using both direct and demographic methodologies. Flemish Cap cod present a unique opportunity to compare the direct and demographic methods which, have never been directly compared before. The Atlantic cod population on the Flemish Cap (NAFO Division 3M) has supported an important commercial fishery for many years (Fernández et al. 2008) and age at 50% maturity (A50) has shown a dramatic decrease from 5 years in 1992 to 2 years old at the beginning of twenty first century (Saborido-Rey and Junquera 1998, 1999). Perhaps more importantly, gonads are examined histologically and recruit spawners and repeat spawners can be separated. This permits the calculation of PMRNs by the direct method as well as by the demographic method. Because direct method is applied to more informative data than demographic method (Heino and Dieckmann 2008), obtaining similar results for both techniques would reinforce the validity of using the demographic method for estimation of PMRNs when there's no possibility of separate recruit and repeat spawner individuals.

Material and methods

A total of 3058 cod ovaries were collected as a part of the European Union survey of the Flemish Cap in June-July from 1990 until 2006, comprising cohorts from 1983 to 2005. Fishing sets were conducted according to a stratified random design ranging in depth from 126 to 749 m, with a haul duration of 30 min (Vázquez et al. 1998). Total length of fish was measured and the otoliths were collected for age determination in the laboratory. Gonads were extracted, fixed, histologically processed and microscopically analyzed following the details provided in Saborido-Rey and Junquera (1998). The survey on Flemish Cap occurs after spawning is completed so that fish in the population are either fish that are immature, fish that are maturing to spawn the next year for the first time, or fish that have spawned already in the current year. These different maturity stages were identified following Saborido-Rey and Junquera (1998): immature females were identified when all the oocytes were in the primary growth stage and no evidence of previous spawning activity was found; recruit spawner females had oocytes in the cortical alveoli stage (CA) or vitellogenic stage (VO) but not postovulatory follicle (POF) nor atretic oocyte stage (AO); repeat spawner females were identified by the presence of CA, VO, POF and/or AO.

Different groupings of the available cohorts were formed in order to create more possibilities for calculating and comparing DemPMRNs and DirPMRNs (Table 1). Cohorts ranging from 1983 to 2005 were combined into six cohort groups. Other groupings combined these cohorts into three, two and one cohort groups.

Table 1

Number of immature, previous spawner and recruit spawner individuals by cohort group

	Immature	Recruit spawners	Previous spawners
Six groups
1983–1987	112	94	280
1988–1990	112	108	164
1991–1993	188	71	397
1994–1996	93	29	104
1997–2001	181	70	390
2002–2005	365	88	154
Three groups
1983–1989	139	132	341
1990–1996	366	170	604
1997–2005	546	158	544
Two groups
1983–1995	502	302	933
1996–2005	549	158	556
One group
1983–2005	1051	460	1489

For maturity ogive estimation, a logistic model was fit to the whole data set. Maturity state was introduced as a binary dependent variable, age and length were considered as continuous explanatory variables and cohort group (Cg) as a factor. Repeat and recruit spawner females were considered as mature individuals.:

where A is fish age, L the fish length and Cg, the cohort group.

With the parameters obtained from equation (1) and with length increments by age, the probability of maturing for each age group, and length was estimated as:

where p(a,l) is the probability of maturing at age a and length l, a−1 is the age previous to a and l–δl is the length at that previous age, δl being the length increment between age a–1 and age a, and o is the fraction mature. These calculations were done for each cohort group.

For the direct method calculations, previous spawners were not included. A logistic model was fit to maturing (recruit spawners) and immature individuals [equation (3)], with maturity state as a binary dependent variable. Age and length were included in the model as continuous variables and the cohort group as a factor:

The meaning of symbols in equation (3) is the same as in equation (1) with the difference that in equation (3) the logit is on the proportion maturing instead of the proportion of mature individuals.

For both the demographic and direct methods, a stepwise selection process was conducted to determine the best model fit for each of these cohort groupings. Variables which when added to the model didn't result in a significant decrease in the residual deviance or when included gave a higher Akaike Information Criteria were rejected from the model.

For the demographic PMRN calculation method, it was necessary to estimate the mean length increment between ages and the maturity ogive to calculate the probability of becoming mature (Barot et al. 2004). Yearly length increments were obtained by calculating the difference in mean size at age between consecutive years:

where δl is the mean length increment resulting from substracting the mean length at age a−1 from the mean length at age a. This process was carried out for every age in every cohort group.

Because of the high correlation observed between age and length in the data, the range of ages used to fit the models was shortened as much as possible. Only ages where both mature (recruit and previous spawners for DemPMRNs and only recruit spawners for DirPMRNs) and immature individuals occurred were considered. Finally for the demographic method 3000 individuals were used, ranging from ages 1 to 7. For the direct method, 1318 individuals were used, from ages 2 to 5.

For demographic methods, it has been established that a minimum number of 50 individuals is necessary by age and cohort for calculating PMRNs and to estimate confidence intervals by means the bootstrap technique (Barot et al. 2004). The same amount of individuals has been also considered in this work for direct method calculations. Therefore, reaction norm midpoints (length for which the probability of maturing is 50%, hereafter Lp50) were calculated for age-cohort group combinations that fulfilled this requirement. For these age-cohort combinations, estimated probabilities of maturing were generated through equation (2). With the demographic method, Lp50 (hereafter DemLp50) was calculated by interpolation as the length where the probability of maturing was 50%. Lp50 for the direct method (DirLp50) was calculated for every cohort group and age studied by solving the equation resulting from fitting of model (3) for the length where the probability of maturing was 50%.

The 95% confidence intervals for both DirLp50 and DemLp50 were calculated by using standard bootstrap techniques (Efron and Tibshirani 1993). A new dataset was created by randomly sampling original data, stratified by length for every year. This is consistent with the original sampling design in the Flemish Cap European Union Survey (Vázquez et al. 1998). With the resampled data, reaction norm midpoints (Lp50) were calculated by means of the process described above. The data set was resampled 999 times and the confidence limits of the reaction norm midpoints were estimated as the 2.5 and 97.5 percentiles of the distribution of the 1000 midpoints of each age and cohort (999 resamplings + 1 original).

Results

Logistic models for the DemPMRNs and the DirPMRNs methods were fit for every cohort grouping option (Tables 2 and 3). The final general logistic model selected as the best fit to the data changed depending on the number of cohort groupings considered. For the DemPMRN estimate the saturated model [equation(1)] was only the best fit when cohorts were divided into six cohort groups (Table 2). When individuals were classified into three cohort groups a logistic model with no three way interaction and no length:cohort group interaction was selected because adding these terms did not significantly improve the model fit. When cohorts were divided into two groups the three way interaction was not included in the best model. When grouping all cohorts into a unique cohort group, the model with the significantly lowest value of residual deviance included length and age and their interaction. For the DirPMRN calculations the complete model including all interactions was selected as the best fit to the data when cohorts were divided into six cohort groups (Table 3). When grouping in three and two cohort groups, no three way interactions were included in the final models. The model with length and age as main effects and their two way interaction had the best fit to the data when grouping cohorts in a single group, from 1983 to 2005.

Table 2

Results of the logistic regressions for the DemPMRN method. Probably of being mature was fit to age, length and cohort group. Models were chosen by stepwise regression

	d.f.	Explained deviance	Residual d.f.	Residual deviance	P(>|Chi|)
Six groups
NULL			2999	3886
Length	1	2020	2998	1866	<0.001
Cohort group	5	129	2993	1737	<0.001
Age	1	60.8	2992	1676	<0.001
Length:cohort group	5	138.2	2987	1538	<0.001
Length:age	1	17.7	2986	1520	<0.001
Cohort group:age	5	37.4	2981	1483	<0.001
Length:cohort group:age	5	15.2	2976	1467	<0.001
Three groups
NULL			2999	3886
Length	1	2020	2998	1866	<0.001
Cohort group	2	73.5	2996	1792	<0.001
Age	1	114.4	2995	1678	<0.001
Length:age	1	129.1	2994	1549	<0.001
Cohort group: age	2	53.1	2992	1496	<0.001
Two groups
NULL			2999	3886
Length	1	2020	2998	1866	<0.001
Cohort group	1	8.2	2997	1858	<0.001
Age	1	59.4	2996	1798	<0.001
Length:cohort group	1	50.8	2995	1747	<0.001
Length:age	1	104.4	2994	1643	<0.001
Cohort group: age	1	38.4	2993	1605	<0.001
One group
NULL			2999	3886
Length	1	2020	2998	1866	<0.001
Age	1	49.3	2997	1816	<0.001
Length:age	1	118.7	2996	1698	<0.001

Table 3

Results from Logistic Regresion for the DirPMRN method. Probably of maturing is fit to age, length and cohort group. Models were chosen by stepwise regression

	d.f.	Explained deviance	Residual d.f.	Residual deviance	P(>|Chi|)
Six groups
NULL			1317	1642
Length	1	421	1316	1221	<0.001
Cohort group	5	29.8	1311	1192	<0.001
Age	1	18.5	1310	1173	<0.001
Length: cohort group	5	28.5	1305	1145	<0.001
Length:age	1	11.9	1304	1133	<0.001
Cohort group:age	5	20.5	1299	1112	<0.001
Length: cohort group:age	5	27.1	1294	1085	<0.001
Three groups
NULL			1317	1642
Length	1	421	1316	1221	<0.001
Cohort group	2	16.8	1314	1205	<0.001
Age	1	37.1	1313	1167	<0.001
Length: cohort group	2	12.1	1311	1155	<0.001
Length:age	1	13.8	1310	1141	<0.001
Cohort group:age	2	20.4	1308	1121	<0.001
Two groups
NULL			1317	1642
Length	1	421	1316	1221	<0.001
Cohort group	1	0.5	1315	1221	<0.001
Age	1	11.6	1314	1209	<0.001
Length: cohort group	1	16.1	1313	1193	<0.001
Length:age	1	26.1	1312	1167	<0.001
Cohort group: age	1	25.9	1311	1141	<0.001
One group
NULL			1317	1642
Length	1	421	1316	1221	<0.001
Age	1	6.83	1315	1214	<0.01
Length:age	1	32.2	1314	1182	<0.001

Lp50 was calculated using the selected models for every grouping, but not for all ages from the initial range used to fit the models. Only those ages for which more than 50 individuals were available in a cohort group were used for calculating the Lp50. Ages that met this criterion are shown in Table 4. For such ages, DirLp50 was estimated using the logistic model fit for each grouping (Table 4). For the demographic method, except for those age-cohort group combinations marked in Table 4, Lp50 was estimated by interpolating the size at which the estimated probability of maturing was 0.5. For those ages marked in Table 4, no interpolation was possible since 0.5 was out of the range of estimated probabilities of maturing. But the maximum probability in this range was so close to 0.5 that extrapolation was considered a reliable way to estimate Lp50.

Table 4

Length (in cm) at which the probability of maturing is estimated to be 50% for every age studied (in years)

Cohort group	Age	N DirPMRN	DirLp 50 (cm)	N DemPMRN	DemLp 50 (cm)	Case
1983–1987	4	52	44.7	60	44.0	1
1983–1987	5	50	47.7	128	43.0	2
1988–1990	3	125	49.5	130	49.0	3
1988–1990	4	81	48.1	169	46.0	4
1991–1993	3	139	46.3	144	44.0	5
1994–1996	3	71	44.4	76	43.0	6
1997–2001	2*	194	47.2	194	48.6	7
2002–2005	2*	345	57.7	347	53.8	8
1983–1989	4	83	48.4	115	45.0	9
1983–1989	5	62	49.0	184	42.0	10
1990–1996	3	325	47.8	340	46.0	11
1990–1996	4*	99	33.6	402	24.6	12
1997–2005	2*	539	52.2	541	53.5	13
1983–1995	3	341	46.9	357	45.0	14
1983–1995	4	182	43.8	517	37.0	15
1996–2005	2*	542	52.1	544	52.2	16
1983–2005	2	646	47.2	648	46.4	17
1983–2005	3	418	46.2	556	43.9	18
1983–2005	4	190	44.1	762	39.3	19

Note: Only those ages and cohort groups for which number of individuals (N DirPMRN and N DemPMRN) was more than 50 were considered for calculations. The ‘Case’ column indicates the order in which the result is presented in Fig. 1.

Cases for which extrapolation was used to calculate DemLp50.

Figure 1

Estimated DemLp50 and DirLp50 with their corresponding 95% confidence intervals calculated by bootstrap resampling. The abscissa axis represents the case study corresponding to all the age and cohort groupings (see Table 4).

DirLp50 and DemLp50 were highly correlated with a Pearson Correlation Coefficient of 0.93. Both estimates of Lp50 showed the same pattern across the age-cohort group combinations (Fig. 1). Confidence intervals for Lp50 calculated for both methods overlapped substantially with the exception of age 4 for the cohort groups 1983–1995, 1990–1996 and 1984–2005 (cases 12, 15 and 19 respectively in Fig. 1). With the exception of age 2 in cohort group 1997–2001, 1997–2005 and 1996–2005, DemLp50 was always less than DirLp50 (Table 4). The mean difference between DirLp50 and DemLp50 estimation was 2.4 cm and the standard deviation of the differences was 2.0 cm.

Discussion

Despite the small difference in the Lp50 between the methods, both captured the process of maturation in the same way showing the same pattern. A high correlation between calculated values of DirLp50 and DemLp50 was obtained and with few exceptions the 95% confidence intervals overlapped in the cases studied. These results lead us to conclude that both the demographic and direct methods of calculating PMRNs are able to capture the process of maturation in a similar way.

The direct method for PMRNs estimation (DirPMRN) considers only immature and recruit spawner individuals for which age and size are known. In this way, with representative samples of the population, the probability of becoming mature for age and size can be easily and reliably calculated by a logistic regression fit to the proportions of mature individuals at age by size class (Heino et al. 2002). The demographic method (DemPMRN) uses less informative data (Heino and Dieckmann 2008), with no recruit spawners identification. The demographic method is based on comparing proportions of mature individuals at age and size at two consecutive time intervals, (two ‘snapshots’ of the same cohort) in order to calculate the expected probability of becoming mature for an age and size. Usual assumptions of this methodology are equal survival and growth rates for immature and mature fishes. Analyses carried out by Barot et al. (2004) examined the impact of deviations from these assumptions on results and showed that modest violations of these assumptions had no significant effects on the estimation of DemPMRN.

It is generally argued that when an individual begins the process of maturation changes in energy allocation are expected which are likely to decrease not only survival but also growth rates (Roff 1992; Stearns 1992, Saborido-Rey and Kjesbu 2007 in press and references therein). There has been no study of changes in mean length between consecutive ages in Flemish Cap cod that could detect a change in growth with maturation. Yoneda and Wright (2005) found that Atlantic cod maturing during their study were initially faster growing, but their allocation of energy into reproduction led to comparable growth rates to those that did not mature. If such changes in Atlantic cod growth after maturation can be extrapolated to Flemish Cap cod, violations of the assumptions of the demographic PMRN method may have occurred in the present study. The almost systematic lower value of DemLp50 could have been produced from such a violation of the assumptions of this method.

The DemPMRN method was developed because in most survey programmes, the maturity state of gonads is evaluated by macroscopic observation, and thus it is impossible to determine if a mature gonad comes from a recruit or a repeat spawner (Saborido-Rey and Junquera 1998). Heino et al. 2002 applied the direct method to northeast Arctic cod, for which the maturity stage of an individual's gonad (determined by macroscopic inspection) was combined with the study of patterns in its otoliths (spawning checks) to distinguish repeat and recruit spawners (Rollefsen, 1933). A variant of the direct method was applied to the chum salmon Oncorhynchus keta (Morita and Fukuwaka 2006), but in this case as individuals were caught on return to their river, determination of maturation status was a trivial issue. Other studies where DirPMRNs were calculated are those on Norwegian spring-spawning herring Clupea harengus (Engelhard and Heino 2004) and smallmouth bass Micropterus dolomieu (Dunlop et al. 2005). Both studies determined the maturation state indirectly, the former by means of patterns of growth on scales, and the latter through discriminant analysis on backcalculated individual growth trajectories. However, these studies employed important assumptions in fish-otolith and fish-scales growth relations and in the determination of recruit spawners using patterns in annual otoliths or scales rings (Francis 1990; Wright et al. 2001 and references therein). Moreover, macroscopic gonad staging entails possible confusion between resting and immature individuals (Ajiad et al. 1999) that might produce inaccurate maturation ogives. For Flemish Cap cod, histological maturity staging began in 1990. Histological analyses of gonads permits the determination of recruit spawners through the identification of different oocytes stages, disentangling recruit spawners from previous spawners and immature individuals with limited error (Zamarro et al. 1993; Saborido-Rey and Junquera 1998). The histology method for determining the maturity stage for Flemish Cap cod was introduced because the European Union Flemish Cap survey is conducted during June-July, when spawning has already concluded for the year and the macroscopic determination of the maturity stage of gonads is impossible. The duration of postovulatory follicles and atretic oocytes in cod ovaries following spawning permits the separation of previous spawners from recruit spawners, which can have only oocytes in vitellogenic and/or cortical alveoli state. Thus, the histological maturity staging avoids the confusion of resting mature individuals with immature ones and, at the same time, allows one to separate recruit spawner individuals from repeat spawners.

Because of the lack of recruit spawner data, studies where PMRNs have been estimated using both the DirPMRN and the DemPMRN method have never been conducted before. The availability of Flemish Cap cod recruit spawners, and at the same time, representative sampling of both immature and mature fish, has made it possible to calculate and compare reaction norms obtained with both DirPMRN and DemPMRN. The results obtained in this work provide a direct test of the two ways of calculating PMRNs. and support the use of the DemPMRN as a reliable way of estimating PMRNs for populations where recruit and repeat spawners can't be separated in the data. Use of the Direct method is recommended when appropriate data are available. Future work will analyze the evolution of PMRNs in Flemish Cap cod and its relation with environmental variables.