Genetic and life‐history changes associated with fisheries‐induced population collapse

Abstract Over the recent years, growing number of studies suggests that intensive size‐selective fishing can cause evolutionary changes in life‐history traits in the harvested population, which can have drastic negative effects on populations, ecosystems and fisheries. However, most studies to date have overlooked the potential role of immigration of fish with different phenotypes as an alternative plausible mechanism behind observed phenotypic trends. Here, we investigated the evolutionary consequences of intensive fishing simultaneously at phenotypic and molecular level in Eurasian perch (Perca fluviatilis L.) population in the Baltic Sea over a 24‐year period. We detected marked changes in size‐ and age‐distributions and increase in juvenile growth rate. We also observed reduction of age at sexual maturity in males that has frequently been considered to support the hypothesis of fisheries‐induced evolution. However, combined individual‐based life‐history and genetic analyses indicated increased immigration of foreign individuals with different life‐history patterns as an alternative mechanism behind the observed phenotypic change. This study demonstrates the value of combining genetic and phenotypic analyses and suggests that replacement or breakdown of locally adapted gene complexes may play important role in impeding the recovery of fish populations.


Introduction
Fishing causes substantial mortality in commercially exploited fish stocks, often exceeding the level of natural mortality (Stokes and Law 2000;Jørgensen et al. 2007;Law 2007). As a consequence, over 80% of the world's fish species are currently fully exploited or overfished (FAO 2011). Moreover, comparison of current population abundance estimates with historical records indicates that heavy exploitation has caused drastic collapses in many commercially harvested marine stocks (Pauly et al. 2002;Myers and Worm 2003;Hutchings and Reynolds 2004;Pauly 2008). Despite various management strategies aimed at reducing fishing pressure, only few populations show signs of recovery in abundance (Stockwell et al. 2003;Hutchings and Reynolds 2004). This suggests that exploitation, and subsequent population recovery, is not simple processes decreasing or increasing biomass, but several environmental, ecological and genetic mechanisms can be involved, such as habitat alterations, shifts in species interactions and changes in population genetic composition as well as the evolution of life-history traits (Hutchings 2005;Enberg et al. 2009;Swain 2011;Kuparinen and Hutchings 2012). Nonetheless, little remains known about the genetic consequences of harvesting, partly because temporal investigations that simultaneously evaluate phenotypic and molecular changes are very scarce ( Arnason et al. 2009;Jakobsd ottir et al. 2011).
Over the last decade, increasing numbers of studies of experimental and wild harvested fish populations have reported significant temporal shifts in life-history traits that have been associated with fisheries-induced selection (Conover and Munch 2002;Jørgensen et al. 2007; Kuparinen and Meril€ a 2007; Allendorf et al. 2008;Allendorf and Hard 2009). At the same time, many studies have also pointed out the difficulties in distinguishing phenotypical responses from true evolutionary changes (Stokes and Law 2000;Kuparinen and Meril€ a 2007;Law 2007;Kuparinen and Meril€ a 2008). For example, faster growth and earlier maturation at a smaller size can be a consequence of fisheries-induced evolution (FIE) or, alternatively, may be explained by a decline in fish density, where reduced competition for food and space leads to increased growth and earlier maturation (Law 2000). Furthermore, only a few of the studies focusing on disentangling genetic and environmental effects of fisheries exploitation have explicitly tested whether the observed phenotypic shifts could be at least partly attributable to population replacements, immigration or emigration (Jakobsd ottir et al. 2011;Hansen et al. 2012). However, while the role of population replacement is a relatively unconsidered concept in fisheries literature, this is not the case for other taxa (e.g. Pergams and Lacy 2008;Mackie et al. 2010;Saltonstall 2002). For example, rapid morphological and genetic changes of white footed mice (Peromyscus leucopus) have been associated with immigration of non-native individuals with a distinct phenotype that displaced the previously resident population due either to selective advantage in the changing habitat or to genetic swamping (Pergams and Lacy 2008). Nevertheless, given that the majority of fish species consist of reproductively isolated populations and genetic differences often exist at a much smaller scale than the dispersal ability of the species (Allendorf et al. 2008;Hauser and Carvalho 2008;Bergek and Bj€ orklund 2009), the influx of non-native individuals with different geno-and phenotypes could be more widespread than previously anticipated. For instance, temporal genetic analysis of North Sea cod (Gadus morhua) demonstrated that a heavily exploited population can be strongly affected by gene flow from neighbouring populations (Hutchinson et al. 2003). Similarly, drastic changes in population genetic composition have been well documented in many fish species, particularly in systems with prolonged and/or strong anthropogenic impacts (Vasem€ agi et al. 2005;Hansen 2008, Palstra andRuzzante 2010) or seasonal variations in abiotic conditions (Crispo and Chapman 2009).
Here, we investigated the consequences of overexploitation and subsequent population recovery of the Eurasian perch (Perca fluviatilis L.) population in Matsalu Bay, in the eastern Baltic Sea (Fig. 1), based on combined temporal genetic and life-history analyses over a 24-year period. The system provided a unique opportunity for studying the effects of intensive fishing as prior to the 1990s, when Estonia was still part of the Soviet Union, only commercial fishermen had rights to fish in the study area with trap and gill nets. However, after the Republic of Estonia reclaimed its independence from the Soviet Union in 1991, fishing became virtually unregulated for a short period of time in many coastal areas of Estonia, which encouraged overexploitation and lead to the subsequent collapse of the Matsalu Bay perch population in the late 1990s ( Fig. 2A; J€ arv 2002; J€ arv et al. 2005). Despite multiple fishing regulations that were introduced in 1999, the Matsalu Bay perch population has not recovered and the population size has remained very low (J€ arv et al. 2005).

DNA extraction and microsatellite genotyping
Total DNA from dried scale and operculum samples was extracted following a salt-based protocol combined with a vacuum-based extraction method as described by Swatdipong et al. (2010a). For ethanol-preserved samples (P€ arnu Bay, 2010), DNA was isolated from fin clips according to Laird et al. (1991). Genetic diversity was initially assessed using 11 microsatellite markers: Svi6, Svi18 and Svi33 (Borer et al. 1999) and PFE01, PFE03, PFE11, PFE12, PFE14, PFE15, PFE19 and PFE22 (Zhan et al. 2009). Data for four loci were subsequently excluded due to low levels of polymorphism (Svi33, PFE14, PFE22) or due to poor amplification (PFE15).

Phenotypic analysis
Possible changes in fish life-histories were investigated by comparing life-history data of perch from the pre-and post-collapse periods in Matsalu Bay. Growth data were available only for the pre-collapse years 1987 and 1990 and for the post-collapse years 2009 and 2010. In these two year groups, changes in age-specific lengths were investigated using linear mixed effect models (LME), where sex and a factor indicating whether the individual belonged to the 'pre-' or 'post-' collapse group, as well as their interactions, were set to fixed effects and the cohort was set to a random effect. These analyses were conducted separately for ages from 1 to 5 years, as no older individuals were present in the post-collapse period, and thus comparison between pre-and post-collapse groups was not possible. Maturity data were available for the pre-collapse years 1987-1992 and post-collapse years 2006-2008 and 2010. Sex-specific maturity ogives (i.e. the age-specific probability of being mature) in pre-and post-collapse years were investigated through a generalized linear model (GLM) with a binomial error structure, having age, length, collapse indicator (factor, similar to that described above) as well as their two-way interactions and the quadratic effects of age and length as effects. Model simplifications were made by stepwise model reduction based on either likelihood ratio tests (LME) or chi-squared tests (GLM) as suggested by Crawley (2007). Whole life-history analyses were conducted using R 2.11.1 (R Development Core Team 2010).

Population genetic analysis
GenePop 4.0 (Rousset 2008) was used to test for deviations from the gametic phase and Hardy-Weinberg (HW) equilibria (5 000 iterations each). To assess levels of genetic diversity over time, allelic richness (A r ; the effect of different sample sizes are taken into account), mean number of alleles (A), observed (H o ) and expected heterozygosity (H e ) were calculated using FSTAT 2.9.3.2 (Goudet 1995) and Microsatellite Toolkit 3.1.1 (Park 2001). To identify potential footprints of genetic bottleneck, the M-ratio measure was calculated following Garza and Williamson (2001). The standard deviation of the M-ratio was calculated by bootstrapping over loci (10 000 replications) as described in Swatdipong et al. (2010b) using the program PopTools 2.7.5 (Hood 2006).
To investigate genetic differentiation between year classes, pairwise F ST values were calculated following Weir and Cockerham (1984) and implemented in FSTAT 2.9.3.2 (Goudet 1995). The significance of allele frequency differences between cohorts was estimated using the genetic differentiation test implemented in GenePop 4.0 (Rousset 2008).
The program Structure 2.3 (Pritchard et al. 2000) was used to detect changes in genetic composition over the period of 24 years based on multilocus genotype information. This approach aims to divide the genotype data into homogenous genetic clusters (K) that are in HW and gametic phase equilibrium. Twenty independent runs were performed for each value of K (K = 1-6), using prior population information, as suggested by Hubisz et al. (2009), with a burn-in period of 100 000 iterations and 100 000 replications. The method developed by Evanno et al. (2005) was used to estimate the most likely number of clusters based on the rate of change in the logarithmic probability of data between successive K-values. In addition to Matsalu Bay samples, individuals from nearby regions (P€ arnu Bay and Moonsund Archipelago Sea) were also included in the analysis. To avoid spurious clustering of the individuals due to incomplete genotype data, individuals with more than two missing loci were excluded from this analysis (693 individuals retained).
The effective population size (N e ) and immigration rate (m) of the Matsalu Bay perch population were estimated using a maximum-likelihood method implemented in the program MLNE 1.0 (Wang and Whitlock 2003). This method requires a user-specified upper limit for N e , which was set at 2 000 after checking that similar results were obtained with higher upper limits for N e . To obtain estimates of N e and m for the pre-and post-collapse periods the data were divided into two groups (from cohorts 1981 to 1990 and from cohorts 1995 to 2008, respectively). We subsequently created a putative donor population (source of immigrants) by pooling all the individual genotypes that clustered together with more than a 0.8 probability as 'foreign' genotypes, identified by Structure (cohorts 1987, n = 3; 1999, n = 1; 2005, n = 6; 2006, n = 41; 2007, n = 20; 2008, n = 7). We also carried out analyses using samples from neighbouring regions (P€ arnu, Virtsu, Paslepa and Haapsalu Bay) as potential source populations (data not shown). The generation length for pre-and post-collapse cohorts was set at 5.3 and 3.3 years, respectively, according to the mean ages of mature individuals. As the intervals between year classes were not integers, all estimates were adjusted according to the equations provided by Wang and Whitlock (2003). More specific estimations of N e and m were also carried out, keeping pre-and postcollapse samples separate. For pre-collapse analysis, the 1981 and 1982 samples were pooled together (0 generation) and the subsequent year classes were analysed separately until 1990 (1.7 generations). Similarly, for post-collapse analysis the 1989 and 1990 samples were pooled together (0 generation) and the subsequent year classes were analysed separately until 2008 (4.76 generations).

Analysis of life-history traits based on multilocus genotype information
To test whether the individuals belonging to distinct genetic clusters identified by Structure differed in their growth rate, the age-specific lengths (ages from one to five) were used in linear mixed models, having sex, cluster, and their interaction as fixed effects, and cohort as a random effect. Similar to the analysis described earlier, model simplification was made by stepwise model reduction based on likelihood ratio tests (LME) using R 2.11.1.

Changes in growth rate and maturation
The mean age and length of perch in the Matsalu Bay population decreased considerably during the study per-iod (Fig. 2B). In the pre-collapse years (1987)(1988)(1989)(1990)(1991)(1992), the mean age and length of perch were 5.1 years and 195.3 mm respectively. However, in the post-collapse years, the mean ages and lengths were considerably reduced (mean age 3.0 years; mean length 164.7 mm). While such a pattern can result from a shift in the demographic structure of the population, analysis of the agespecific length showed that the growth rates of individuals had changed over time. One-year-old perch in 2009 and 2010 were on average 17.5 mm longer (Fig. 2C) (likelihood ratio (L) = 16.3; P < 0.001) than fish of the same age in the pre-collapse years (1987 and 1990), while 2-year-old perch were 10.9 mm longer in the post-collapse years (L = 11.8, P < 0.001). As shown in Fig. 2C, most of the difference in juvenile length-at-age can be attributed to lower growth rates during the first year of age. No significant differences in length at the ages of 3, 4 or 5 years were found between the pre-and post-collapse groups (group indicator could be reduced from the model with P > 0.05).
In the pre-collapse samples, 40% of 2-year old, 37% of 3year old and 10% of 4-year old males were immature, while in the post-collapse group all the studied males at the ages of two (n = 22), three (n = 54) and four (n = 33) were mature. In males, the probability of being mature (maturity ogive) was found to differ significantly between the preand post-collapse groups (v 2 =À29.9, P < 0.001) as well as to depend on the age of the individual (v 2 =À17.2, P < 0.001; slope estimate was 0.8193 (z = 3.862, P < 0.001)). In females, on the other hand, the probability of being mature did not differ between the pre-and postcollapse groups (v 2 =À0.262, P = 0.608), but was significantly affected by age (v 2 =À146.2, P < 0.001; slope estimate was 1.9370 (z = 7.754, P < 0.001)).

Deviations from Hardy-Weinberg and gametic phase equilibria
Hardy-Weinberg (HW) equilibrium testing over all loci indicated that five out of 14 cohort samples from Matsalu Bay (1987, 1995, 2006, 2007 and 2008) deviated significantly from HW equilibrium proportions (P < 0.05 before Bonferroni correction). However, after applying the Bonferroni correction, deviations in three recent cohorts (2006-2008) remained significant ( Table 1). The linkage disequilibrium test indicated that in the post-collapse cohorts (2006)(2007)(2008) a large proportion of locus pairs significantly deviated from the gametic phase equilibrium (P < 0.05; 7-10 locus pairs of 21 comparisons after the Bonferroni correction; Table 1). In contrast, just one locus pair deviated from gametic phase disequilibrium in eleven of the Matsalu Bay cohorts sampled before the population collapse.

Temporal patterns of genetic diversity
Genetic diversity measures (H o , H e , A, A r ) from Matsalu Bay showed considerable fluctuations over the period of 24 yearsthe largest changes, however, were associated with the more recent post-collapse cohorts born in 2006-2008 (Fig. 3). Both mean expected heterozygosity (H e ) and allelic richness (A r ) were lowest in the 2006 cohort, whereas the mean number of alleles (A) and observed heterozygosity (H o ) had the smallest values in the 2008 year class (Fig. 3A,B). M-ratio estimates, on the other hand, showed no marked decrease in post-collapse cohorts ( Table 1), suggesting that a sharp decrease in population size did not result in obvious signs of genetic bottleneck.

Bayesian clustering and analysis of life-history traits based on multilocus genotype information
Bayesian clustering analysis showed that the temporal samples from Matsalu Bay and nearby regions can be separated into two distinct genetic groups (Fig. 4). These two clusters were identified using the DK method (Evanno et al. 2005) which had highest mean log likelihood of À9886.41 for K2 (À10134.20 = K1, À9941.86 = K3, À9931.12 = K4, À10012.40 = K5 and À10116.90 = K6). This did not distinguish populations based on their geographical proximity but instead, suggested that the majority of post-collapse individuals spawning in Matsalu Bay (cohorts 1999-2001, 2005-2008) have distinct multilocus genotypes compared to cohorts sampled earlier from the same location. In addition to the post-collapse samples, a few individuals in the 1986 and 1987 cohorts also possessed multilocus genotypes that resembled post-collapse samples (Fig. 4).
Analysis of the age-specific length in post-collapse cohorts (taking into account the Bayesian clustering information) showed that the individuals belonging to distinct genetic clusters also differed in their growth rates. The length of 1-year-old fish depended on sex, inferred cluster and the interaction of sex and inferred cluster (P < 0.05 in all model reductions). The length of 1-year-old females carrying the local Matsalu Bay multilocus genotype was on an average 22 mm smaller than for fish carrying a foreign genotype in the post-collapse sample (SL 70 mm and 92 mm respectively). Similarly, local Matsalu Bay males were smaller than the foreign individuals of the same age (SL 76 mm and 80 mm respectively). At the age of 2 years, the length did not differ between sexes, but perch carrying the local Matsalu genotype were still 16 mm smaller than the individuals of foreign origin (likelihood ratio = 15.2, P < 0.001). At the ages of 3, 4 and 5 years, however, no significant differences between the fishes of local and foreign origin were found (P > 0.05 in all model reductions). In contrast to the analysis of the whole data set (Fig. 2C), only a minor increase in mean length was observed at the age of one year among local fish in the post-collapse years, when individuals with non-native genotypes were excluded (Fig. 2D) (likelihood ratio = 5.9, P < 0.015; after Bonferroni correction, P > 0.05). At the ages of 2-5 years, no differences in length between pre-and post-collapse samples with local genotype were found (P > 0.05 in all model reductions).

Discussion
Both theoretical and empirical studies suggest that intensive fishing can cause changes in the life-history traits of the harvested population, yet distinguishing possible evolutionary changes from plastic environmental variations is far from trivial (e.g. Stokes and Law 2000;Heino and Godø 2002;Kuparinen and Meril€ a 2007;Nussl e et al. 2009;Uusi-Heikkil€ a et al. 2010;Kuparinen et al. 2011). In this study, we detected marked changes in size-and age-distributions and increases in juvenile growth rates and reductions in age at sexual maturity in males, parameters which have frequently been considered as evidence supporting a FIE scenario. We were not able to estimate probabilistic maturation reaction norms for post-collapse data as all the fish at age two or older were already mature in contrast to pre-collapse data, where ca 40% of 2-and 3-year old were immature. However, based on combined individual-based life-history and genetic analysis, the findings of this study provide empirical evidence that temporal phenotypic trends can also arise from increased immigration of foreign individuals with different life-history patterns.
The results of the present study are in accord with recent findings in Atlantic cod, demonstrating that fishing can differentially target genetically and phenotypically distinct components of a population, thus leading to phenotypic changes associated with a change in the population composition ( Arnason et al. 2009;Jakobsd ottir et al. 2011). Such genetic shifts would remain unnoticed without the use of molecular genetics tools (Hutchinson et al. 2003) and therefore this study, similar to others ( Arnason et al. 2009;Jakobsd ottir et al. 2011), demonstrates the value of combined genetic and life-history analyses when studying the effects of population recovery and harvesting-induced selection. When analysing the possible causes behind the changes in life-history traits, it is therefore important to rule out population replacement, emigration and immigra-tion scenarios as alternative mechanisms behind phenotypic change (Jakobsd ottir et al. 2011;Hansen et al. 2012). The dispersal of individuals may also be non-random in relation to genotype, individual state and environmental factors (Bowler and Benton 2005;Edelaar and Bolnick 2012). Nevertheless, further research is needed to explicitly test for the presence of non-random dispersal and gene flow in perch in Matsalu Bay. Currently, there exists a rather limited number of temporal genetic studies of exploited marine fishes (e.g. Hutchinson et al. 2003;Nielsen and Hansen 2008;Arnason et al. 2009;Larsson et al. 2010;Palstra and Ruzzante 2010), and it is possible that many extinctions of local fish populations have gone unnoticed and/or undocumented because of population replacements through immigration (Vasem€ agi et al. 2001(Vasem€ agi et al. , 2005. Future genetic studies using archived material and new genomic tools could shed light onto this critical but underrated issue (Nielsen and Hansen 2008;Palstra and Ruzzante 2010).
This study also provided insights into effective population size associated with the collapse and recovery of a targeted population. In the Matsalu Bay perch population, a drastic reduction in the effective population size was observed, as N e in the post-collapse period (cohorts 1995-2008) was seven times smaller than in the pre-collapse period (cohorts 1981-1990). Contrary to expectations, this sharp decrease in N e did not result in a rapid loss of genetic diversity. This can be explained by the increased influx of non-native immigrants, which compensated for the expected loss of genetic diversity among the remnant native population. This increased influx of immigrants is not surprising, given that a reduced N e can make the local populations more sensitive to immigration, potentially leading to genetic swamping and/or loss of local adaptations (Hutchinson et al. 2003;Pergams and Lacy 2008;Allendorf et al. 2008;Mackie et al. 2010;Edelaar and Bolnick 2012).
before the overexploitation of the Matsalu Bay perch population (e.g. 1986: m = 0.276, 95% CI = 0.157-0.404; Fig. 4). It is important to keep in mind, however, that the estimates of m are obtained by pooling immigrant genotypes identified with Structure to establish allele frequencies of a putative donor population. Therefore, estimated migration rates should be regarded with caution until identification of the actual source population of immigrants. Analysis of additional four sampling sites indicated that these geographically close coastal areas (Moonsund Archipelago Sea and P€ arnu Bay) most likely do not represent the source population(s) of immigrants that were found in Matsalu Bay. We therefore suggest that the immigrant genotypes observed in our study area might originate from near-shore littoral perch from the River Kasari delta area that flows into Matsalu Bay, while the vast majority of genotypes from the pre-collapse samples belong to the more open-water pelagic form. Hence, the observed phenotypic and genotypic changes may be related to expansion or decrease of the spatial shift in the spawning areas of these populations. The occurrence of genetically and phenotypically distinct but geographically close populations in Eurasian perch is also supported by earlier research. For example, lake and river populations may exchange many migrants due to their proximity (Hendry et al. 2002), also perch from littoral and pelagic zones within the same lake can exhibit significant morphological and growth differences (Svanb€ ack and Ekl€ ov 2002; Bartels et al. 2012) and evidence for stable genetic differentiation at a small spatial scale have been found both in lakes (Gerlach et al. 2001) and in the Baltic Sea (Bergek and Bj€ orklund 2009). As a result, we reiterate that despite the uncertainty regarding the origin of the immigrants, migration remains the most plausible mechanism responsible for the sharp change in genetic composition and marked increase in the frequency of observed Hardy-Weinberg and genotypic equilibrium deviations in the post-collapse samples collected from Matsalu Bay.
Similar to many earlier studies that focused on the consequences of overfishing, this study describes a rapid collapse of a local population most likely as a result of over-exploitation, followed by a slow recovery despite of the reductions in fishing pressure (Pauly et al. 2002;Myers and Worm 2003;Hutchings and Reynolds 2004;Olsen et al. 2004Olsen et al. , 2005Pauly 2008). Intriguingly, this reflects the current status of many fish stocks that, in spite of reduced fishing pressure, struggle to recover (Hutchings 2000;Stockwell et al. 2003;Hutchings and Reynolds 2004). It has been suggested that the sharp decrease in the abundance and the slow recovery of perch in the Moonsund Archipelago might be at least partly related to the growing cormorant (Phalacrocorax carbo sinensis) colony (Vetemaa et al. 2010), but as only small proportion of the cormorant diet consists of perch (6-10%), increased predation by cormorants most likely does not play a major role in the slow recovery of perch in Matsalu Bay. As an alternative explanation, it is possible that the increase of non-adaptive alleles and the breakdown of locally adapted assemblages may have had a negative effect on the recovery of the Matsalu Bay perch population. This is in line with an increasing number of studies in both marine and freshwater fishes that show marked population-level differences in ecologically important traits, suggesting that adaptive population differentiation is likely to be more widespread than previously thought (Hauser and Carvalho 2008). On the other hand, dispersal and gene flow can either promote or constrain adaptive divergence through either genetic or demographic routes (Lenormand 2002;Garant et al. 2007). For example, immigration of maladapted individuals may decrease the mean fitness (Storfer 1999), while immigration can have a healing effect that counteracts inbreeding depression in small populations (Ingvarsson and Whitlock 2000). Future research, combining genomic tools and temporal monitoring of life-history traits of the exploited fish stocks, is clearly needed to further understand the importance of local adaptation in the field of fisheries science and conservation.

Data Archiving Statement
Data deposited in the Dryad repository: doi:10.5061/dryad. td237