Qiuyun Ma, Yan Jiao, Can Zhou, Yiping Ren

aCollege of Fisheries, Ocean University of China, Qingdao, 266003, China

bDepartment of Fish and Wildlife Conservation, Virginia Polytechnic Institute and State University, Blacksburg, VA, 24061, USA

cLaboratory for Marine Fisheries Science and Food Production Processes, Qingdao National Laboratory for Marine Science and Technology, Qingdao, 266000, China

dCollege of Marine Sciences, Shanghai Ocean University, Shanghai, 201306, China

Keywords:

ABSTRACT Sexual and spatio-temporal variations have been observed in the life history parameters of many aquatic species and their causes have been related to harvesting pressure and environmental changes. This study aims to explore sexual, spatial and temporal variation in the growth and maturity through weight-at-length, length-at-age, and maturity-at-length relationships for Lake Erie Walleye (Sander vitreus) as a case to test some hypotheses. Hypotheses on whether harvest pressure and environmental changes (both local and global scale) caused the temporal changes of these life history traits were further diagnosed. Sexual and spatio-temporal variations in these life history traits were formulated using mixed-effects models. Our study found that geographic basin, sex,year and cohort all have substantial effects on the growth and maturity of Walleye based on survey data from 1989 to 2015. Multiple factors including water supply of Lake Erie, temperature, fishing pressure of Walleye, and global climate factors were found to correlate with the temporal variations of growth and maturity of Walleye significantly. Our findings should contribute to the future interpretation of Walleye life history variations and population dynamics. The methodology constructed in this study could be applied to explore the heterogeneity and impacting factors for other species in aquatic ecosystems.

1.Introduction

Individual growth and maturity traits influence our understanding of population dynamics and fisheries management decisions (Chen et al.,2011; Watson et al., 2012). Individual growth patterns most commonly described by the weight-at-length, length-at-age, and maturity is studied through maturity-at-length or at-length relationships. Weight-at-length relationships are used to describe the dimensionality of growth (Keys,1928). The von Bertalanffy growth model, analyzing length-at-age relationships, is widely used in fisheries to compare the growth patterns among populations or species (Ricker, 1975; Chen et al., 1991). Since the growth of individuals is one important component of population dynamics, the length-at-age growth model is widely used in many stock assessment models, such as the yield-per-recruitment model and age/length-structure models (Beverton & Holt, 1957; Doll et al., 2017;Jiao et al., 2005; Methot, 1989). Maturity-at-length relationships, which describe the probability of maturity at a certain age usually by a logistic curve (Nash et al., 2010; Quinn & Deriso, 1999), are important considerations in the research on reproduction, recruitment, and population dynamics, for example, a defined maturity-at-length relationship is required to estimate the spawning stock size in stock-recruitment or spawner-per-recruitment models (Brodziak et al., 2014; Sun et al.,2005).

Sexual and spatio-temporal variations are commonly observed and have received much attention in modelling studies (Gertseva et al.,2010; Jiao et al., 2010; Quist et al., 2002). Yellow Croaker (Larimichthys polyactis) along the Chinese coast in high latitudes tend to gain more weight than those at the same length in lower latitude (Liu et al., 2012;Ma et al., 2017). Sauger (Sander canadensis) and Walleye (Sander vitreus)in the Ottawa River, Canada, were found to have sex specific life history patterns with females longer at age and mature later than the males(Haxton, 2015). Using generalized linear mixed models for Albacore Tuna (Thunnus alalunga) across the South Pacific, Farley et al. (2014)found that the proportion of mature females at length varied significantly with latitude and time of year.

Fig. 1.Survey stations in different basins of Lake Erie Ontario partnership gill net survey (1989–2015).

To investigate the variations of growth and maturity, mixed effects models have become of high interest, because of their advantage to account for multiple sources of heterogeneity (e.g. spatial variations,sexual variations) in growth and maturity (Aggrey, 2009; Cnaan et al.,1997; Farley et al., 2014). This kind of approach provides a unifying framework for disparate statistical methods, which includes both fixed effect parameters and random effect parameters (Gelman & Hill, 2006;Helser & Lai, 2004; Jiao et al., 2010; Thorson & Minto, 2015). In the mixed effects model, the fixed effects relate to the entire population, and the random effects pertain to the value of the parameter fitted for individuals or subgroups with heterogeneity (Bernal-Rusiel et al., 2013;Cnaan et al., 1997; Thorson & Minte-Vera, 2016). Therefore, the random effects are estimated to describe the distribution of parameters among multiple sources of heterogeneity, which can avoid the problem of overparameterization (Aggrey, 2009; Jiao et al., 2010; Schinckel et al.,2005).

Driving or correlated factors that have been found to influence growth and maturity include fishing pressure, population density, local environmental factors (e.g., water temperature) and large-scale climate oscillations (Ficke et al., 2007; Kuparinen & Meril¨a, 2007). Busch et al.(1975) revealed that water temperature is an important factor determine the year-class strength for Walleye, since water warming could extend the incubation period, increasing the exposure of eggs to worse environment. Climate oscillations such as North Atlantic Oscillation(NAO) and El Ni˜no were revealed to potentially influence the freshwater ecosystems, through the influence on its physics, hydrology, and chemistry (Ficke et al., 2007; George et al., 2004; Straile et al., 2003). By reducing population abundance, fishing pressure can generate density-dependent compensatory responses in individual growth and maturity (Kuparinen & Meril¨a, 2007; Trippel, 1995).

Relationships between the potential factors and life history traits were usually explored separately by previous studies (Ficke et al., 2007;Kuparinen & Meril¨a, 2007). However, these factors (e.g. fishing pressure, local environment or global climate factors) tend more likely to influence the growth, maturity and mortality of species simultaneously and correlated in many aquatic ecosystems. For example, Lake Erie is a productive ecosystem but is heavily influenced by water policy (Clean Water Act of 1972 and Great Lakes Water Quality Agreement), invasive species, fishing behavior (including policy on fishing intensity) and climate changes (Bates et al. 2015a, b; Ludsin et al., 2001; Stepien et al.,2005; Xu et al., 2018). In this case, it’s difficult to figure out a couple of impacting factors to interpret the variations of life history traits for some important species.

Lake Erie Walleye supports one of the largest freshwater fisheries in North America with the commercial sector exclusive to Ontario and the recreational sector largely in four U.S. states. Additionally, as a primary predator, Walleye influences the lake ecosystem dynamics (Kershner et al., 1999; Knight & Vondracek, 1993). Its research and understanding,especially on its life history, population structure and dynamics, are extremely important for the optimization of its management. Lake Erie Walleye is managed through five management units (MU) with various size limit and recommended allowable harvest from statistical catch-at-age stock assessment models (Barton 2011; Kayle et al., 2015;WTG, 2016). The MU1, MU2 and MU3 for Lake Erie Walleye are corresponding to the west, west central and east central basins, while the south portion of MU4 is the Pennsylvania ridge, and the remaining portion of MU4 and MU5 in Ontario waters is the east basin of Lake Erie geographically. A recent tagging study in Lake Erie documented different movement rates and natural mortality rates among Walleye populations from different basins (Vandergoot & Brenden, 2014), suggesting that Walleyes from different basins may develop a unique set of life history traits, including growth and maturation (Muth & Wolfert,1986; Vandergoot & Brenden, 2014; Wolfert, 1969). As the foundation of the stock assessment, if the growth and maturation of Lake Erie Walley vary among different basins, different sex or other aspects, such variation needs to be considered in the population dynamics modelling and management decision making (Lorenzen, 2016).

This study aims to explore sexual, spatial and temporal variations in the growth and maturity of Lake Erie Walleye, and diagnose whether multiple factors (harvest pressure and environmental changes) caused the temporal changes in these life history parameters. First, mixed effects models were used to estimate growth and maturity variations among sex, basin, year and cohort (birth year). Second, the temporal variations (inter-annual and cohorts) of growth and maturity were analyzed to investigate their relationships with the fishing intensity and environmental factors in different spatial scales.

2.Materials and Methods

2.1.Data sources

We used fish age, length, body mass, and maturity data from the Lake Erie Ontario partnership gill net survey (Ontario Ministry of Natural Resources and Forestry, 2016). The survey covers 5 basins (west (W),west central (WC), east central (EC), Pennsylvania ridge (PR) and east basins (E)) from August–November and follows a stratified-random design (Fig. 1). A total of 49,490 Walleye individuals were sampled during the 27 years of the Lake Erie surveys (1989–2015) (Table 1 and Fig. 2). For each landed fish, its weight (g), total length (cm), sex,maturity, and age (in years) were recorded at the time of sampling(Ontario Ministry of Natural Resources and Forestry, 2016). The weight of Walleye collected ranged from 18 to 5431 g, with a mean of 749 g,while the total length varied from 12.9 cm to 78.7 cm with a mean of 39.5 cm (Fig. 3). The age of Walleye collected was from 0 to 24,frequently age 0 to 7 (Fig. 3).

Table 1 The sample size among different basins and sexes for Walleye Sander vitreus in Lake Erie from the Ontario partnership survey.

2.2.Models

A power law weight-at-length model (Keys, 1928) with lognormal error distribution was used equation (1).

whereWis the wet weight of an individual fish (g),Lis the total length(cm),ais the condition coefficient, andbis the allometric growth coefficient. Before parameter estimation, the weight-at-length model was log-transformed as

The von Bertalanffy equation (von Bertalanffy, 1938) with a lognormal error distribution was used to model Walleye growth

wheret(year) is the age of Walleye,L∞(cm) is the asymptotic length;t0is the age at which length is zero, andK(1/year) is the exponential rate of approach to the asymptotic length.

The maturity-at-length was modeled as a generalized linear model with a logit link function

Fig. 2.The sample size (# of fish sampled) among different years and cohorts for Walleye Sander vitreus collected from Lake Erie from Ontario partnership survey (1989–2015).

Fig. 3.The distribution of weight, length and age for Walleye Sander vitreus collected from the Lake Erie Ontario partnership gill net survey (1989–2015).

Table 2 Akaike Information Criterion (AIC) values and mean squared error (MSE) from cross-validation of alternative growth and maturity models (Scenario 1) for Lake Erie Walleye Sander vitreus

whereMtis the sexual maturity status of individuals at aget,smeasures the slope of the curve, andL50is the length at which 50% individuals are mature.

In the mixed effects models, sexual, spatial, inter-annual and cohort variation are considered as random variations of parameters (ln(a),K,and -s×L50(intercept)). Different combinations of variation are considered for all the weight-at-length, length-at-age and maturity-atlength relationships in this study. Among them, the models without variations considered from other factors are called global models (M1,equations (2)–(5)), in which only fixed effects are considered.

The mixed effects models (M2-M8, equations ) considered random effects from Sex (M2,RE=REs), Basin (M3,RE=REb), Sex +Basin (M4,RE=REs+REb), Year (M5,RE=REy), Cohort (M6,RE=REc), Sex +Basin +Year (M7,RE=REs+REb+REy), Sex +Basin +Cohort (M8,RE=REs+REb+REc), respectively (Table 2).

Because the cohort coverages in the earlier and latest years were very limited, for example only one specimen sampled was born in 1979 and all specimen born in 2015 were age 0, two alternative model scenarios were applied to estimate the random effects for cohorts while retaining the years with limited coverage of the cohort data. Scenario 1 modeled all the cohorts including those with low samples in the earlier and late years; scenario 2 used global tendency for cohorts born at the beginning and latest period of the sampling years, i.e., in 1979–1985, and 2012–2015, without considering their cohort random effects. Random effects estimate for cohorts born in 1986–2011 would be compared between these two scenarios. By comparing these 2 model scenarios we would be able to find out whether low sample size in the earlier and late years causes noise rather than true signals in the estimate, and whether the random effect during 1986–2011 is robust to the low sample years and cohorts. If no significant differences are found between them especially for the years of 1986–2011, random effects of cohorts born in 1986–2011 from scenarios 1 would be used in the following analysis.

Model comparison was based on Akaike Information Criterion (AIC)(Akaike, 1974; Burnham & Anderson, 2001) and a ten-fold cross--validation (Kohavi, 1995), in which the mean squared error (MSE) between the predicted and observed values for testing data was estimated.Lower AIC and MSE values indicate a better model. All the computations were conducted in R, using packageslme4andnlme(version: R x64 3.3.2) (Bates et al. 2015a, b, Team R. 2016; Pinheiro et al., 2017).

2.3.Driving or correlated factors considered

As one of the most important fisheries in Lake Erie, whether Walleyeexperience life history great changes caused by harvest, and whether other biotic and abiotic factors influence its life history, have been of high concern for Walleye fishery. Based on the estimated temporal random effects on growth and maturity from the selected models, which can best reveal the temporal random effects (years 1989–2015 and cohorts 1986–2011, respectively), we attempted to investigate the possible effects of biotic and abiotic factors in both the local (local biological and environmental factors in Lake Erie) and global (large climate ocean oscillation that often impact in a large spatial scale or even globally such as NAO etc.) scales (Table 3, Appendix Fig. 1).

Table 3 The loadings of original variables to the first 5 principal components (PC) for the biotic and abiotic factors.

Biotic factors contain the chlorophyll A concentration (the primary productivity of the ecosystem, available from Ocean Color Web(NASA-GSFC, 2017)), the abundance of Yellow Perch (Perca flavescensi,the prey and competitor of Walleye (Forney, 1974)) and Walleye (to test for density dependent growth and maturation), and fishing intensity(fishing mortality for all age 2+individuals, sport fishery effort and commercial fishery effort). The fishery data and abundance were extracted from the most recent stock assessment reports (WTG, 2016;YPTG, 2017).

Abiotic factors were the environmental factors in both local and larger scales. Hydrologic data of Lake Eire were obtained from the Great Lakes Environmental Research Laboratory, including water temperature, precipitation, lake evaporation, ice cover days, ice melt date, and the water supply (precipitation +runoff - evaporation) (NOAA GLERL 2017). The annual averages across the whole lake of the above factors in Lake Erie were used in this study. Large scale climate indices, i.e.Southern Oscillation Index (SOI), Pacific Decadal Oscillation (PDO),NAO and Pacific-North American Pattern (PNA), were obtained from National Oceanic and Atmospheric Administration (NASA-NOAA,2017). The annual averages of these indices were used and normalized in this study before the following analysis.

Fig. 4.Growth and maturity curves fitted from the global model M1 for Lake Erie Walleye Sander vitreus. Data were in dots; red curves are estimated functions. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5.The random effects for cohorts born in 1986–2011 estimated from two scenarios in the two mixed effects models (M6 and M8, considering random effects from cohorts and random effects from sexes, basins and cohorts, respectively). The red lines represent the line y =x. See subsection Models in the Materials and Methods section for details of these two scenarios. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6.Random effects on parameter ln(a) in weight-at-length relationships among different models of Lake Erie Walleye Sander vitreus

These multiple biotic and abiotic factors share substantial correlation with each other, and such correlations complicate the interpretation of the result if we calculate the relationship of each factor with the growth and maturity characteristics of Lake Erie Walleye. Therefore, we first used principal component analysis (PCA) to convert these correlated variables into uncorrelated principal components (PC), which was computed by the functionprcomppackagestatsin R (version: R x64 3.3.2) (Hotelling, 1933; Team, 2016). In the PCA, the loadings of these factors are their coefficients in the linear combination predicting the PCs.

The selected principal components, which explained most variance(≥80%) of these factors, were used as the predictors for the inter-annual or cohort parameter (lna,K, andL50, respectively) in a linear regression.Their correlations were assumed to be significant whenPvalues, derived from the linear regressions between the parameter and all five principal components, were less than 0.05. Incorporating their loadings to the principal components, the possible influence of these factors was evaluated on the growth and maturity of Lake Erie Walleye.

3.Results

Based on the pooled data, the global models (M1) without random effects were constructed for the weight-at-length, length-at-age and maturity-at-length relationships (Fig. 4). In the weight-at-length relationship of Lake Erie Walleye, parameterln(a)was estimated as −5.747(95% CI =−5.758, −5.736). The estimate of parameterbwas 3.295(95% CI =3.292, 3.298). In the von Bertalanffy growth model,L∞,Kandt0were estimated as 61.514 cm (95% CI =61.313, 61.717), 0.370 year−1(95% CI =0.367, 0.373) and −0.964 year (95% CI =−0.972,−0.956), respectively. The estimate of parameterL50in the maturity-atlength model was 41.1 (95% CI =39.4, 42.9), while thesestimate was 0.211 (95%CI =0.206, 0.215).

The random effects for cohorts born in 1986–2011 estimated from two model scenarios were the same (Fig. 5), and the results from pairedt-test revealed thatP-values were 1. Therefore, the mixed effects models considered random effects of cohorts born in 1986–2011 from scenarios 1. Both the AIC and MSE from cross-validation revealed that mixed effects models (M2-M8) were better than the global models (M1) for the growth and maturity of Walleye (Table 2). Mixed effects models with multiple random effects (i.e. M4 with sex and basin; M7 with sex, basin,and year; M8 with sex, basin and cohort) performed better than those models with single random effect (i.e. M2 with sex, M3 with basin, M5 with year and M6 with cohort). Residuals from three global models(Appendix Figs. 2–4) and all the other mixed effects models, revealed evident sexual, spatio-temporal variation of Walleye growth and maturity, and matched the estimated random effect in parameters (see below). The models with random effects of year (M5) fitted better (with lower AIC and MSE values) than the models with random effects of the cohort (M6) for the weight-at-length and maturity-at-length relationships, but the AIC and MSE values of M5 is a little higher than that of M6 for the length-at-age relationships. For weight-at-length and maturityat-length relationships, M7 with random effects of sexes, basins and years performed best, while M8 with random effects from sexes, basins and cohorts fitted best for length-at-age relationships.

Female Walleye are longer in length than the male individuals at the same age, and mature at much longer size. The males gained a little more (about 0.035%) weight than the females at the same length(Fig. 6). The sexual variation in the length-at-age and maturity-at-length relationship were considerable (Figs. 7 and 8). The female is about 9.51% (average among different ages) longer than the male at the same age. At length 37.7–38.5 cm, 50% proportion of male individuals mature, with 6.3 cm longer for female individuals (L50=44.0–44.8 cm).

Fig. 7.Random effects on parameter K in length-at-age relationships among different models of Lake Erie Walleye Sander vitreus

There were large variations among the regions of Lake Erie to the weight-at-length, length-at-age, and maturity-at-length relationships of Walleye (Figs. 6–8). The main difference was between the central basins(west central basin and east central basin) and other basins (west basin,east basin and Pennsylvanian Ridge). Walleye in the central basins weighted more at the same length and longer in length when at the same age, and matured at much smaller size, than those in other basins. Based onL50, individuals in the central basins (37.5–39.7 cm) matured approximately 4.0 cm smaller than others (41.0–44.5 cm).

The inter-annual and cohort variations of growth and maturity traits were complicated (Figs. 6–8). The random effects estimates ofln(a)exhibited a general decline trend, with the highest values in 1991 and the lowest value in 2014 (Fig. 6).Kestimates kept in relatively low value before 2008 with exceptions in years 1990 and 1991; but were relatively high in the recent decade (Fig. 7). The length of 50% maturity has an overall increasing trend with the highest in 2008, followed by 2009(Fig. 8).

The cohorts born in 1986–1993 were with highera, however, the cohorts born in the next 20 years had relatively lowera(Fig. 6). TheKestimates for the cohorts born during 1986–2003 fluctuated, and then increased from 2003 with a little decrease since 2009 (Fig. 7). The random effects ofL50estimates for cohorts born during 1986–1993, kept increasing from the minimum to the maximum value (Fig. 8). The estimates ofL50were relatively stable and low for cohorts born in recent 2 decades.

Results of PCA revealed that 80% of variance of the 12 local factors and the 4 global environmental factors was explained by the top 5 principal components (Table 4), while the first and second principal component (PC1 and PC2) explained 24% and 21% respectively. The loadings of precipitation and the water supply of Lake Erie to PC1 were positive (0.40 and 0.39, respectively), while the loadings of the fishing mortality, commercial and sport fishery effort for Walleye were negative to PC1 (−0.39, −0.35 and −0.29 respectively, Table 3). The loadings of the abundance of Walleye and Yellow Perch in Lake Erie to PC2 were 0.31 and 0.40 respectively, while the ice melt date, ice cover days, and the NAO were −0.42, −0.38 and −0.30 respectively (Table 3).

Linear regressions suggested that survey year and birth year specificln(a),KandL50of Lake Erie Walleye were influenced significantly (P<0.05) by different principal components. PC1 significantly influenced all the life history parameters in different survey years or cohorts, exceptln(a)in different survey years (Figs. 9 and 10). The significant influences of the principal components to the life history parameters in different survey years were positive (Fig. 9). The influences onln(a)of different cohorts from the PC2 and PC4 in the birth year were significantly positive, while those from the PC1 and PC3 were significantly negative(Fig. 10).

When PC1 is higher, representing higher water supply and lower fishing pressure on Lake Erie Walleye fishery, Walleye grow faster,mature at smaller size and become slimmer. When the PC2 is higher,indicating higher abundances of Walleye and Yellow Perch and less ice cover, Walleye born in those years tend to have better body condition shown as higherln(a). The higher value of PC3, meaning higher water temperature in both Lake Erie and the western coast of North American(positive PDO), lower sport fishing effort, and lower abundance of Walleye, would cause Walleye to grow faster but slimmer. When the PC4(the mixture of the global factors) is higher, meaning that below normal average temperature in the western United States but above normal average temperature in the eastern United States (positive SOI, negative PDO, positive NAO and negative PNA), and higher sport fishing effort,Lake Erie Walleye has better body condition (higherln(a)).

Fig. 8.Random effects on intercept in maturity-at-length relationships among different models of Lake Erie Walleye Sander vitreus. Random effects were applied to intercept (-s ×L50) in the model fitting for maturity-at-length, and need to be divided by (-s) to attain the variations for L50.

Table 4 The importance of the first 5 principal components (PC) from the principal component analysis (PCA) for the biotic and abiotic factors.

4.Discussions

This study indicated that the growth and maturity of Walleye in Lake Erie are influenced by multiple factors. It revealed substantial sexual,temporal and spatial variation in the growth and maturity of Lake Erie Walleye. Male Walleye grow slower but mature at smaller size than the females, and Walleye individuals in the central basins grow faster and mature at smaller size. The water supply, the fishing pressure of Walleye, and the global climate factors significantly influenced the growth and maturity of Walleye during the study period. This study should contribute to the future interpretation of the growth and maturation variation and population dynamics for this important species.

Compared with the previous records of estimates (0.00617(0.00439–0.00865) fora, 3.15 (3.05–3.25) forb) of Walleye in the FishBase (Froese et al., 2014), our estimatesawere lower but the estimate ofbwere higher. Due to the negative correlation of the estimatedaandb, these parameters must be compared together jointly instead of separately. The estimated weight of Walleye in this study was consistent with those of the records given the same length. The estimateL∞was 61.5 cm, much smaller than the maximum total length 78.7 cm, and previous studies documented in the FishBase showed a wide range ofL∞(54.4–104.0 cm). The variation range forKdocumented in the FishBase was 0.16–0.47/year, covering the estimate 0.370/year in this study. The length at 50% mature for Walleye was 41.1 cm, and according to the length-at-age and weight-at-length relationships, the age and weight at 50% were 2.02 year and 663 g, respectively. There is a positive relationship amongln(a) andL50for Lake Erie Walleye (Figs. 6 and 8),indicating the consistent trend of individual growth condition and maturity. The relationship amongln(a) andKis not very clear,i.e.negative among sexual, positive among regions, complicated among survey and born years (Figs. 6 and 7).

When the random effects were applied to the growth and maturity model, only one parameter in each model was selected, i.e.ain weightat-length relationship,Kin length-at-age relationship andL50in the maturity-at-length model. One reason for this was the correlation among parameters inside these models, for example there are a high negative correlation betweenaandbin the modelW=aLband strong correlations forKwithL∞andt0in von Bertalanffy model (Gallucci & Quinn,1979; Pilling et al., 2002; Helser & Lai, 2004). Another reason is that some parameters were used to exhibit the characteristics for the species rather than the individuals’ growth, such asb; additionally, the extrapolation parameter from the available data such ast0, was difficult to interpret. Similar with this study, in the research for the growth variations of 25 exploited marine fishes, the random effects from year,age, or cohort were applied to parametersaandK, whilebwas set to be 3, andL∞andt0were calculated by the maximum and minimum values of length and age for computational convenience (Thorson &Minte-Vera, 2016). Actually, instead of the life history parameters, the mean weight or length at age can be used to evaluate the growth variation, such as in the study of the temporal and spatial variability in size-at-age for Largemouth BassMicropterus salmoidesand Pacific HalibutHippoglossus stenolepis(Helser & Lai, 2004; Kruse et al. 2016).Hence, alternative approaches are available to conduct such kind of research.

Fig. 9.The scatter plots between estimated life history parameters in different survey years (M7) and the principal components (PC) that have significant relationships (P <0.05) with the life history parameters. P values were derived from the linear regressions between the parameter and PCs.

The sexual variation of fish growth and maturity have been confirmed by many studies (Dougherty et al., 2007). The female Walleye in Ottawa River exhibited lowerKbut much longerL∞than the male(Haxton, 2015); the studies for 73 Walleye populations in Ontario revealed that the females matured later and attained a larger asymptotic length with lower mortality than the male Walleye (Lester et al., 2000).The undeveloped Walleye population in Big Crooked Lake, Wisconsin exhibit much later 50% maturity age (3.89–4.88) of the female during 1997–2003 (Schueller et al., 2005). In this study, the female walleye was indicated to grow larger at the same age and mature at smaller size than the males, which was consistent with previous studies about Walleye in other regions. Scientists have found the existence of sexual size dimorphism (Darwin, 1871), and the female-biased sexual size dimorphism was associated with the natural selection for large size and fecundity(Legrand & Morse, 2000), while the male sexual size dimorphism was attributed to scramble competition preferred small males with increased speed and agility (Sandercock, 2001). This study confirmed the results of sexual size dimorphism for Lake Erie Walleye (Henderson et al., 2003).

In Lake Erie, there are three geographical basins: west, central(divided into west central and east central in this paper) and east basins(containing Pennsylvania Ridge and east basin in the paper). The west basin is the shallowest among Lake Erie, only 7.6 m average deep, followed by the central basin (water depth 10–25 m). The east basin is the deepest one, with water depth mostly among 25–60 m. Different depth and latitude lead to different temperature among basins, while the depth and temperature are important factors for species habitat. The Pennsylvania Ridge, shallower than both sides, separates the central and east basin, being an obstacle for abiotic and biotic components interaction between them. In this case, these basins in Lake Erie are markedly diversity in habitat and trophic state. Walleye collected in different basins might have distinct population characteristics and life history variation,which was also indicated by a tagging study (Vandergoot and Brendon 2014). Stepien and Faber (1998) found that the spawning stocks in the west basin of Lake Erie are genetically distinguishable from the stock in the east basin. In addition, the Walleye stocks in the west and east basins, are demonstrated to have markedly different abundance, mortality and distribution patterns (WTG, 2016). The results from this study indicated that Walleyes in the central basin grow faster and mature at smaller size, meaning that the habitat in central basins might be able to promote the individual growth and maturity of Walleye.

The results of model selection revealed that the mixed effects model with random effects from survey year was better than the mixed effects model with random effects from cohort for the weight at length and maturity at age relationships, but worse for the length at age relationships. This phenomenon indicated that the early development of Walleye can affect the fish growth in the whole life, while the condition and maturity can be influenced more easily by the environment, food availability and population density of the observing year, than the early development. Age at 50% maturity is mainly dominated by the maturity status of age 1 and 2 fishes, which is another reason that maturity is influenced heavier by observation year rather than birth year.

The principal component analysis and linear regression indicated that water temperature, water supply, fishing pressure and global climate factors have a substantial influence on the growth and maturity of Lake Erie Walleye. Multiple studies have well documented that temperature is one of the most important factors, which the growth of fish depends on (Brander, 1995; Neuheimer et al., 2011; P¨ortner et al. 2001).In this study, several factors, related to temperature, such as water temperature, ice coverage, and evaporation, exhibit obvious effects on the growth and maturity of Lake Erie Walleye. By reducing population abundance, the fishing pressure can generate density-dependent compensatory responses in individual growth and maturity (Kuparinen& Meril¨a, 2007; Trippel, 1995). Although global climate changes do not directly influence Lake Erie walleye, such as NAO and El Nino, previous studies have found that they have potential impacts on the freshwater ecosystems, through the influence on its physics, hydrology, chemistry,and biology (Ficke et al., 2007; George et al., 2004; Straile et al., 2003).This study revealed that oscillations in Pacific and Atlantic (SOI, PDO,NAO and PNA) have significant correlations with the growth of Lake Erie Walleye, especially with its body condition.

Fig. 10.The scatter plots between estimated life history parameters in different cohorts (M8) and the principal components (PC) that have significant relationships(P <0.05) with the life history parameters. P values were derived from the linear regressions between the parameter and PCs.

Except for the variables analyzed in this paper, gear selectivity may be another possible source of variation of Lake Erie Walleye’s growth and maturity. To cover samples with the range of fish length groups as wide as possible, gill nets used in the Lake Erie Ontario partnership survey include 14 kinds of mesh sizes (from 32 mm to 152 mm) with different survey effort (Ontario Ministry of Natural Resources and Forestry, 2016). Sample sizes could be different to a degree among sexes,regions, years and cohorts. Additionally, in this long-term period(1989–2015), there might be sampling bias resulting from different survey vessels, stations coverage, gear etc. Ages were determined by either otolith or scale reading, and some technique improvement in these three decades may lead to some uncertainties. These potential uncertainties for Lake Erie Walleye growth and maturity were not within the scope of this study, but they deserve further exploration by some simulation or sensitivity analysis in the future.

The variations among years and cohorts in the growth and maturity of Lake Erie Walleye may need to be considered in its stock assessment and fishery management, because the growth and maturity are important input or sub models in many stock assessment models, which contribute to the stock status evaluation and overall management decision making (Lorenzen, 2016; Thorson et al., 2015). For example,the sexual variations of weight and maturity could influence the calculation of spawning stock biomass which is further used in the stock recruitment and spawner per recruitment models.

Our study suggests that multiple factors influence the variation of growth and maturity, and the stock assessment should consider such variations. The methodology constructed in this study could be applied to explore the heterogeneity and multiple impacting factors for other species in aquatic ecosystems, especially for those commercially important species in Great Lakes, which are under similar impacting factors as Lake Erie Walleye.

Declaration of competing interest

The authors declare that there is no conflicts of interest.

CRediT authorship contribution statement

Qiuyun Ma: Conceptualization, Formal analysis, Methodology,Writing - original draft, Writing - review & editing, Visualization,Software, Funding acquisition. Yan Jiao: Conceptualization, Resources,Validation, Writing - review & editing, Supervision, Funding acquisition. Can Zhou: Software, Writing - review & editing. Yiping Ren:Resources, Writing - review & editing, Validation, Supervision.

Acknowledgement

We gratefully acknowledge Ontario Ministry of Natural Resources and Forestry and Ontario Commercial Fisheries’ Association for providing the partnership survey data. We thank the Department of Fish and Wildlife Conservation at Virginia Polytechnic Institute and State University (Virginia Tech) to host Q. Ma to work on this project under the supervision of Y. Jiao. We also thank the China Scholarship Council who provided the funding for Q. Ma to work at Virginia Tech. This research was also partially supported by a grant “Verification of natural mortality estimation of Walleye in Lake Erie based on integrated Bayesian statistical catch-at-age model” to Y. Jiao at Virginia Tech by the Great Lakes Fishery Commission (459013).

Appendix

Appendix Fig. 1.The temporal trend of the biotic and abiotic factors which may influence the life history of Walleye from 1986 to 2015. Chlorophyll A concentration in years 1987–1995 are unavailable, using the average value of all else years.

Appendix Fig. 2.Variation of ln(weight) residuals in weight-at-length relationships from the global model M1. Basin: west (W), west central (WC), east central (EC),Pennsylvania ridge (PR) and east basins (E). The previous 6 cohorts are 1979, 1980, 1982, 1983, 1984 and 1985, and the next cohorts are sequential without breakpoint.

Appendix Fig. 3.Variation of ln(length) residuals in length-at-age relationships from the global model M1. Basin: west(W), west central (WC), east central (EC),Pennsylvania ridge (PR) and east basins (E). The previous 6 cohorts are 1979, 1980, 1982, 1983, 1984 and 1985, and the next cohorts are sequential without breakpoint.

Appendix Fig. 4.Variation of maturity residuals in maturity-at-length model from the global model M1. Basin: west(W), west central (WC), east central (EC),Pennsylvania ridge (PR) and east basins (E). The previous 6 cohorts are 1979, 1980, 1982, 1983, 1984 and 1985, and the next cohorts are sequential without breakpoint.