Received: 20 April 2016

|

Revised: 26 April 2017

|

Accepted: 28 April 2017

DOI: 10.1111/mec.14179

ORIGINAL ARTICLE

Rapid, broad-scale gene expression evolution in experimentally harvested fish populations Silva Uusi-Heikkil€ a1

| Tiina S€ avilammi1 | Erica Leder1,2,3 | Robert Arlinghaus4,5 |

Craig R. Primmer1 1

Department of Biology, University of Turku, Turku, Finland 2

Natural History Museum, University of Oslo, Oslo, Norway 3

Department of Aquaculture, Institute of Veterinary Medicine and Animal Sciences, Estonian University of Life Sciences, Tartu, Estonia 4

Department of Biology and Ecology of Fishes, Leibniz-Institute of Freshwater Ecology and Inland Fisheries, Berlin, Germany 5 Division of Integrative Fisheries Management, Department of Crop and Animal Sciences, Faculty of Life Sciences, Humboldt-Universit€at zu Berlin, Berlin, Germany

Correspondence Silva Uusi-Heikkil€a, Division of Genetics and Physiology, Department of Biology, University of Turku, Turku, Finland. Email: [email protected]. Present address Craig R. Primmer, Department of Biosciences, and Biotechnology Institute, University of Helsinki, Helsinki, Finland Funding information € ; LeibnizAcademy of Finland; Koneen S€a€atio Gemeinschaft; AXA Research Fund; ICES

Abstract Gene expression changes potentially play an important role in adaptive evolution under human-induced selection pressures, but this has been challenging to demonstrate in natural populations. Fishing exhibits strong selection pressure against large body size, thus potentially inducing evolutionary changes in life history and other traits that may be slowly reversible once fishing ceases. However, there is a lack of convincing examples regarding the speed and magnitude of fisheries-induced evolution, and thus, the relevant underlying molecular-level effects remain elusive. We use wild-origin zebrafish (Danio rerio) as a model for harvest-induced evolution. We experimentally demonstrate broad-scale gene expression changes induced by just five generations of size-selective harvesting, and limited genetic convergence following the cessation of harvesting. We also demonstrate significant allele frequency changes in genes that were differentially expressed after five generations of sizeselective harvesting. We further show that nine generations of captive breeding induced substantial gene expression changes in control stocks likely due to inadvertent selection in the captive environment. The large extent and rapid pace of the gene expression changes caused by both harvest-induced selection and captive breeding emphasizes the need for evolutionary enlightened management towards sustainable fisheries. KEYWORDS

adaptation, captive breeding, fisheries-induced evolution, gene expression, size selection, transcriptome sequencing

1 | INTRODUCTION

has been studied less intensively (Fraser, Moses, & Schadt, 2010; McKinney et al., 2015). Because gene expression can be regulated at

Gene expression changes play an important role in adaptive evolution

many steps (from DNA to protein), and the regulatory networks gov-

and have been a topic of great interest in recent years (Jones et al.,

erning gene expression are complex and modular, gene expression

2012; Manceau, Domingues, Mallarino, & Hoekstra, 2011; Papakostas,

changes can be more labile than DNA sequence changes (Huminiecki

Vasemagi, Himberg, & Primmer, 2014). A number of studies have doc-

& Wolfe, 2004). Therefore, it has been suggested that gene expression

umented that the role of gene expression changes in adaptation of

changes may be particularly important for fostering (or constraining)

phenotypic traits is influenced strongly by just one or few genes with

genetic adaptation to strong selection pressures over short timescales

large effect in a broad range of species (Manceau et al., 2011; Miller

(Jones et al., 2012; Papakostas et al., 2014).

et al., 2007; Shapiro et al., 2004). However, the extent to which adap-

Human activities, such as habitat fragmentation, harvesting and

tations in polygenic traits are mediated by gene expression changes

domestication, are affecting countless animal and plant species

3954

|

© 2017 John Wiley & Sons Ltd

wileyonlinelibrary.com/journal/mec

Molecular Ecology. 2017;26:3954–3967.

€ UUSI-HEIKKILA

|

ET AL.

3955

around the globe by changing the environment in which they occur

populations, particularly in freshwater environments (Lorenzen, Bev-

and creating new adaptive landscapes (Palumbi, 2001; Sullivan, Bird,

eridge, & Mangel, 2012). However, in many cases, whether or not

& Perry, 2017). For example, wildlife harvesting not only exceeds

hatchery-reared fish actually help population recovery has been ques-

natural mortality rates (e.g., Jørgensen et al., 2007), but it is also

tioned (Araki, Cooper, & Blouin, 2007; Ryman & Laikre, 1991). This is

selective in terms of sex, age, size and a range of other phenotypes

mainly because unintentional selection pressures brought about by

(Allendorf & Hard, 2009; Arlinghaus et al., 2017; Kuparinen & Festa-

artificial breeding and captive rearing can induce genetic and pheno-

Bianchet, 2017). A classic example of selective wildlife harvesting is

typic changes that are maladaptive in the natural environment leading

fishing, where the largest individuals, cohort after cohort, are selec-

to high postrelease mortality and low reproductive success (Araki

tively removed from the population, often at very high rates. The

et al., 2007; Christie, Ford, & Blouin, 2014; Lorenzen et al., 2012).

resulting selection against large body size is predicted to favour fast

While the fitness consequences of captive breeding on wild popula-

life histories characterized by the evolution of earlier maturation at

tions have been demonstrated, much less is known about the molecu-

smaller size and higher energy investment into reproduction, which

lar basis of inadvertent selection in captivity (but see Christie, Marine,

reduces postmaturation growth (Heino, Dıaz Pauli, & Dieckmann,

Fox, French, & Blouin, 2016). A better understanding of the broad-

2015; Stearns, 1992). Despite the accumulating empirical evidence

scale effects of captive breeding at the gene expression level, and

from observational phenotypic data collected in the wild (Jørgensen

identifying the genes upon which domestication selection acts could

et al., 2007), there has been a long controversy on the question of

enable more appropriate recovery programmes and management

whether phenotypic changes in time series are truly evolutionary or

practices for exploited and endangered fish populations.

merely reflect plastic life history responses (Browman, Law, & Mar-

We studied changes in gene expression at the genomewide level

shall, 2008; Hilborn & Minte-Vera, 2008) and what is the rate of

in experimental fish populations subjected to five generations of

fisheries-induced evolution (Andersen & Brander, 2009). Further-

either size-selective (large-harvested fish) or random (random-har-

more, the impact of fisheries-induced selection on gene expression

vested fish) harvesting followed by a six generation no-harvest per-

(how many genes and how fast expression changes) is unknown.

iod simulating a fishing moratorium. We sought to address three

This is a major shortcoming for understanding the molecular mecha-

questions (Figure 1): (i) How rapidly and to what extent do experi-

nisms involved in harvest-induced evolution.

mental fish populations respond to size-selective exploitation at the

Another question, which has been rarely studied empirically, is

gene expression level? (ii) How do the experimentally exploited fish

whether or not populations can recover from harvest-induced evolu-

populations respond to the cessation of harvesting? and (iii) To what

tionary changes (e.g., following a fishing moratorium) and if so, at

extent does unintentional selection induced by captive breeding

what pace? Size-selective harvesting can be a more intensive selec-

result in a response at the gene expression level?

tion force than natural selection (e.g., predators, competitors; Darimont et al., 2009). Therefore, populations might not adapt to the cessation of size-selective harvesting as rapidly as they adapt to the fisheries selection (Conover, Munch, & Arnott, 2009; Enberg, Jorgensen, Dunlop, Heino, & Dieckmann, 2009). Also, intensive selec-

2 | MATERIAL AND METHODS 2.1 | Experimental design

tion together with high mortality (and genetic drift) might deplete

We used zebrafish (Danio rerio) caught from the wild (in West Ben-

population genetic diversity (Marty, Dieckmann, & Ernande, 2015;

gal, India in 2006) as a model species to study transcriptome-wide

Pinsky & Palumbi, 2014). This can reduce the adaptive potential of a

gene expression changes in response to five generations of size-

population and hamper its recovery (Allendorf, England, Luikart,

selective harvesting (Uusi-Heikkil€a et al., 2015), followed by six gen-

Ritchie, & Ryman, 2008). There are empirical examples of wild fish

erations of no harvesting, and in response to up to 11 generations

populations not fully recovering demographically after the cessation

of captive breeding. The founder fish were bred randomly for one

of fishing (Neubauer, Jensen, Hutchings, & Baum, 2013) potentially

generation to enable acclimation to captive conditions and reduce

due to harvest-induced decline in genetic variation (Swain, Sinclair, &

possible parental effects before the selection experiment started.

Hanson, 2007); however, no studies have examined the recovery

Briefly, we subjected replicated (two replicates per harvest treat-

pattern at the molecular level. Understanding the extent and mecha-

ment, each consisted of 450 individuals) wild-origin zebrafish popula-

nisms of evolutionary responses in exploited populations is crucial

tions to selection for small body size (hereafter referred to as “large-

for the design of management tools aimed at sustainable exploitation

harvested”, which is the treatment most similar to fisheries), while

of natural biological resources.

we harvested the control lines randomly with respect to body size

Many fisheries are experiencing overexploitation and conse-

(“random-harvested”). For the large-harvested fish, we applied a 75%

quently temporary moratoria or no-take reserves are being estab-

per-generation harvest rate (based on standard length) and removed

lished to help fish stocks to replenish. As a moratorium can be

the largest 75% of individuals (i.e., 338 individuals) from the popula-

considered as a rather extreme management practice to facilitate the

tion. This is consistent with intensive lethal capture fisheries (Lewin,

recovery of an overexploited fish stock, supplementary stocking (i.e.,

Arlinghaus, & Mehner, 2006), thus the large-harvested fish repre-

the release of hatchery-reared fish to supplement the natural popula-

sented a harvest scheme occurring in highly exploitive capture fish-

tion) is more commonly applied to restore and conserve wild

eries. In the random-harvested line, the same proportion (75%) of

3956

€ UUSI-HEIKKILA

|

ET AL.

Generation

F2

F5

F11

Large-harvested

Random-harvested

F I G U R E 1 A conceptual figure of the experimental design and the research questions. The effects of size-selective harvesting on gene expression were studied by comparing the expression patterns of large- and random-harvested fish after five generations of harvesting (red line). The effects of a no-harvest period on gene expression were studied in two ways, first by comparing the expression patterns of large- and random-harvested fish after six generations of no harvesting (dashed orange line) and second, by comparing the overall expression changes (from F2 to F11: thick solid orange line) to expression changes during the harvesting period (from F2 to F5: thin solid orange line) among the large-harvested fish. The effects of captive breeding on gene expression were studied by comparing expression patterns of the F2- and F11generations among the random-harvested fish and comparing the changes in both harvest treatments during the whole experimental period (from F2 to F11; grey lines) [Colour figure can be viewed at wileyonlinelibrary.com]

individuals was harvested in each generation, but randomly with

sampling, fish had been reared and bred in the laboratory for three

respect to body size. In both harvest treatments, the remaining 25%

generations (and harvested for two generations). F2-generation sam-

of females and males (i.e., 112 individuals) were used for reproduc-

ples were the earliest samples available from the selection experi-

tion (Uusi-Heikkil€a et al., 2015). The harvesting regime was con-

ment that were suitable for RNA analyses. Liver was chosen as a

ducted for five generations (F1–F5) and then halted for an additional

target organ for analyses because of its important role in various

six generations (F6–F11), hereafter referred to as the “no-harvest per-

major physiological processes in fish, such as skeletal and soft tissue

iod.” During the no-harvest period, the experimental populations

growth, lipid, protein and carbohydrate metabolism, energy storage,

consisted of 110–120 individuals. As fish were not harvested during

maturation, reproduction and development (Devlin, Sundstrom, &

this period, all fish were allowed to reproduce. Large- and random-

Muir, 2006). Additionally, many hormones and hormone receptors

harvested fish were reared in a common-garden environment, includ-

involved in the regulation of growth and ovarian function (oocyte

ing similar feeding regimes and rearing densities across generations

growth and ovulation), energy homoeostasis and metabolism,

(Uusi-Heikkil€a et al., 2015), thus the potential for different environ-

immunological functions and various behaviours, such as aggression,

mental conditions to induce differences in gene expression was mini-

feeding and foraging behaviour, are synthesized or expressed in the

mized. See Appendix SI Materials and Methods: Zebrafish rearing

liver (Bjornsson et al., 2002). Differences in maturation status (i.e.,

conditions and harvesting for more details about fish rearing, feeding,

whether fish are immature, maturing or mature) likely lead to differ-

harvesting and breeding.

ences in expression level of maturation-associated genes given the metabolic and hormonal changes in the liver during maturation (Ng,

2.2 | RNA sampling

Tam, & Woo, 1986; Soengas, Barciela, & Aldegunde, 1995). The unusually highly expressed maturation-related genes could lead to

The liver transcriptomes were characterized for a total of 48 individ-

significantly lower coverage of other genes and bias the expression

uals which included four fish from each treatment replicate sampled

estimates of these genes, thus substantially complicating the detec-

in each of three generations: F2- (large- or random-harvested for

tion of differentially expressed genes of other biological processes.

two generations), F5- (large- or random-harvested for five genera-

Therefore, the liver sampling schedule did not follow the age (days

tions) and F11-generations (six generations of no harvesting following

post fertilization; DPF) of the experimental fish but rather aimed for

five generations of large or random harvesting). Before the first liver

sampling the fish at similar developmental stage. The fish matured,

€ UUSI-HEIKKILA

|

ET AL.

3957

and consequently were sampled, at different ages in each generation

To detect differentially expressed genes between the harvest

(F2 at age 116 DPF, F5 at age 69 DPF and F11 at age 97 DPF). How-

treatments, treatment replicates and among generations, a design

ever, there were no significant age differences between the random-

matrix including all explanatory variables (harvest treatment, genera-

and large-harvested fish within the F2-, F5- or F11-generations

tion, their interaction and treatment replicate) was created using

(Table S1). A detailed description of the RNA sequencing can be

contrasts between the explanatory variables. A linear model was

found in Appendix SI Materials and Methods: Sampling, sequencing

then fitted using the design matrix, and comparisons for main and

and processing the sequence data; RNASeq read alignment.

interaction effects between harvest treatments and among generations were extracted. All the explanatory variables were treated as fixed effects. When identifying differentially expressed genes, the p-

2.3 | Gene expression quantification

value (Benjamini-Hochberg adjusted for multiple testing) threshold

The effects of size-selective harvesting were assessed by comparing

was set to 0.05. In addition, differentially expressed genes were

the gene expression profiles of the large-harvested and random-har-

identified using the empirical false discovery rate (FDR). Differences

vested fish following five generations of harvesting (Figure 1). As all

in absolute fold change (i.e., nondirectional) of differentially

treatment replicates had by then likely experienced similar (inadver-

expressed genes between the harvest treatments and harvesting

tent) selection pressures associated with captive rearing, the ran-

periods were tested with a Wilcoxon signed-rank test as the differ-

dom-harvested treatment can be considered as a control with

ences were not normally distributed. To visualize differences in over-

respect to size-selective harvesting. To assess the effect of six gen-

all gene expression patterns across generations, we conducted a

erations of no harvesting following five generations of harvesting,

principal component analysis (PCA) on all scaled transcripts. Three

we compared the gene expression profiles of large- and random-har-

most important principal components were analysed separately for

vested treatments at the F11-generation (Figure 1). We additionally

each generation using a linear model with the principal component

studied the change in the gene expression profile of the large-har-

as a response variable and harvest treatment, treatment replicate,

vested fish by comparing the changes in gene expression from the

individual standard length and wet mass as explanatory variables. In

F2-generation to F11-generation to changes from the F2- to F5-gen-

the linear model,

eration (Figure 1). By doing so, our aim was to assess whether there

y
a þ b þ c þ d

were genes where the expression effects of size-selective harvesting remained following the six generation no-harvest period.

y was the principal component, a was the harvest treatment, b was

A second comparison aimed to assess the effects of captive

the treatment replicate, c was the standard length (mm), and d was

breeding. The random-harvested individuals were not intentionally

the wet mass (g). All explanatory variables were treated as fixed

selected for any obvious phenotypic trait, but it likely adapted to the

effects.

laboratory environment during the 11 generations of laboratory rear-

We used a hypergeometric test to study explicitly whether dif-

ing in addition to being affected by the genetic drift. Therefore,

ferentially expressed genes overlapped between certain comparisons

comparing differences in gene expression profiles within the ran-

(e.g., between the harvest treatments across generations) more than

dom-harvested treatment across generations (F2 vs F11) can provide

expected by chance. More information about the empirical FDR

an insight into the magnitude of gene expression changes induced

and testing for gene expression differences can be found in

by captive breeding (Figure 1). To obtain a more conservative esti-

Appendix SI Materials and Methods: Test for gene expression differ-

mate of the gene expression changes induced by captive breeding,

ences.

changes in gene expression profiles across generations in large- and

Functional

enrichment

analysis

was

conducted using

the

random-harvested fish were compared and genes where expression

Database for Annotation, Visualization and Integrated Discovery

changed in both harvest treatments were identified (Figure 1).

(DAVID v6.7) to identify gene ontology (GO) terms that were overrepresented in the genes expressed differently between har-

2.4 | Test for gene expression differences

vest treatments or between different generations of the same treatment (Huang, Sherman, & Lempicki, 2009). Although the

Read count files were analysed in R using the limma-voom (v2.9.8)

knowledge of overrepresented GO terms is helpful, it may only

linear model as implemented in the R/Bioconductor package (Anders

partially

& Huber, 2010). Read counts were normalized for the analysis using

explored genes (Antonov, Dietmann, Rodchenkov, & Mewes,

a “remove unwanted variation” (RUVg) normalization strategy (Risso,

2009). Therefore, the functional analysis was extended to study

Ngai, Speed, & Dudoit, 2014). The RUVg approach uses negative

protein–protein interaction (PPI) networks. We submitted a gene

control genes that are assumed to be not differentially expressed

list to Ingenuity Pathway Analysis (IPA) to study the most

among the treatments of interest (Risso et al., 2014). Therefore, to

important interaction networks of the differentially expressed

normalize expression levels, comparisons between harvest treat-

genes and identify their relevant upstream regulators. More

ments in all generations were combined, and the 5,000 least differ-

details of identification of GO terms and gene networks can

entially expressed genes were selected as the normalization gene

be found in Appendix SI Materials and Methods: Functional

set. Library sizes were normalized using the

EDGER

package.

resolve

the

molecular

mechanisms

enrichment and gene network analysis.

relevant

to

the

3958

€ UUSI-HEIKKILA

|

2.5 | Gene expression variance

ET AL.

individuals (resulting in excess of false positives). Although eQTLs for a certain gene can be found in any chromosome (trans-eQTLs), we

The variance in expression level for each gene was estimated across

focused on local (cis) eQTLs which map to the approximate location

the four individuals in each treatment replicate (interindividual vari-

(i.e., < 1M bp) of the gene.

ance) using RUVg-normalized and scaled read counts (Risso et al.,

We further compared the allele frequency changes (from F2-gen-

2014). The median of the gene-specific variances across all genes for

eration to F5-generation) between the large- and random-harvested

each treatment replicate was then calculated. Gene expression vari-

fish (i) across all 58,217 SNPs, (ii) in 14,500 gene-associated SNPs

ance within each generation and treatment replicate was boot-

(one association per SNP extracted by the

strapped 10,000 times with 1,000 random genes in each bootstrap.

lowest p-value) in differentially expressed genes and in genes that

MATRIXEQTL

based on the

After that, the change in the estimated median gene expression

were not differentially expressed (detected at the F5-generation),

variance from the F2-generation to F5-generation and from the

and finally (iii) in SNPs assigned as eQTLs and not assigned as

F5-generation to F11-generation was calculated. We further tested

eQTLs. The observed average allele frequency changes were com-

whether the difference in gene expression variance was significant

pared to a permuted allele frequency change distribution. The per-

within the harvesting (from F2 to F5) and the no-harvest periods (from

muted distribution consisted of 10,000 mean values, which were

F5 to F11) between the harvest treatments using a paired t-test.

sampled from a data set where we pooled the allele frequency changes of (i) all SNPs in random- and large-harvested fish, (ii) gene-

2.6 | SNP calling and annotation

associated SNPs occurring in differentially expressed genes and in genes that were not differentially expressed separately for both har-

All individuals were used for single nucleotide polymorphism (SNP)

vest treatments, and (iii) all SNPs assigned and not assigned as an

identification and calling. Only SNPs with a minimum allele frequency

eQTL separately for both harvest treatments.

of at least 0.1 were included. We used SAMtools (v.0.1.19) and the corresponding BCFtools to extract SNPs within a coverage range of 10 – 10,000, and filtered out SNPs within 10 bp distance from indels.

2.8 | Estimating the contribution of genetic drift

In addition to the default parameters in SAMtools, additional filtering

To determine the relative contribution of neutral (drift) and adaptive

steps were conducted: removal of potential indels, SNPs with an

(selection) processes during gene expression evolution, the variation

overall locus quality of lower than 999 (in SAMtools, a quality of 999

in gene expression level explained by the replicates (drift) and the

refers to a very good variant quality), and SNPs occurring outside

harvest treatment (selection) was estimated using MANOVA. For

chromosomes (mitochondrial and unassigned scaffold regions). SNPs

each generation, the amount of variation in gene expression

where at least 43 (of 48) individuals in each generation were success-

explained by the two models was compared: one with harvest treat-

fully genotyped with at least 159 coverage (i.e., at least 43 individuals

ment and treatment replicate as explanatory variables (replicate

had at least 15 reads covering a valid SNP position) were retained.

nested within harvest treatment) and one with only treatment repli-

Following these filtering steps, 58,217 SNPs remained. To visualize

cate as an explanatory variable. The input data for the MANOVA

differences in allele frequency patterns across generations, we con-

analysis were the three first principal components, which explained

ducted a PCA on SNPs retained in the analysis and with all individuals

most of the variation in expression count data between the harvest

having a valid genotype (21,181 SNPs; Fig. S1).

treatments and generations. All analyses were conducted in R (version 3.1.3; R Core Team

2.7 | Association between gene expression and SNP allele frequency changes To study the association between gene expression and SNP allele frequency changes in order to demonstrate a genetic basis for changes in gene expression levels, we first conducted an expression quantitative locus (eQTL) analysis. eQTLs are regions of the genome containing DNA sequence variants (e.g., SNPs) that influence the

2016) packages LME4, BASE STATS, LIMMA, RUVSEQ, EDGER and MATRIXEQTL.

3 | RESULTS 3.1 | Broad-scale differences in gene expression and reduction of gene expression variance after five generations of size-selective harvesting

expression level of one or more genes. To identify SNPs that are

Expression differences were detected between the large- and ran-

located close to genes and then link variation in expression values to

dom-harvested treatments in 509 transcripts (2.80% of all expressed

SNP genotypes, we used the MATRIXEQTL R package (Shabalin, 2012).

genes) at the F2-generation (Figure 2a). Three generations later, the

We performed linear regression of RUVg-normalized expression val-

number of differentially expressed genes had increased by more than

ues and SNPs with generation and treatment replicate as covariates

ninefold, with expression differences between large- and random-

(included to account for population stratification). SNPs with more

harvested fish in 4,310 genes (23.7%; Figure 2a). The absolute fold

than four genotypes missing were filtered out and gene expression

change of these differentially expressed genes was significantly

values were quantile normalized to fit a standardized normal distri-

higher among large-harvested (1.06  0.02; median  95% confi-

bution to avoid the overdispersion effect resulting from outlier

dence

interval)

than

random-harvested

fish

(0.92  0.01;

€ UUSI-HEIKKILA

|

ET AL.

3959

W = 1,0146,000, p < .001; Fig. S2A). Fold change was also signifi-

(0.125; Figure 5a). The average allele frequency change in gene-

cantly higher during the harvesting period (1.06  0.02; F2–F5) than

associated SNPs was generally higher in differentially expressed

during the no-harvest period (0.80  0.010; F5–F11) among the

genes than in genes that were not differentially expressed (p < .001;

large-harvested fish (W = 1,0567,000, p < .001; Fig. S2B). An indi-

Figure 5b). Although this was the case among both harvest treat-

vidual-level PCA revealed that the gene expression differences were

ments, the average SNP allele frequency change was higher in large-

explained by both selection and drift, as variation in principal compo-

harvested than in random-harvested fish in both differentially

nent 1 (PC 1) was explained by harvest treatment (F = 22.5,

expressed genes (large-harvested 0.174; random-harvested 0.152)

p < .001), in principal component 2 (PC 2) by harvest treatment

and in genes that were not differentially expressed (large-harvested

(F = 80.1, p < .001) and treatment replicate (F = 25.2, p < .001), and

0.159; random-harvested 0.128; Figure 5b). Overall, we identified

in principal component 3 (PC 3) by harvest treatment (F = 13.3,

221 eQTLs in 140 genes. eQTLs were present in 24 differentially

p = .004). A PCA figure visualizes that the gene expression profiles

expressed genes (of 4,300). The average allele frequency changes of

of the large-harvested and random-harvested individuals were not

SNPs not assigned as eQTLs were well within the permuted distribu-

markedly diverged at the F2-generation (Figure 3a). However, after

tions in both harvest treatments (Figure 5c) but unlike in random-

an additional three generations of size-selective harvesting (F5-gen-

harvested fish (p = .721), the average allele frequency change of

eration), the gene expression profiles of large- and random-harvested

eQTL SNPs in large-harvested fish (0.177) was significantly higher

fish had clearly diverged and formed distinct clusters (Figure 3b).

than in SNPs not assigned as eQTLs (0.140; p < .001; Figure 5c).

The differentially expressed genes between large- and randomharvested individuals at the F5-generation contributed to various major biological processes, such as translation and transcription (p < .001), lipid and steroid biosynthesis (p = .003) and energy meta-

3.2 | Gene expression response to a no-harvest period

bolism (p = .034; Table S2). The central genes of the three most sig-

After six harvest-free generations that followed five generations of

nificant gene networks were ELAVL1, EBP and MSMO1 (Fig. S3A-C;

size-selective or nonselective harvesting, differences in gene expres-

Table S3A). In addition, five putative upstream regulators were impli-

sion between large- and random-harvested fish had eroded to some

cated (Table S3B).

extent, but there still remained a large number of differentially

Reductions in gene expression variance between individuals

expressed genes (3,171 genes, 17.4%; Figure 2a). Variation in PC 1

within harvest treatments were evident in both size-selective (large-

was weakly explained by harvest treatment (F = 3.62, p = .078) and

harvested) and nonselective (random-harvested) harvesting regimes

treatment replicate (F = 3.77, p = .074). Variation in PC 2 was

during the harvesting period (Figure 4a; Fig. S4). The reduction in

explained by none of the variables, and in PC 3 by harvest treatment

variance was particularly clear in one of the replicates of both har-

(F = 45.6, p < .001), individual standard length (F = 14.7, p = .003)

vest treatments while in the other replicate the variance reduction

and wet mass (F = 11.9, p = .005). A consistent gene expression dif-

was less substantial (Figure 4a).

ference between the transcriptome profiles of harvest treatments in

The average allele frequency change during the harvesting period

the F11-generation was not evident in the two first principal compo-

in more than 58,000 SNPs was significantly (p < .001) higher among

nents (Figure 3c), but it was evident in the third principal component

large-harvested fish (0.140) compared to random-harvested fish

(Figure 3f). It appears that the effect of harvest treatment on gene

(a)

Between harvest treatments

F2 509

Large-harvested

Random-harvested

142

4 310

F11 826

3 171

111

(b)

Between generations

F5

F 2 vs. F 5

F 5 vs. F11

F2 vs. F11

3 112

2 564

3 007

414

801

924

2 083

5 764

4 978

F I G U R E 2 Number of differentially expressed genes between the harvest treatments and between generations within the harvest treatments. (a) Number of differentially expressed genes between large- and random-harvested fish in each generation (bold) and genes in common between generations (italics). (b) Number of differentially expressed genes between generations within a harvest treatment (bold) and genes that are shared between harvest treatments for the same generation comparison (italics) [Colour figure can be viewed at wileyonlinelibrary.com]

€ UUSI-HEIKKILA

| (b) 50 0

−50

0

50 0 −50

PC 2 (9.88%)

(c)

50

(a)

−100

−50

0

50

100

150

ET AL.

−50

3960

−100

−50

0

50

100

150

−100

−50

0

50

100

150

−50

0

50

100

150

PC 1 (16.9%) (e) 50 0

0

50

−50

−50

−50

0

PC 3 (8.47%)

(f)

50

(d)

−100

−50

0

50

100

150

−100

−50

0

50

100

150

−100

PC 2 (9.88%) F I G U R E 3 A principal component analysis (PCA) based on transcript abundances of all 18 192 expressed genes. PC 1 and PC 2 contrasted in (a) the F2-generation, (b) the F5-generation and (c) the F11-generation. PC 2 and PC 3 contrasted in (d) the F2-generation, (e) the F5generation, and (f) the F11-generation. Red diamonds represent large-harvested and grey circles random-harvested treatments. Different shades depict individuals from the two replicates within both harvest treatments [Colour figure can be viewed at wileyonlinelibrary.com]

(a)

Change in variance %

60

Large− harvested

Random− harvested

(b) 60

40

40

20

20

0

0

−20

−20

−40

−40

−60

−60

Large− harvested

Random− harvested

F I G U R E 4 Change in gene expression variance during the (a) harvesting (from F2 to F5) and (b) no-harvest period (from F5 to F11) in both harvest treatment replicates. There was no significant difference in gene expression variance between the harvest treatments during the harvesting period (t = 0.486, p = .627). During the no-harvest period, the difference in variance between the harvest treatments was significant (t = 618.1, p < .001). Error bars show the 95% confidence intervals. Different shades of red represent the two replicates within large-harvested fish and grey within random-harvested fish [Colour figure can be viewed at wileyonlinelibrary.com]

|

ET AL.

Frequency

(a)

0 500 1500 2500 3500

€ UUSI-HEIKKILA

Random−harvested

Frequency

(b)

0 500 1500 2500 3500

Frequency

0.13

Large−harvested

0.14

0.12

eQTL

0.11

0.12

0.13

non eQTL

0.15

0.16

SNPs in other differentially SNPs expressed genes

other SNPs

0.11

(c)

0.12

0 500 1500 2500 3500

0.11

0.14

0.15

0.16

0.17

0.14

0.18

SNPs in differentially expressed genes

0.17

non eQTL

0.13

3961

0.18

eQTL

0.15

0.16

0.17

0.18

Permuted mean allele frequency change

F I G U R E 5 The distribution of permuted SNP allele frequency changes compared to those observed (marked by arrows) in large- and random-harvested fish (from F2-generation to F5-generation). (a) Across all SNPs, (b) across gene-associated SNPs in differentially expressed genes (solid arrows) and in genes that are not differentially expressed, referred to as “other SNPs” (dashed arrows), and (c) across SNPs assigned as eQTLs (solid arrows) and SNPs not assigned as eQTLs (dashed arrows). Red bars and arrows represent large-harvested and grey bars and arrows random-harvested fish [Colour figure can be viewed at wileyonlinelibrary.com] expression differences was weaker and that the random processes

period were enriched for biological processes associated with protein

played a more important role following the no-harvest period (F11)

transport and localization (p = .015), and for genes in the insulin sig-

compared to following the harvesting period (F5).

nalling pathway (p = .045; Table S5).

Many of the genes differing between the harvest treatments were

Expression variance during the no-harvest period (F5–F11) had

enriched for biological processes related to RNA processing and meta-

differing patterns of change for the size-selected and nonselected

bolism (p < .001), protein catabolism (p = .004), ribosome biogenesis

harvest treatments compared to patterns during the harvesting per-

(p = .012) and nitrogen compound metabolism (p = .016; Table S4).

iod (F2–F5). Interestingly, gene expression variance decreased by

Central genes of the key networks were HNF1A, ESR1 and MDM2

25% in large-harvested fish after the cessation of harvesting, but

(Fig. S5A-C; Table S3C). Also, several significant upstream regulators

increased by 34% in random-harvested fish (Figure 4b; Fig. S4).

were implicated (Table S3D). Out of the more than 3,000 differentially expressed genes between the harvest treatments in the F11-generation only around a quarter (826 genes) were in common with genes

3.3 | Genetic drift

differentially expressed between the large- and random-harvested

In the F5-generation, the MANOVA model with the harvest treat-

lines at the F5-generation (p < .001; Figure 2a).

ment and the treatment replicate (full model) as explanatory vari-

A more conservative approach to study the effect of no-harvest

ables explained substantially more of the variation (54%) in the three

period on a size-selectively exploited population is to compare the

principal components as opposed to a model with only a treatment

gene expression profiles of the large-harvested fish across genera-

replicate as an explanatory variable (reduced model; 28%; Table S6),

tions (i.e., from F2 to F11). In that comparison, we identified 3,007

suggesting that differences in gene expression were better explained

differentially expressed genes (Figure 2b) out of which 1,159 (38.5%)

by the effects of both selection and genetic drift rather than by

were the same genes that were affected by size selection during the

genetic drift alone.

harvesting period (p < .001). These over 1,000 genes were largely different compared to the differentially expressed genes shared between the harvesting and no-harvest period in the random-harvested fish (75 genes in common; p = .006). The genes that were

3.4 | The effects of captive rearing on gene expression

differentially expressed in large-harvested fish during the harvesting

To identify gene expression changes potentially induced by captive

period and remained differentially expressed after the no-harvest

rearing, we compared the expression patterns of the random-

3962

€ UUSI-HEIKKILA

|

ET AL.

harvested fish at the F2- and F11-generations, and identified 4,978

expression between wild and farmed Atlantic salmon (Salmo salar)

differentially expressed genes (27.4%). These genes contributed to

subjected to four to seven generations of selection aimed at increas-

biological processes such as DNA metabolic processes (p = .001),

ing growth rate revealed a difference in less than 2% of genes (Debes,

(p = .029)

catabolism

Normandeau, Fraser, Bernatchez, & Hutchings, 2012; Roberge, Einum,

(p = .029; Table S7). The central genes of the three most significant

Guderley, & Bernatchez, 2006). One potentially important difference

networks

S6A-C;

between our study and the domestication studies is the fact that

Table S3E). More than 20 upstream regulators were implicated

those studies focused on gene expression changes in the brain, rather

(Table S3F). Out of the almost 5,000 differentially expressed genes

than the liver, and this could have limited the number of genes differ-

between the F2- and F11-generations in random-harvested fish, 924

ing significantly in their expression (Jobling, Hollox, Hurkes, Kivisild, &

were in common with those differentially expressed between the F2-

Tyler-Smith, 2014). A notable exception is a recent study on fish

and F11-generations in large-harvested fish (p < .001; Figure 2b).

domestication that demonstrated large expression differences similar

Approximately 3% of these genes were associated with insulin sig-

to what we found between wild and hatchery-reared steelhead trout

nalling pathway (p = .006).

(Oncorhynchus mykiss) after only one generation of hatchery rearing

oxidative

phosphorylation were

CHAF1B,

TRIP13

and

and

protein

CYB5R2

(Fig.

(Christie et al., 2016). The authors identified 723 differentially expressed genes (4.6%) between wild and first-generation hatchery

4 | DISCUSSION

offspring reared in common environment.

Our experimental approach to understand the molecular conse-

ments have documented divergence in gene expression among differ-

In addition to domestication studies, other artificial selection experiquences of harvest-induced evolution revealed that broad-scale

ently selected lines. For instance, foxes (Vulpes vulpes) selected for

changes in gene expression and SNP allele frequencies arose after

aggressive or tame behaviour in a long-term breeding experiment dif-

just five generations of size-selective harvesting. This suggests that

fered in their gene expression in 335 genes (Kukekova et al., 2011), and

size-selective harvesting can rapidly induce large gene expression

the fraction of genes that diverged under directional selection on body

changes in exploited populations, and the extent of such changes

weight in chicken (Gallus domesticus) was shown to be on average 13%

cannot be explained by genetic drift alone. A no-harvest period

(Resnyk et al., 2015). Selection for improved residual feed intake

resulted in the expression divergence of a largely different set of

resulted in gene expression differences in less than 4% of all the

genes indicating that although selection and drift also occur in the

expressed genes in beef cattle (Bos taurus; Weber et al., 2016) and in

absence of size-selective harvesting, the targets of these processes

pigs (Sus scrofa; Vincent et al., 2015), and in 41 genes in chickens with

are different and “recovery” may be unlikely to occur, even under a

divergent feed efficiency (Yi et al., 2015). While aggression, body

fishing moratorium. Exploited populations, particularly in freshwaters,

weight and feeding behaviour are also polygenic traits similar to body

are often supplemented with hatchery-reared individuals. Therefore,

length on which selection was operating on in our study, the larger gene

we also studied the rate and magnitude of gene expression changes

expression differences we observed due to size-selective harvesting

induced by captive breeding. We demonstrated substantial gene

may be partly related to the differences in tissues sampled (e.g., brain,

expression changes in random-harvested fish suggesting that

muscle or intestine vs liver), differences in gene expression analyses and

although not selected for any obvious phenotypic trait, captive

therefore in sensitivity to detect differentially expressed genes (microar-

breeding and rearing alone might induce large unintentional genetic

ray vs RNASeq) or to the fact that many of these studies were con-

changes in fish populations.

ducted on domesticated animals which may have already lost genetic diversity. Nevertheless, our results demonstrate that size-selective har-

4.1 | Gene expression differences after five generations of size-selective harvesting

vesting (and not drift alone) can rapidly induce gene expression changes on a broader scale than reported in many of the previous domestication studies or selection experiments conducted on vertebrates. Further

After five generations of size-selective harvesting, the number of dif-

research is required to identify the specific mechanisms causing the

ferentially expressed genes between large- and random-harvested fish

gene expression differences. However, given the common-garden envi-

encompassed approximately 24% of all expressed genes investigated,

ronment, which minimized environmental differences between the har-

although relatively few differences (<3%) were present early in the

vest treatments, along with substantially higher average allele frequency

selection experiment at the F2-generation. During the harvesting per-

changes in differentially expressed genes compared to genes that were

iod, the magnitude of change in differentially expressed genes was

not differentially expressed in both harvest treatments (Figure 5b) and

also higher among large-harvested than random-harvested fish

in eQTL SNPs in large-harvested fish compared to SNPs not assigned as

(Fig. S2A). Such broad-scale gene expression differences are unusual

eQTLs (Figure 5c), it is likely that at least some of the expression differ-

compared to many earlier studies investigating gene expression differ-

ences we report have a genetic basis.

ences under domestication selection. For example, a comparison of

Earlier research investigating the effects of harvesting on life his-

brain gene expression levels between domesticated vs wild dogs, pigs

tory traits in the same selection lines as used in the current study

or rabbits revealed expression differences in less than 1% of all

showed that while individuals from large- and random-harvested

expressed genes (Albert et al., 2012). Similarly, differences in gene

treatments had similar prematuration growth rates, large-harvested

€ UUSI-HEIKKILA

|

ET AL.

3963

individuals reached a significantly lower adult body size than ran-

17%, but the differences in gene expression profiles were still sub-

dom-harvested individuals (Uusi-Heikkil€a et al., 2015), exhibited a

stantial after this period and expression differences remained in

lower condition factor (Fig. S7) and a lower age-specific maturation

genes that are linked to growth-related processes at the physiologi-

probability (Uusi-Heikkil€a et al., 2015). Although the individuals anal-

cal and behavioural level. Differentially expressed genes detected

ysed in the present study were immature and there were no differ-

between the harvest treatments after the no-harvest period were

ences in body size between the large- and random-harvested fish at

enriched in biological processes that are important in energy release

the time of sampling (Table S1), gene expression profiles can never-

and cell growth. Among the most highly expressed regulators,

theless potentially differ at earlier developmental phases between

hypocretin receptor 2 (HCRTR2) is known to be associated with

zebrafish lines selected for differing growth patterns (Amaral & John-

appetite and feeding behaviour (De Lecea et al., 1998).

ston, 2012), for instance when fish begin to allocate energy to matu-

The genes that were differentially expressed between the har-

differentially

vest treatments at the F11-generation were generally different to

expressed genes between large- and random-harvested fish were

those that differed at the F5-generation: only a quarter (26%) of the

enriched in biological processes that are related to energy allocation

differentially expressed genes between the large-harvested and ran-

and growth, albeit at a very general level. Thus, these differences

dom-harvested stocks were identical between the harvesting and

could relate to differences in adult body size which have previously

the no-harvest period. Thus, the proportion of genes that were

been found to evolve in response to size-selective mortality (Uusi-

affected by size-selective harvesting had decreased while a new set

Heikkil€a et al., 2015).

of genes, likely affected by another selection pressure and/or drift,

ration

and

postmaturation

growth.

Indeed,

the

Protein turnover is a fundamental biological process related to

were identified. The random-harvested fish were under selection for

somatic growth (Reeds, 1988). We detected a significant proportion of

captive rearing during the entire experiment, while the large-har-

the differentially expressed transcripts between large- and random-

vested fish experienced an intensive size-selection pressure for five

harvested fish to be enriched for processes associated with protein

generations (in addition to selection for captive rearing) and after

synthesis and metabolism, such as transcription and translation.

that, selection solely for captive rearing that likely began to favour

Another important process differing significantly between the harvest

very different characteristics than size selection (e.g., large, fecund

treatments was oxidative phosphorylation, which is the metabolic

and aggressive individuals; Roberge et al., 2006; Devlin, Sakhrani,

pathway forming ATP that stores and supplies energy, for instance for

Tymchuk, Rise, & Goh, 2009). This might suggest that the large-har-

protein and lipid metabolism (Alberts et al., 2002). Harvest treatments

vested fish were set on a different evolutionary trajectory. In addi-

also differed in expression of genes associated with steroid synthesis.

tion, the magnitude of the expression of differentially expressed

Steroid hormones are associated with many major biological functions

genes was significantly lower during the no-harvest period than dur-

in fish, including feeding behaviour, aggression, stress and oocyte mat-

ing the harvesting period among the large-harvested fish (Fig. S2B).

uration (Nagahama & Yamashita, 2008; Wingfield et al., 1998). Ster-

Indeed, the lack, or only very slow recovery, of phenotypes after a

oids are also known to promote growth in fish by enhancing anabolic

phase of intensive harvesting has been repeatedly demonstrated in

processes (Matty, 1986). Central genes in key gene networks were

individual-based eco-genetic models (Dunlop, Eikeset, & Stenseth,

associated with hypoplasia, growth failure and biosynthesis of choles-

2015; Enberg et al., 2009; Marty et al., 2015).

terol, which serves as a precursor for the biosynthesis of steroid hor-

Even after six generations of no harvesting (i.e., at F11), signifi-

mones. Furthermore, genes that were differentially expressed and up-

cant differences in expression levels remained in one-third of the

regulated during the harvesting period in large-harvested fish (Fig. S8)

genes observed to have differential expression patterns after the five

were tightly associated with maturation (Table S8). It could be specu-

generation harvesting period. These approximately 1,000 genes were

lated that differences between harvest treatments found earlier in

enriched for biological functions such as protein transport and local-

maximum body size, exploration behaviour (Uusi-Heikkil€a et al., 2015)

ization, and insulin signalling pathway. One could argue that these

and condition factor between the harvest treatments are consistent

genes might have been under selection to captive rearing but this

with the observed differences in expression of genes that can be

seems unlikely because out of these, only 75 genes were in common

broadly associated with feeding, circadian rhythms, somatic growth

with the genes that were differentially expressed among random-

and maturation. However, it is good to keep in mind that while our

harvested fish in the equivalent generation-level comparison.

data represent interesting signals and gene ontology enrichments, fol-

Although the differentially expressed genes were mostly different

low-up studies are required to establish more direct links between the

between the harvesting and no-harvest period, a small group of

gene expression differences observed here, and the phenotypic differ-

reproduction- and maturation-related genes were up-regulated

ences observed previously (Uusi-Heikkil€a et al., 2015).

among large-harvested fish during the harvesting period and downregulated during the no-harvest period (Fig. S8). Large-harvested fish

4.2 | Gene expression differences following a no-harvest period

have been shown to invest a relatively high amount of energy into reproduction (Uusi-Heikkil€a et al., 2015), which could at least partly explain the up-regulation of these genes during the harvesting per-

The number of differentially expressed genes between the harvest

iod. Selection pressure for reproductive investment likely relaxed

treatments decreased during the no-harvest period from 24% to

after the harvesting was halted. Admittedly, age could have

3964

€ UUSI-HEIKKILA

|

ET AL.

confounded the between-generation comparisons as the fish were

€ tterer, 2015). If decanalization was stronger (Chen, Nolte, & Schlo

sampled in each generation at different age (in days). However, the

among random-harvested fish, this could have led to increased gene

expression of maturation-related genes might not have been con-

expression variance compared to the large-harvested fish. However,

founded by the sampling design as fish were always sampled at the

this remains speculative as quantifying decanalization among harvest

same developmental stage (i.e., they were all immature). Although

treatments was beyond the scope of this study.

expecting the populations to reverse back to the early-harvest state

Gene expression variation represents a source of variability that

might not be entirely realistic, at least not in the current laboratory

can improve fitness in varying environments and under stressful con-

setting where the fish were under selection for captive rearing, both

ditions or varying selection pressures (Papakostas et al., 2014;

above mentioned approaches (i.e., comparison between harvest

Whitehead & Crawford, 2006). Although sometimes considered as

treatments and generations) generally lead to the same conclusion: a

costly noise (Wang & Zhang, 2011), interindividual gene expression

component of the effects of size-selective harvesting still remained

variation has been shown to contribute substantially to physiological

despite six generations of no harvesting.

performance among individuals, thus it can be biologically relevant (Li, Liu, Kim, Min, & Zhang, 2010). Therefore, loss of gene expression

4.3 | The effect of harvesting and a no-harvest period on gene expression variance During the harvesting period, gene expression variance decreased in both harvest treatments. This could be due to reduced population

variation could be detrimental for exploited fish populations because it can reduce adaptive potential and hamper their recovery.

4.4 | Gene expression response to captive rearing

sizes but it is possible that also random-harvested fish were under

Today, many declining wild populations, especially freshwater fishes,

selection pressure due to high harvest rate or intrinsic fecundity

are supplemented with captive-reared individuals and this can have

selection (Uusi-Heikkil€a et al., 2015). However, the reduction in gene

many detrimental ecological and genetic effects on existing popula-

expression variance was not entirely consistent within the harvest

tions (Laikre, Schwartz, Waples, & Ryman, 2010; Lorenzen et al.,

treatments (i.e., one of the treatment replicate in both harvest treat-

2012). The genetic and ecological concerns include direct genetic

ments showed less reduced variance) possibly due to the limited

effects caused by introgression or hybridization, genetic changes in

number of biological replicates.

hatchery stocks brought about by selection and drift, and lowered sur-

We show not only that gene expression variance was reduced

vival and reproductive success of hatchery-reared individuals (Araki

after five generations of size-selective and nonselective harvesting,

et al., 2007; Christie et al., 2014). Introgression of genetic material

but it continued to decrease among large-harvested fish during the

from captive-reared fish that are maladapted in the wild may indeed

no-harvest period which mimicked a harvest moratorium. However,

cause negative fitness effects in the wild populations, alter the gene

among random-harvested fish gene expression variance increased

pools of local stocks and negatively affect population productivity

during the no-harvest period. This might suggest that the combined

(Chilcote, Goodson, & Falcy, 2011). Although it is known that hatchery

effects of selection and drift were stronger in the large-harvested

fish may have lower relative reproductive success (Araki et al., 2007;

treatment and/or the response to captive rearing was relaxed in the

Christie et al., 2014) and survival rate (Lorenzen et al., 2012) than

random-harvested treatment, but not in the large-harvested line

their wild counterparts, the systematic functional genetic effects of

after the cessation of harvesting. Thus, it is possible that large- and

captive rearing remain unclear (but see Christie et al., 2016).

random-harvested fish responded differently to captive rearing. Fur-

We showed that nine generations of rearing and breeding in cap-

ther, it has been suggested that adaptive (life history) evolution can

tivity affected the expression of a large number of genes. Over 27%

be very rapid during early phases of selection but it may also cease

of all expressed genes investigated were differentially expressed

rapidly (Reznick, Shaw, Rodd, & Shaw, 1997). It is plausible that the

between the F2- and F11-generations in random-harvested fish.

large-harvested fish experienced two types of selection: first size-

However, as the fish at the F2-, F5- and F11-generations were sam-

selection favouring small body size and then selection to captive

pled at different ages due to the fact that age at maturity was evolv-

breeding potentially favouring large body size (Devlin et al., 2009;

ing during the experiment one cannot rule out the possibility that

Heath, Heath, Bryden, Johnson, & Fox, 2003). These two opposite

age differences could have also contributed to between-generation

selection forces could have decreased gene expression variance in

comparisons within harvest treatments. Most earlier studies that

large-harvested fish compared to the random-harvested fish. Finally,

have shown gene expression changes induced by captive rearing

based on the expression (Figure 3a-c) and SNP PCA figures (Fig. S1),

have focused on a relatively small set of candidate genes (Debes

it appears that during the early selection, the individuals from the

et al., 2012; Roberge et al., 2006), and few studies thus far have

different harvest treatments were already genetically different, but

examined this at a genomewide level (Christie et al., 2016). Our

their gene expression profiles were not markedly different. Thus, at

results demonstrate that captive rearing affected biological processes

the F2-generation they seemed to exhibit strong canalization which

that can be broadly related to growth, such as energy production,

broke down at the F5-generation potentially due to selection, and

protein catabolism, and fatty acid metabolism, despite random-har-

gene expression remained decanalized at the F11-generation.

vested fish were not directly selected for growth or large (or small)

Decanalization could uncover hidden gene expression variance

body size. Our findings of the functional genetic effects of captive

€ UUSI-HEIKKILA

|

ET AL.

3965

rearing are broadly in line with salmon domestication studies (Chris-

hybrids (Christie et al., 2014), which has been shown to reduce pop-

tie et al., 2016; Debes et al., 2012; Roberge et al., 2006; Tymchuk,

ulation productivity in salmonids (Chilcote et al., 2011). Therefore, it

Sakhrani, & Devlin, 2009), although many of the earlier studies com-

can be questioned whether hatchery-reared fish should be used

pared wild fish and fish under intense growth selection and there-

intensively for enhancing wild populations that are still naturally

fore cannot be directly compared to our results.

reproducing as hybridization might have far-reaching consequences

The key gene networks of the differentially expressed genes

in terms of adaptive capacity and genomic integrity of wild popula-

were associated with RNA repair, post-transcriptional modification

tions without necessarily increasing fisheries yield (Lorenzen et al.,

and cholesterol biosynthesis (Table S3E). Some of the significant

2012).

upstream regulators are known to be involved in reproduction through mediating oestrogen and progesterone production, steroidogenesis and cell proliferation (Table S3F). Thus, adaptation to captive

4.5 | Conclusions

environment may have not only affected growth-related but also

We have shown that both size-selective harvesting and captive rear-

maturation-related processes. In fact, we showed that maturation-

ing can induce rapid and substantial changes in gene expression pro-

and reproduction-related genes were up-regulated among the ran-

files in experimentally exploited fish populations. Gene expression

dom-harvested fish (from F5 to F11), and we have shown earlier that

profiles did not fully converge after the cessation of harvesting. Our

the random-harvested fish have a higher age-specific maturation

results thus suggest that the evolutionary response to size-selective

probability than large-harvested fish (Uusi-Heikkil€a et al., 2015). This

harvesting can be broad, rapid and potentially difficult to reverse. This

is in agreement with other studies showing that fish held in captivity

can be undesirable for fishing because phenotypes (and genotypes)

tend to mature later than their wild conspecifics (Debes & Hutch-

favoured most by natural selection and by fishers (e.g., the large fish)

ings, 2014). Thus, it is possible that in captivity, traits under selection

are removed. It can also be harmful from the evolutionary perspective

are related to body condition (body fat content), growth, and poten-

because of reduced gene expression variance, thus potentially reduced

tially maturation. Although we do not assess fitness consequences

adaptive potential of the exploited populations. Reduction in gene

directly in this study, earlier research has demonstrated the negative

expression variance could thus be a factor potentially contributing to

fitness consequences for wild populations of captivity-induced modi-

the lack of genetic (and phenotypic) recovery of size-selectively

fication of reproductive traits (reviewed in Christie et al., 2014).

exploited fish populations in the wild. Hence, our results reinforce the

A second approach for identifying genes potentially affected by

recommendation of applying evolutionary principles to management

captive rearing was to compare the gene expression profile diver-

and promote management that maintains large and diverse breeding

gence of large-harvested fish between the F2- and F11-generations

populations to foster the full range of phenotypes and genotypes that

to that of random-harvested fish between the same generations.

natural selection can act upon (Schindler et al., 2010). Inadvertent

The comparison of gene expression profiles of F2- and F11-genera-

selection due to captive rearing alone also resulted in large changes in

tions in large-harvested fish revealed 3,000 differentially expressed

gene expression profiles, suggesting that the use of hatchery-reared

genes and in random-harvested fish almost 5,000 genes (Figure 2b).

fish for supplementing wild populations might affect the adaptive

Out of these two sets of differentially expressed genes, only 924

capacity and genomic integrity of wild populations.

genes were in common between the two harvest treatments. These approximately 900 genes potentially include those predominantly affected by selection for captive rearing, although it is possible that

ACKNOWLEDGEMENTS

some of the remaining almost 4,000 genes are as well, but have

We thank Karena Kuntze, Asja Vogt, Marcus Ebert, Sylvia Werner,

been differently affected by the alternative selection regimes. A sig-

Sarah Becker, Yvonne Klaar, Theresa Arlt, Julie Menard and David

nificant proportion of the 900 differentially expressed genes in com-

Lewis for fish husbandry, care taking and data collection; Sanna Pau-

mon between the harvesting treatments were enriched within insulin

sio for assistance in RNA extraction; Henrik Zwadlo for technical

signalling pathway. Insulin is tightly associated with growth, thus

assistance and Christian Wolter and Thomas Meinelt for designing

indirectly with gonad development. The insulin (and insulin-like

and helping with the selection experiment. This research used com-

growth factor) signalling pathway has also been shown to play a

puting resources of the CSC-IT Center for Science, Espoo, Finland.

major role in the control of longevity in vertebrates (van Heemst,

Funding for this study was received through the AXA Research

2010) and in the balance of the functioning of the immune system.

Fund, ICES, and Kone Foundation to SUH, Adaptfish and BTypes

Traits mediated by the insulin signalling pathway, such as growth,

grants by the Leibniz Community to RA, and Academy of Finland to

maturation and disease resistance, are known to be important traits

CRP. Finally, we thank three reviewers for excellent suggestions and

for artificial breeding programmes (Gjedrem & Thodesen, 2005).

constructive feedback.

These results add an important functional genetic element to the genetic concerns of stocking: the divergent gene expression profiles of captive-reared individuals reported here and elsewhere (Christie

DATA ACCESSIBILITY

et al., 2016) likely contribute to the reduced fitness of captive-reared

The raw data from each library are available at the NCBI Sequence

individuals in the wild (Araki et al., 2007), as well as captive-wild

Read Archive (SRA) under the Accession no. SRP105243.

3966

|

AUTHOR CONTRIBUTION S.U.H., C.R.P., E.L. and R.A. designed the study; S.U.H. collected the data, S.U.H., T.S. and C.R.P. analyzed the data; S.U.H., C.R.P., T.S., E.L. and R.A. wrote the manuscript.

REFERENCES Albert, F. W., Somel, M., Carneiro, M., Aximu-Petri, A., Halbwax, M., Thalmann, O., . . . P€ a€abo, S. (2012). A comparison of brain gene expression levels in domesticated and wild animals. Plos Genetics, 8, e1002962. Alberts, B., Johnsson, A., Lewis, J., Raff, M., Roberts, K., & Walter, P. (2002). Molecular biology of the cell. New York, USA: Garland Science. Allendorf, F. W., England, P. R., Luikart, G., Ritchie, P. A., & Ryman, N. (2008). Genetic effects of harvest on wild animal populations. Trends in Ecology and Evolution, 23, 327–337. Allendorf, F. W., & Hard, J. J. (2009). Human-induced evolution caused by unnatural selection through harvest of wild animals. Proceedings of the National Academy of Sciences, 106, 9987–9994. Amaral, I. P. G., & Johnston, I. A. (2012). Experimental selection for body size at age modifies early life-history traits and muscle gene expression in adult zebrafish. The Journal of Experimental Biology, 215, 3895–3904. Anders, S., & Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biology, 11, R106. Andersen, K. H., & Brander, K. (2009). Expected rate of fisheries-induced evolution is slow. Proceedings of the National Academy of Sciences, 106, 11657–11660. Antonov, A. V., Dietmann, S., Rodchenkov, I., & Mewes, H. W. (2009). PPI spider: A tool for the interpretation of proteomics data in the context of protein-protein interaction networks. Proteomics, 9, 2740– 2749. Araki, H., Cooper, B., & Blouin, M. S. (2007). Genetic effects of captive breeding cause a rapid, cumulative fitness decline in the wild. Science, 318, 100–103.  s, J., Klefoth, T., Monk, C. T., Arlinghaus, R., Laskowski, K. L., Alo € der, A. (2017). Passive gear-induced timidity Nakayama, S., & Schro syndrome in wild fish populations and its potential ecological and managerial implications. Fish and Fisheries, 18, 360–373. € rnsson, B. T., Johansson, V., Benedet, S., Einarsdottir, I. E., Hildahl, J., Bjo € nsson, E. (2002). Growth hormone endocrinology Agustsson, T., & Jo of salmonids: Regulatory mechanisms and mode of action. Fish Physiology and Biochemistry, 27, 227–242. Browman, H. I., Law, R., & Marshall, C. T. (2008). The role of fisheriesinduced evolution. Science, 320, 47. € tterer, C. (2015). Temperature stress mediates Chen, J., Nolte, V., & Schlo decanalization and dominance of gene expression in Drosophila melanogaster. PLoS Genetics, 11, e1004883. Chilcote, M. W., Goodson, K. W., & Falcy, M. R. (2011). Reduced recruitment performance in natural populations of anadromous salmonids associated with hatchery-reared fish. Canadian Journal of Fisheries and Aquatic Sciences, 68, 511–522. Christie, M. R., Ford, M. J., & Blouin, M. S. (2014). On the reproductive success of early-generation hatchery fish in the wild. Evolutionary Applications, 7, 883–896. Christie, M. R., Marine, M. L., Fox, S. E., French, R. A., & Blouin, M. S. (2016). A single generation of domestication heritably alters the expression of hundreds of genes. Nature Communications, 7, 10676. Conover, D. O., Munch, S. B., & Arnott, S. A. (2009). Reversal of evolutionary downsizing caused by selective harvest of large fish. Proceedings of the Royal Society B: Biological Sciences, 276, 2015–2020. Darimont, C. T., Carlson, S. M., Kinnison, M. T., Paquet, P. C., Reimchen, T. E., & Wilmers, C. C. (2009). Human predators outpace other agents of trait change in the wild. Proceedings of the National Academy of Sciences, 106, 952–954.

€ UUSI-HEIKKILA

ET AL.

De Lecea, L., Kilduff, T. S., Peyron, C., Gao, X. B., Foye, P. E., Danielson, P. E., . . . Sutcliffe, J. G. (1998). The hypocretins: Hypothalamus-specific peptides with neuroexcitatory activity. Proceedings of the National Academy of Sciences, 95, 322–327. Debes, P. V., & Hutchings, J. A. (2014). Effects of domestication on parr maturity, growth, and vulnerability to predation in Atlantic salmon. Canadian Journal of Fisheries and Aquatic Sciences, 71, 1371–1384. Debes, P. V., Normandeau, E., Fraser, D. J., Bernatchez, L., & Hutchings, J. A. (2012). Differences in transcription levels among wild, domesticated, and hybrid Atlantic salmon (Salmo salar) from two environments. Molecular Ecology, 21, 2574–2587. Devlin, R. H., Sakhrani, D., Tymchuk, W. E., Rise, M. L., & Goh, B. (2009). Domestication and growth hormone transgenesis cause similar changes in gene expression in coho salmon (Oncorhynchus kisutch). Proceedings of the National Academy of Sciences, 106, 3047–3052. Devlin, R. H., Sundstrom, L. F., & Muir, W. M. (2006). Interface of biotechnology and ecology for environmental risk assessments of transgenic fish. Trends in Biotechnology, 24, 89–97. Dunlop, E. S., Eikeset, A. M., & Stenseth, N. C. (2015). From genes to populations: How fisheries-induced evolution alters stock productivity. Ecological Applications, 25, 1860–1868. Enberg, K., Jorgensen, C., Dunlop, E. S., Heino, M., & Dieckmann, U. (2009). Implications of fisheries-induced evolution for stock rebuilding and recovery. Evolutionary Applications, 2, 394–414. Fraser, H. B., Moses, A. M., & Schadt, E. E. (2010). Evidence for widespread adaptive evolution of gene expression in budding yeast. Proceedings of the National Academy of Sciences, 107, 2977–2982. Gjedrem, T., & Thodesen, J. (2005). Selection. In T. Gjedrem (Ed.), Selection and breeding programs in aquaculture (pp. 89–110). Dordrecht: Springer. Heath, D. D., Heath, J. W., Bryden, C. A., Johnson, R. M., & Fox, C. W. (2003). Rapid evolution of egg size in captive salmon. Science, 299, 1738–1740. van Heemst, D. (2010). Insulin, IGF-1 and longevity. Aging and Disease, 1, 147–157. Heino, M., Dıaz Pauli, B., & Dieckmann, U. (2015). Fisheries-induced evolution. Annual Review of Ecology, Evolution, and Systematics, 46, 461–480. Hilborn, R., & Minte-Vera, C. V. (2008). Fisheries-induced changes in growth rates in marine fisheries: Are they significant? Bulletin of Marine Science, 83, 95–105. Huang, D. W., Sherman, B. T., & Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protocols, 4, 44–57. Huminiecki, L., & Wolfe, K. H. (2004). Divergence in spatial gene expression profiles following species-specific gene duplications in human and mouse. Genome Research, 14, 1870–1879. Jobling, M., Hollox, E., Hurkes, M., Kivisild, T., & Tyler-Smith, C. (2014). Human evolutionary genetics. New York: Garland Science. Jones, F. C., Grabherr, M. G., Chan, Y. F., Russell, P., Mauceli, E., Johnson, J., . . . Kingsley, D. M. (2012). The genomic basis of adaptive evolution in threespine sticklebacks. Nature, 484, 55–61. Jørgensen, C., Enberg, K., Dunlop, E. S., Arlinghaus, R., Boukal, D. S., Brander, K., . . . Rijnsdorp, A. D. (2007). Managing evolving fish stocks. Science, 318, 1247–1248. Kukekova, A., Johnson, J. L., Teiling, C., Li, L., Oskina, I. N., Kharlamova, A. V., . . . Acland, G. M. (2011). Sequence comparison of prefrontal cortical brain transcriptome from a tame and an aggressive silver fox (Vulpes vulpes). BMC Genomics, 12, 482. Kuparinen, A., & Festa-Bianchet, M. (2017). Harvest-induced evolution: Insights from aquatic and terrestrial systems. Philosophical Transactions of the Royal Society B, 372, 20160036. Laikre, L., Schwartz, M. K., Waples, R. S., Ryman, N., & Ge M Working Grp (2010). Compromising genetic diversity in the wild: Unmonitored large-scale release of plants and animals. Trends in Ecology and Evolution, 25, 520–529. Lewin, W., Arlinghaus, R., & Mehner, T. (2006). Documented and potential biological impacts of recreational fishing: Insights for

€ UUSI-HEIKKILA

|

ET AL.

management and conservation. Reviews in Fisheries Science, 14, 305–367. Li, J. J., Liu, Y., Kim, T., Min, R. Q., & Zhang, Z. L. (2010). Gene expression variability within and between human populations and implications towards disease susceptibility. PLoS Computational Biology, 6, e1000910. Lorenzen, K., Beveridge, M. C. M., & Mangel, M. (2012). Cultured fish: Integrative biology and management of domestication and interactions with wild fish. Biological Reviews, 87, 639–660. Manceau, M., Domingues, V. S., Mallarino, R., & Hoekstra, H. E. (2011). The developmental role of agouti in color pattern evolution. Science, 331, 1062–1065. Marty, L., Dieckmann, U., & Ernande, B. (2015). Fisheries-induced neutral and adaptive evolution in exploited fish populations and consequences for their adaptive potential. Evolutionary Applications, 8, 47–63. Matty, A. J. (1986). Nutrition, hormones and growth. Fish Physiology and Biochemistry, 2, 141–150. McKinney, G. J., Hale, M. C., Goetz, G., Gribskov, M., Thrower, F. P., & Nichols, K. M. (2015). Ontogenetic changes in embryonic and brain gene expression in progeny produced from migratory and resident Oncorhynchus mykiss. Molecular Ecology, 24, 1792–1809. Miller, C. T., Beleza, S., Pollen, A. A., Schluter, D., Kittles, R. A., Shriver, M. D., & Kingsley, D. M. (2007). cis-regulatory changes in kit ligand expression and parallel evolution of pigmentation in sticklebacks and humans. Cell, 131, 1179–1189. Nagahama, Y., & Yamashita, M. (2008). Regulation of oocyte maturation in fish. Development, Growth & Differentiation, 50, S195–S219. Neubauer, P., Jensen, O. P., Hutchings, J. A., & Baum, J. K. (2013). Resilience and recovery of overexploited marine populations. Science, 340, 347–349. Ng, T. B., Tam, P. L. L., & Woo, N. Y. S. (1986). Sexual maturation in the black seabream, Mylio macrocephalus Teleostei, Sparidae: Changes in pituitary gonadotropes, hepatocytes and related biochemical constituents in liver and serum. Cell Tissue Research, 245, 207–213. Palumbi, S. R. (2001). Evolution - Humans as the world’s greatest evolutionary force. Science, 293, 1786–1790. Papakostas, S., Vasemagi, A., Himberg, M., & Primmer, C. R. (2014). Proteome variance differences within populations of European whitefish (Coregonus lavaretus) originating from contrasting salinity environments. Journal of Proteomics, 105, 144–150. Pinsky, M. L., & Palumbi, S. R. (2014). Meta-analysis reveals lower genetic diversity in overfished populations. Molecular Ecology, 23, 29–39. R Core Team (2016). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. http:// www.R-project.org/ Reeds, P. J. (1988). Regulation of protein metabolism. In J. F. Quirke, & H. Schmid (Eds.), Control and regulation of animal growth (pp. 25–44). Pudoc: Wageningen. Resnyk, C. W., Chen, C., Huang, H., Wu, C. H., Simon, J., Le Bihan-Duval, E., . . . Cogburn, L. A. (2015). RNA-Seq analysis of abdominal fat in genetically fat and lean chickens highlights a divergence in expression of genes controlling adiposity, hemostasis, and lipid metabolism. PLoS ONE, 10, e0139549. Reznick, D. N., Shaw, F. H., Rodd, F. H., & Shaw, R. G. (1997). Evaluation of the rate of evolution in natural populations of guppies (Poecilia reticulata). Science, 275, 1934–1937. Risso, D., Ngai, J., Speed, T. P., & Dudoit, S. (2014). Normalization of RNA-seq data using factor analysis of control genes or samples. Nature Biotechnology, 32, 896–902. Roberge, C., Einum, S., Guderley, H., & Bernatchez, L. (2006). Rapid parallel evolutionary changes of gene transcription profiles in farmed Atlantic salmon. Molecular Ecology, 15, 9–20. Ryman, N., & Laikre, L. (1991). Effects of supportive breeding on the genetically effective population size. Conservation Biology, 5, 325–329. Schindler, D. E., Hilborn, R., Chasco, B., Boatright, C. P., Quinn, T. P., Rogers, L. A., & Webster, M. S. (2010). Population diversity and the portfolio effect in an exploited species. Nature, 465, 609–612.

3967

Shabalin, A. A. (2012). Matrix eQTL: Ultra fast eQTL analysis via large matrix operations. Bioinformatics, 28, 1353–1358. Shapiro, M. D., Marks, M. E., Peichel, C. L., Blackman, B. K., Nereng, K.  nsson, B., . . . Kingsley, D. M. (2004). Genetic and developmental S., Jo basis of evolutionary pelvic reduction in threespine sticklebacks. Nature, 428, 717–723. Soengas, J. L., Barciela, P., & Aldegunde, M. (1995). Variations in carbohydrate metabolism during gonad maturation in female turbot (Scophthalmus maximus). Marine Biology, 123, 11–18. Stearns, S. C. (1992). The evolution of life histories. Oxford: Oxford University Press. Sullivan, A. P., Bird, D. W., & Perry, G. H. (2017). Human behavior as long-term ecological driver of non-human evolution. Nature Ecology & Evolution, 1, 0065. Swain, D. P., Sinclair, A. F., & Hanson, J. M. (2007). Evolutionary response to size-selective mortality in an exploited fish population. Proceedings of the Royal Society, 274, 1015–1022. Tymchuk, W., Sakhrani, D., & Devlin, R. (2009). Domestication causes large-scale effects on gene expression in rainbow trout: Analysis of muscle, liver and brain transcriptomes. General and Comparative Endocrinology, 164, 175–183. Uusi-Heikkil€a, S., Whiteley, A. R., Kuparinen, A., Matsumura, S., Venturelli, P. A., Wolter, C., . . . Arlinghaus, R. (2015). The evolutionary legacy of size-selective harvesting extends from genes to populations. Evolutionary Applications, 8, 597–620. feu, C., Gilbert, H., & Lefaucheur, Vincent, A., Louveau, I., Gondret, F., Tre L. (2015). Divergent selection for residual feed intake affects the transcriptomics and proteomic profiles of pig skeletal muscle. Journal of Animal Science, 93, 2745–2758. Wang, Z., & Zhang, J. (2011). Impact of gene expression noise on organismal fitness and the efficacy of natural selection. Proceedings of the National Academy of Sciences, 108, E67–E76. Weber, K. L., Welly, B. T., Van Eenennaam, A. L., Young, A. E., PortoNeto, L. R., Reverter, A., & Rincon, G. (2016). Identification of gene networks for residual feed intake in Angus cattle using genomic prediction and RNA-Seq. PLoS ONE, 11, e0152274. Whitehead, A., & Crawford, D. L. (2006). Variation within and among species in gene expression: Raw material for evolution. Molecular Ecology, 15, 1197–1211. Wingfield, J. C., Maney, D. L., Breuner, C. W., Jacobs, J. D., Lynn, S., Ramenofsky, M., & Richardson, R. D. (1998). Ecological bases of hormone-behavior interactions: The “emergency life history stage”. American Zoologist, 38, 191–206. Yi, G., Yuan, J., Bi, H., Yan, W., Yang, N., & Qu, L. (2015). In-depth duodenal transcriptome survey in chickens with divergent feed efficiency using RNA-Seq. PLoS ONE, 10, e0136765.

SUPPORTING INFORMATION Additional Supporting Information may be found online in the supporting information tab for this article.

How to cite this article: Uusi-Heikkil€a S, S€avilammi T, Leder E, Arlinghaus R, Primmer CR. Rapid, broad-scale gene expression evolution in experimentally harvested fish populations. Mol Ecol. 2017;26:3954–3967. https://doi.org/ 10.1111/mec.14179