Received: 5 April 2017

|

Revised: 9 June 2017

|

Accepted: 21 June 2017

DOI: 10.1111/mec.14234

ORIGINAL ARTICLE

Gene expression stasis and plasticity following migration into a foreign environment Brian K. Lohman1

| William E. Stutz2 | Daniel I. Bolnick1

1

Department of Integrative Biology, University of Texas at Austin, Austin, TX, USA

Abstract Selection against migrants is key to maintaining genetic differences between popu-

2

Office of Institutional Research, Western Michigan University, Kalamazoo, MI, USA

lations linked by dispersal. However, migrants may mitigate fitness costs by proactively choosing among available habitats, or by phenotypic plasticity. We

Correspondence Brian K. Lohman, Department of Integrative Biology, University of Texas at Austin, Austin, TX, USA. Email: [email protected]

previously reported that a reciprocal transplant of lake and stream stickleback (Gasterosteus aculeatus) found little support for divergent selection. Here, we revisit that experiment to test whether phenotypic plasticity in gene expression may have helped migrants adjust to unfamiliar habitats. We measured gene expression

Funding information Howard Hughes Medical Institute (DIB); David and Lucille Packard Foundation (DIB)

profiles in stickleback via TagSeq and tested whether migrants between lake and stream habitats exhibited a plastic response to their new environment that allowed them to converge on the expression profile of adapted natives. We report extensive gene expression differences between genetically divergent lake and stream stickleback, despite gene flow. But for many genes, expression was highly plastic. Fish transplanted into the adjoining habitat partially converged on the expression profile typical of natives from their new habitat. This suggests that expression plasticity may soften the impact of migration. Nonetheless, lake and stream fish differed in survival rates and parasite infection rates in our study, implying that expression plasticity is not fast or extensive enough to fully homogenize fish performance. KEYWORDS

convergence, gene expression, migration, plasticity, stickleback

1 | INTRODUCTION

mating success (Hereford, 2009). This selection against migrants is a key to maintaining genetic differences between populations linked

What happens to an organism when it moves into a new habitat?

by dispersal (Lenormand, 2002; Nosil et al., 2005). Yet, migrants may

Populations in disparate environments commonly exchange migrants.

evade selection in two ways. First, they can proactively choose

These migrant individuals are exposed to unfamiliar abiotic condi-

among available habitats to avoid environments to which they are

tions and biotic communities to which their phenotypes may be

mismatched (Edelaar & Bolnick, 2012; Edelaar, Siepielski, & Clobert,

poorly suited (Hereford, 2009; Kawecki & Ebert, 2004; Lenormand,

2008). Or, migrants may plastically alter one or more phenotypic

2002). Migrants can be maladapted to their new habitat because

traits to acclimate to a new habitat (Davidson, Jennions, & Nicotra,

they inherited alleles that were selectively favoured in their native


pez-Maury, 2011; Ghalambor, McKay, Carroll, & Reznick, 2007; Lo

range but are untested by selection in their new habitat (Nosil,

Marguerat, & B€ahler, 2008).

Vines, & Funk, 2005). Or, migrants’ traits may have been shaped,

Plasticity is most frequently measured as change in phenotype in

during ontogeny, by their native environment (Davis & Stamps,

response to an environmental change. Reciprocal transplants or com-

2004; Stamps & Davis, 2006). Either way, migrants’ poor fit to their

mon garden experiments have been successful at partitioning the rel-

new habitat may frequently result in reduced survival, fecundity or

ative contributions of heritability vs. plasticity for a myriad of

Molecular Ecology. 2017;26:4657–4670.

wileyonlinelibrary.com/journal/mec

© 2017 John Wiley & Sons Ltd

|

4657

4658

|

LOHMAN

ET AL.

ecologically relevant traits (Conover & Present, 1990; Pfennig, 1992;

each habitat, (ii) differences in gene expression associated with being

Schlichting & Pigliucci, 1998; West-Eberhard, 2003). A limitation of

moved from one habitat to another and (iii) convergence in expres-

this literature, however, is a tendency to focus on readily measured

sion profiles between native and transplanted individuals.

phenotypic traits (e.g., morphology and size), which may not be the

The threespine stickleback fish (Gasterosteaus aculeatus) offers an

most crucial traits for migrants’ performance and fitness, and which

opportunity to study plasticity of both phenotypes and gene expres-

may not be representative of plasticity for other more subtle yet

sion. Across Vancouver Island, British Columbia, there are many

important traits (e.g., immunity, physiological homeostasis).

replicate pairs of lake and stream stickleback (Hendry, Taylor, &

Gene expression profiling offers a much broader approach to

McPhail, 2002; Thompson, Taylor, & McPhail, 1997). These parap-

assay the response of an individual to both abiotic and biotic stres-

atric lake and stream populations are typically genetically and pheno-

sors. Phenomics is often limited to the study of a few morphological

typically divergent (Eizaguirre et al., 2011; Feulner et al., 2015;

traits that are identified a priori. In contrast, transcriptomics casts a

Reusch, Wegner, & Kalbe, 2001; Roesti, Hendry, Salzburger, & Ber-

broader net across many possibly relevant traits, although the choice

ner, 2012; Weber, Bradburd et al., 2017). These phenotypic differ-

of tissue to obtain RNA still places some constraints on the trait

ences persist to some degree in constant laboratory settings

space being studied. Transcriptomics thus allows for more agnostic

indicating there are heritable differences (Oke et al., 2016) [for other

discovery of relevant traits or pathways. However, because of the

common garden studies of phenotypic plasticity see (Berner et al.,

substantial cost of transcriptomic analyses, there are few studies of

2011; Jiang, Peichel, Ling, & Bolnick, 2017; Kalbe & Kurtz, 2006;

transcriptome-wide plasticity in natural settings, and most of these

Scharsack, Kalbe, Harrod, & Rauch, 2007)]. Adjoining lake and stream

have very limited biological replication (Todd, Black, & Gemmell,

environments differ in both abiotic and biotic conditions including

2016). Other studies have achieved higher replication (and thus

flow regime, oxygen concentration, habitat structure, resource avail-

power) by testing for plasticity of just a few candidate genes. For

ability, prey composition and parasite communities (Berner et al.,

example, Stutz, Schmerer, Coates, and Bolnick (2015) showed that

2009; Kaeuffer, Bolnick, Hendry, & Peichel, 2012; Lenz, Eizaguirre,

stickleback fish transplanted between lakes converged strongly to

Rotter, Kalbe, & Milinski, 2013; Stuart et al., 2017). The magnitude

resemble the immune gene expression profile (for seven candidate

and direction of environmental differences between a lake and its

genes) of natives of their new environment, indicating strong plastic-

outlet stream effectively predict the direction of phenotypic differ-

ity (Stutz et al., 2015). But, is this plasticity particular to immune

entiation between lake and stream resident stickleback (Stuart et al.,

genes, or is it representative of gene expression across the tran-

2017). The implication, invoked by many studies of lake and stream

scriptome? We expect that the whole transcriptome may respond in

stickleback (summarized in Weber, Bradburd et al., 2017), is that

one of four general patterns: (i) a large stress response, (ii) a lack of

environmental differences drive divergent selection on lake and

response, suggestive of a tolerance strategy, or (iii) a plastic

stream stickleback.

response that may ameliorate selection (note that response 1 or 2

To test for this inferred selection, multiple studies have trans-

may be either adaptive or nonadaptive, depending on subsequent

planted lake and stream stickleback into their native and neighbour-

changes in fitness). (iv) Plastic responses that allow immigrants to

ing habitats, measuring whether residents systematically outperform

converge on the resident phenotype may be adaptive if residents

immigrants in a variety of measures (survival, growth, infection; Bol-

are locally adapted.

nick, 2004; Bolnick & Stutz, 2017; Hanson, Moore, Taylor, Barrett, &

Here, we describe how the stickleback transcriptome responds

Hendry, 2016; Hendry et al., 2002; Kaufmann, Lenz, Kalbe, Milinski,

to an unfamiliar environment. We recently conducted a reciprocal

& Eizaguirre, 2017; Moser, Frey, & Berner, 2016; Scharsack, Kalbe

transplant of threespine stickleback, moved between adjacent but

et al., 2007). However, these experiments yielded surprisingly incon-

ecologically very different lake and stream habitats (Bolnick & Stutz,

sistent evidence for divergent selection (summarized in extended

2017; Stuart et al., 2017). The lake and stream populations are

data of (Bolnick & Stutz, 2017)). Why is divergent selection rarely

divergent with respect to morphology (Berner, Grandchamp, & Hen-

observed (but see (Hendry et al., 2002; Kaufmann et al., 2017),

dry, 2009; Oke et al., 2016; Stuart et al., 2017), genomic SNPs

despite evidence of phenotypic divergence? Several recent studies

(Weber, Bradburd, Stuart, Stutz, & Bolnick, 2017) and immune gene

discuss the possibility of habitat choice helping to maintain lake–

allele frequencies (MHC IIb, Stutz and Bolnick, 2017). These differ-

stream differences (Berner & Thibert-Plante, 2015; Bolnick et al.,

ences strongly suggest that there is divergent selection, yet our

2009; Jiang, Torrance, Peichel, & Bolnick, 2015; Weber, Bradburd

reciprocal transplant experiment yielded little signal of local adapta-

et al., 2017). Another possibility is that plasticity mitigates selection

tion (Bolnick & Stutz, 2017). We hypothesized that plasticity may

against migrants. Here, we use a reciprocal transplant experiment

help migrants ameliorate the ill effects of dispersal into a foreign

that found negligible support for divergent selection, to also test for

neighbouring habitat (Oke et al., 2016). If so, we would expect trans-

plasticity. We measured both physical traits (e.g., change in mass

planted fish’s transcriptomes to converge on the expression patterns

over a given period or the value of an ecologically relevant trait) and

of the native population. Here, we present a test of this prediction

gene expression profiles via TagSeq (Lohman, Weber, & Bolnick,

by analysing sticklebacks’ transcriptomic response to transplantation

2016) for a large number of transplanted individuals. Using this data,

between adjoining lake and stream habitats. Specifically, we tested

we tested whether migrants’ gene expression shifts to more closely

for (i) baseline differences in gene expression between natives of

resemble expression by the native population in their new habitat,

LOHMAN

|

ET AL.

4659

suggesting a role for expression-mediated phenotypic plasticity by

stream) contained three fish. Half the cages received a 1:2 ratio of

migrants.

lake:stream fish, the other half of the cages received a 2:1 ratio.

There is ample evidence for phenotypic plasticity in ecologically

Thus, a total of 240 fish were transplanted, with 60 in each of the

relevant traits in stickleback. For example, previous experiments

four treatments detailed below. Within each cage, the three fish

reared stickleback from different habitats in a common garden set-

were uniquely marked with dorsal spine clips to facilitate identifica-

ting (laboratory aquaria), and fed them alternative diets to test for

tion. After 8 weeks, Bolnick and Stutz recaptured the caged stickle-

plasticity in feeding morphology (Day & McPhail, 1996; Svanb€ack &

back. As a control for the effect of caging, Bolnick and Stutz also

Schluter, 2012). These studies measured body shape, gill raker and

collected wild uncaged fish from both lake and stream at the conclu-

gape traits that are both readily measured and clearly relevant to

sion of the experiment, from habitat immediately adjoining the

foraging. Life-history traits also show plasticity in stickleback, includ-

cages. Hereafter, here we refer to uncaged fish as the “wild” group,

ing breeding size, clutch size, egg size and relative clutch mass (Baker

all fish recovered from cages are “transplanted.” Within the trans-

& Foster, 2002). Finally, prior studies have examined plasticity in

planted fish, we distinguish between “natives” (same origin and desti-

gene expression (Gibbons, Metzger, Healy, & Schulte, 2017; Leder

nation habitats) and “immigrants” (different origin and destination).

et al., 2014; Robertson, Bradley, & MacColl, 2016; Wang et al.,

At the conclusion of the field experiment, Bolnick and Stutz euth-

2014). One such study focused on expression of two candidate

anized the collected fish with an overdose of MS-222. Fish were

genes for osmoregulation and salinity tolerance (McCairns & Ber-

weighed, measured and dissected to remove head kidneys (“prone-

natchez, 2010). A larger, whole-transcriptome approach suggested

phros”) which were stored in RNAlater (Ambion) for subsequent RNA

that the invasion of freshwater and thermal tolerance drove the evo-

extraction and expression analysis. Head kidney was chosen because

lution of gene expression plasticity (Morris et al., 2014). However,

it is the major site of hematopoiesis and the site of an immune

while these studies of gene expression plasticity have sought to

response (Fischer, Koppang, & Nakanishi, 2013; Fischer et al., 2006;

answer how the transcriptome may respond to a novel environment,

Scharsack, Kalbe, Derner, & Millinski, 2004; Scharsack, Koch, & Ham-

they have been carried out in the laboratory and do not account for

merschmidt, 2007). After dissection, specimens were preserved in

the diverse stressors of the wild. We therefore tested whether

ethanol for later dissection to enumerate parasites by complete dis-

migrants between lake and stream habitats indeed exhibit a strong

section and examination under dissecting microscope (including body

plastic response to their new environment that allows them to con-

cavity, external surface and all organs). Morphological features were

verge on the gene expression profile of the native population. To

measured with digital callipers (pelvic width is the width of the pelvic

the extent that native gene expression is honed by previous natural

girdle, body depth is the distance from the base of the first dorsal

selection, it is reasonable to suspect that such convergence reflects

spine and the anterior point of the pelvic girdle, gape width is the

adaptive transcriptional plasticity.

distance between mouth corners). Because of its importance in defence against parasites, Bolnick and Stutz (2017) sequenced MHC

2 | METHODS 2.1 | Sample acquisition

IIb exon 2 from all caged fish, using DNA from prerelease spine clips. MHC IIb was sequenced and data were analysed as described in Stutz and Bolnick (2014). A previous clinal survey of stickleback from this lake and stream revealed population differences in MHC IIb

We analyse data from a reciprocal transplant experiment using stick-

allele frequencies and significant associations between MHC alleles

leback from Roberts Lake and Stream (Vancouver Island, British

and the prevalence of particular parasites (Stutz & Bolnick, 2017).

Columbia, Canada), whose fitness effects were previously reported

The lake–stream transplant experiment revealed that transplanted

by Bolnick and Stutz (2017). Prior studies have documented differ-

foreign fish accumulated higher parasite infections than native fish,

ences between these populations, with respect to genotype and

but parasites apparently exploited native MHC genotypes (Bolnick &

phenotype (Weber, Bradburd et al., 2017; Stutz & Bolnick, 2017;

Stutz, 2017; Stutz & Bolnick, 2014). Here, we use the MHC data to

Berner et al., 2009, among others). Wild-caught stickleback were

test for correlations between genotype and transcriptome profile.

trapped, weighed, measured for length, individually marked with unique spine clips and then placed in cylindrical wire cages. Lake cages and stream cages both received a total of 60 lake fish and 60 stream fish. Each cage was ~1.6 m in diameter, placed in 1 m deep

2.2 | RNAseq library preparation, sequencing and bioinformatics

water and sealed to the substrate to prevent escape. Cages were

Following total RNA extraction (Ambion AM1830), we built 96 indi-

made of wire mesh that allowed free flow of water and movement

vidual TagSeq libraries according to Lohman et al. (2016). We

of prey items. In the lake, cages were situated along the shoreline

selected fish to construct a fully balanced design with 16 individuals

approximately 150 m from the outlet stream. In the stream, cages

in each treatment. We selected individuals that had been housed in

were placed 150 m downstream from the lake. In an effort to

the same cage, as available. We prioritized cage over sex ratio result-

reduce the influence of gene flow on stream genotypes, stream fish

ing in more males than female (48 vs. 39 in the final data set). How-

for the experiment were gathered from 1.5 km downstream of the

ever, the sex ratios with transplanted fish are nearly even: foreign

lake outlet. Each of the 80 enclosures (40 in the lake and 40 in the

transplants:

15/14

and

native

transplants:

16/14

(males/

4660

|

LOHMAN

ET AL.

females). TagSeq libraries were sequenced on the HiSeq 2500 with

size of 30 genes. We used dynamic tree cut and merged modules

1x100V4 chemistry at the Genome Sequencing and Analysis Facility

with greater than 80% similarity, producing a total of 11 modules.

at the University of Texas at Austin, generating an average of ~5 M

We plot the FDR-corrected Pearson correlation coefficient between

raw reads per sample. This read depth is appropriate for TagSeq (be-

module eigengenes and trait values. We have assumed that lake and

cause sequencing effort is targeted at the 30 end of the mRNA) and

stream fish have similar gene expression networks and combined

has been shown to be successful (Dixon et al., 2015; Kenkel & Matz,

both to generate a single coexpression network. However, merging

2016; Lohman et al., 2016; Meyer, Aglyamova, & Matz, 2011; Meyer

lake and stream fish may confound correlations between modules

et al., 2009; Wright, Aglyamova, Meyer, & Matz, 2015).

and traits with an effect of origin. To ensure that this hidden vari-

Raw reads were processed (removal of adapter contamination,

able problem did not obscure our results, we also generated popula-

poly-A and PCR duplicates followed by quality filtering, n—20)

tion-specific networks (thereby eliminating the hidden variable) with

according to the iRNAseq pipeline (Dixon et al., 2015; Lohman et al.,

identical parameters and recalculated correlations between modules

2016; Meyer et al., 2011), producing a total of 19,556 genes. The

and traits.

stickleback genome contains 20,787 genes, so we conclude our read depth was sufficient, especially considering that we sampled a single tissue at a single time point. We further filtered these genes by removing all genes for which the mean among all samples was less

2.3.2 | What is the effect of being transplanted into a novel environment?

than one, resulting in 9,748 genes for further analysis. Due to a

We tested for changes in gene expression of transplanted (caged)

machine error during the HiSeq run, BaseSpace was unable to con-

fish as function of origin habitat, destination habitat and the inter-

vert cycle 35 to a base call, and thus base 35 is N in every read. We

action between origin and destination. Using our estimated gene

adjusted for this by adding the –n option to all calls to fastx_clipper

network, we calculated the FDR-corrected Pearson correlation

in the iRNAseq pipeline. Mapping with Bowtie2 should not be influ-

coefficient between module eigengenes and traits unique to the

enced by this error (~53.3% alignment rate (min = 14.07%,

transplant design (treatment, origin, destination, delta mass and

max = 61.3%), postquality filtering, adaptor trimming and poly-A

delta length). A main effect of origin indicates stable gene expres-

removal). GO enrichment was conducted according to Wright et al.

sion differences between native lake vs. native stream fish. These

data-

expression differences can be stable because they are heritable, or

base and following previously described procedures (Dixon et al.,

because they are environmentally induced only early in ontogeny

2015; Lohman, Steinel, Weber, & Bolnick, In review).

but remain canalized in adults, which we used for this experiment.

(2015) using transcriptome annotation built with the

UNIPROTKB

A main effect of destination indicates genes that respond plasti-

2.3 | Statistical analysis

cally to recently experienced environments. An interaction between origin and destination would indicate ecotype differences in how

We analysed gene expression using a series of linear models in

they respond to a given environment. Such interactions could

DESeq2 (Love, Huber, & Anders, 2014), limma (Ritchie et al., 2015)

entail G*E effects on expression, but we point out they could also

(R Development Core Team 2007). All raw p values were

arise from ecotype differences in the extent of canalization of

and base

R

multiple test corrected (10% FDR, Benjamini–Hochberg). We sought

early plasticity.

to estimate three effects.

2.3.1 | What are the differences between wild fish from Roberts Lake and Stream?

2.3.3 | How well do immigrants converge on the expression profile of natives? We conducted a PCA of expression of all genes in all fish and then

We tested for differences in gene expression between wild

selected only transplanted fish and used the leading 57 PC axes (ex-

(uncaged) fish from Roberts Lake vs. Roberts Stream by modelling

plaining 90% of the total variance) for subsequent linear discriminant

gene count as a function of origin (lake or stream). We tested for

analysis. The original expression matrix has too much collinearity for

GO enrichment within the main effect of origin with a Mann–Whit-

LDA. Dropping higher-order PCA axes reduces this collinearity,

ney U via GO_MWU (Dixon et al., 2015). We used weighted gene

enabling LDA. This approach is sometimes called DAPC (Jombart,

coexpression network analysis (WGCNA; Langfelder & Horvath,

Devillard, & Balloux, 2010; Kenkel & Matz, 2016). We plotted these

2008) to estimate correlations between suites of coexpressed genes

results in LDA space, adding vectors connecting each ecotype’s

and traits, including morphology, parasite burdens and genotypes

expression at home to the same ecotype’s expression in the foreign

(e.g., MHC allelic diversity). WGCNA is an unbiased, data-driven

habitat. These vectors represent the magnitude and direction of

method to cluster groups of genes with similar expression patterns.

expression plasticity along DAPC axes. Convergence in expression

We removed batch effects and normalized counts using limma

would result in an angle of 180° between the vectors for lake and

(Ritchie et al., 2015) before starting WGCNA. We followed the tuto-

stream ecotypes. Moreover, we compared the lengths of these vec-

rial of Langfelder and Horvath (2008), and constructed a signed net-

tors to evaluate whether lake and stream ecotypes are equally plas-

work with a soft thresholding power of 6 and a minimum module

tic.

LOHMAN

|

ET AL.

4661

Lastly, if plasticity effectively recreates lake–stream differences,

correlated with infection by nematodes (r = .22, p < .04, Figure 3b)

then we would expect that genes that are more highly expressed in

and genes in the red module are correlated to infection by any spe-

lake natives would also be more highly expressed in fish transplanted

cies of Proteocephalus (Figure 3c, r = .28, p < .01). In addition to

into the lake. This can be tested by measuring the correlation, across

population-independent correlations between modules and traits, we

genes, between the origin effect sizes and destination effect sizes

also observe correlations which are population specific. There are a

estimated in analysis (2) above. Adaptive plasticity generating con-

larger number of significant correlations between infection by para-

vergence on the native expression profile should result in a positive

sites and modules in this population-specific setting, suggesting that

correlation.

lake and stream fish may respond differently to different kinds of parasites (see Supporting information).

3 | RESULTS 3.1 | What are the differences between wild fish from Roberts Lake and Stream?

3.2 | What is the effect if being transplanted into a novel environment? Focusing next on transplanted (caged) stickleback, we observed sig-

Our linear model revealed that 647 genes were differentially

nificant effects of both origin and destination for many genes (507

expressed between wild Roberts Lake and Roberts Stream stickle-

and 111, respectively when p < .05, see Supporting information for

back (Wald, p < .1 after 10% FDR, or 306 when p < .05). GO analy-

full list, after 10% FDR). Here, the effect of origin represents geno-

sis showed that these genes are enriched both for a variety of

type effects that persisted after transplantation (because the effect

categories including (but not limited to) genes regulating macrophage

of transplantation is averaged). Approximately 94% of the genes with

differentiation (biological processes, Mann–Whitney U, p < .05 after

significant (p < .1, after 10% FDR) origin effect in transplanted fish

10% FDR correction, Figure 1), and genes involved in the MHC

were also significantly different between wild fish ecotypes. This

Class II protein complex (cellular components, Mann–Whitney U,

overlap of origin effects in caged and wild fish suggests that stickle-

p < .1 after 10% FDR). Both of these GO groups have known func-

back exhibit realistic lake–stream expression differences when placed

tions in parasite defence and have been previously implicated in the

in lake or stream cages.

response to selection and parasite prevalence in stickleback (Bolnick

Destination effects represent plasticity that was independent

& Stutz, 2017; Eizaguirre, Lenz, Kalbe, & Milinski, 2012; Lohman

of genotype (genotype effects are averaged in our model). Most

et al., In review). Previous studies revealed that Roberts Lake and

notably, this list of genes includes hsp90 (lower in fish trans-

Stream stickleback populations harbour significantly different para-

planted into the stream, Wald, p  .001 after 10% FDR), a stress

site communities (Bolnick & Stutz, 2017), with corresponding differ-

response protein which has been studied in many different ani-

ences in MHC Class II allele frequencies (Stutz & Bolnick, 2017).

mals (Queitsch, Sangster, & Lindquist, 2002; Rutherford & Lind-

In addition to gene-by-gene linear modelling, we also tested for

quist, 1998). In addition, stat1 (lower in fish transplanted into the

correlations between modules of coexpressed genes and various

stream, Wald, p < .09 after 10% FDR) was also significantly differ-

traits, including morphology, infection by parasites and MHC Class

ent between fish transplanted in alternate environments. This tran-

IIb genotype (Figure 2). We found morphology to be correlated with

scription factor has a rich history of study for its critical role in

many different modules, each with modest correlation but highly sig-

multiple signalling cascades throughout the immune system (Mur-

nificant p values. It is noteworthy that all modules except the tur-

phy, 2011).

quoise module have a negative correlation with morphology

We found only 10 genes whose expression depended on the

(coexpressed gene modules are given arbitrary colour names). There

interaction of origin and destination (Wald, p < .1 after 10% FDR, or

are correlations between MHC allele number and several modules,

4 when p < .05, Figure 4; Table S1, Fig. S1). Such interactions can

including greenyellow, blue, magenta and pink. Interestingly, MHC

loosely be interpreted as genotype by environment interactions (e.g.,

allele number and two measures of parasite diversity have equal

genetic differences in plasticity), with the caveat that we are study-

strength but opposite sign in their correlation to the greenyellow

ing wild-caught fish, so we cannot infer heritable differences with

module (Figure 3a). This is consistent with previous experimental

certainty. Of these 10 genes, two candidates are possibly involved in

and theoretical data that animals with more diverse MHC genotypes

defence against parasites: cyp24a1, a cytochrome p450 variant

should have fewer parasites (Wegner, Kalbe, Kurtz, Reusch, & Milin-

(Annalora et al., 2010), and dhx58, an antiviral gene about which lit-

ski, 2003). Finally, we also considered linear discriminant axes of

tle is known (Leavy, 2012). In both cases, lake natives have higher

MHC II genotypes from a prior analysis of these same fish. We find

expression in the lake than do stream fish, but decreased expression

that LDA1 and LDA3 of MHC II are correlated with turquoise and

when moved into the stream. Stream fish have higher expression in

purple modules. The turquoise expression module is also correlated

their native habitat, but only higher than foreign lake fish for dhx58.

with fish origin (r = .36, p  .001), so these correlations are likely a

Furthermore, it is noteworthy that for all 10 interaction genes, the

result of differences in MHC genotype between the two ecotypes.

genes are more highly expressed in lake than in stream fish (all in

Infection by several functional groups of parasites is significantly cor-

lake cages). And, for all 10 genes, the lake natives decrease expres-

related with particular modules. For instance, the purple module is

sion when moved into the stream (Figure 4).

4662

|

LOHMAN

16/102 peptidyl−lysine modification

37/212 peptidyl−amino acid modification 19/119 histone modification

p < 0.01

67/257 membrane region

p < 0.05

3/6 microvillus membrane

p < 0.1

114/470 plasma membrane part

43/239 chromatin organization 33/175 cell division

2/5 ripoptosome 20/70 synapse

3/8 regulation of macrophage differentiation

76/328 cell junction

21/58 positive regulation of neuron differentiation

1/7 MHC class II protein complex

7/11 homophilic cell adhesion via plasma membrane adhesion molecules

52/240 extracellular space

22/76 sensory perception

9/44 extracellular matrix component

51/212 system process

54/289 extracellular region

36/137 neurological system process

22/89 extracellular matrix

8/25 neuromuscular process

40/393 ribonucleoprotein complex

10/26 pigment metabolic process

11/98 ribosome

7/14 heme metabolic process

6/63 ribosomal subunit

9/13 tetrapyrrole biosynthetic process

17/66 cytosolic part

12/22 tetrapyrrole metabolic process

6/9 hemoglobin complex

6/9 heme biosynthetic process

29/228 transferase complex

7/19 pigment biosynthetic process

30/217 cellular amide metabolic process 26/190 peptide metabolic process 80/495 organonitrogen compound metabolic process

ET AL.

42/330 catalytic complex 38/300 nucleoplasm part 17/155 nuclear body 3/12 DNA repair complex

57/308 organonitrogen compound biosynthetic process

F I G U R E 1 GO analysis results on Lake vs. Stream wild fish. Blue terms are underexpressed while red terms are overexpressed relative to the lake baseline. p values are Mann–Whitney U. Dendrograms indicate similarity of GO groups. Left group is from the biological processes cluster while the right group is from cellular components

We used DESeq2 to estimate caging effects by comparing wild

into foreign or native environment) and origin are both correlated

fish to natives within each environment. We found a moderate num-

with the turquoise module. In contrast, destination is only weakly

ber of differentially expressed genes; 35 genes were differentially

correlated to the pink and magenta modules. The red module has a

expressed between wild lake fish and caged lake natives. Somewhat

negative correlation to origin and a positive correlation to change in

more genes (79) were differentially expressed between wild stream

length over the course of the experiment. Change in mass is corre-

fish and stream natives (all Wald, p < .1 after 10% FDR, or 19 and 52

lated with both the magenta and purple modules. Interestingly, there

when p < .05, respectively. See Supporting information for full list).

is no overlap between change in mass and change in length. This dif-

There are very few notable differences due to caging in lake geno-

ference suggests a change in condition within individuals (Figure 5).

types. Lake natives have higher expression of cyp24a1 than wild lake fish (log2 fold change = 3.8, p = .065 after 10% FDR correction). Lake transplants also have higher expression of ebf4, an early B-cell factor (log2 fold change = 4.3, p = .049 after 10% FDR correction)

3.3 | How well do immigrants converge on the expression profile of wild controls?

than wild lake fish. In contrast, when we make the same contrast but

We tested for convergence between natives and immigrants in the

in stream genotypes, almost all differentially expressed genes (76 of

entire expression profile. Within a bivariate discriminant function

79 passing p < .1 after 10% FDR correction) exhibit a pattern of

space, we found that LDA1 separates fish by origin (lake vs. stream,

lower expression in transplants than in wild fish (see Supporting infor-

explains 86% of variance). LDA2 separated fish based on their trans-

mation for full list of genes and statistics). Stream transplants have

plant destination (explains 10% of variance). LDA3 roughly separates

lower expression of immune genes with known function including the

native/non-native status (explains 3.5% of variance, Figure 6;

complement system (complement 3, 8 and 9), a leucocyte-derived

Fig. S3). We plotted a vector from the mean of each resident eco-

chemotaxin (lect2l) and three fibrinogen genes (alpha, beta and

type at home, to the mean expression of the same ecotype when

gamma). In addition, two coagulation factor genes are lower in natives

moved into a new environment. The vector showing the expression

(factor 13 and 7i) than wild stream fish. The cage effect for stream

change of lake fish is almost in exactly the opposite direction from

fish is partially confounded, however, with genotype. The stream

the expression change of stream fish (~180°, visually). In each case,

transplants were from 1.5 km downstream of the cage site, whereas

fish moved into a new habitat converged on the expression profile

the wild fish were collected among the cages, 100 m downstream

of their new neighbours along LD2 (but not along LD1 or LD3). Lake

from the lake. So, differences between wild stream and transplanted

fish moved into the stream actually overshot the stream expression

stream fish may be genetic rather than exclusively a plastic response

profile, resulting in a much larger reaction norm vector than stream

to caging. There were almost no genes (only 2) that showed signifi-

fish moved into the lake (LD2 ~origin + destination + origin:destina-

cant effects of caging in both the lake and in the stream, indicating

tion) and found a significant effect of the interaction (p  .001).

that there is no generic transcriptomic response to caging (Fig. S2).

Because of this overshooting, both the lake-to-stream migrants and

Our coexpression analysis of transplanted fish revealed signifi-

stream-to-lake migrants were significantly different (for LD2) from

cant correlations between traits unique to this subset of fish and

the resident “target.” We conclude that immigrant stickleback par-

modules of gene expression. For example, treatment (transplanted

tially converge on native expression profiles after emigration to a

LOHMAN

|

ET AL.

4663

Module−trait relationships Lake vs Stream MEblack

−

−

−

−

−

0.19

−

−

−

−

−

−

−

−

−

−

−

− −0.27 − −0.27

−0.25 −

−

−

−

−

0.24

−

−

−

−

−

−

−

−

−

−

−

0.24 −0.23 − −0.24

−

−

−

−

−

−

−

−

0.18

−

−

−

−

−

−

−

0.27 −0.23 − −0.25

MEyellow

−0.26 − −0.23 −

−

−

−

−

−

−

−

−

−

−

−

−

−

−

− −0.27 − −0.27

MEpurple

−0.19 −

MEpink

MEmagenta

−0.3 −0.2

−

−

−

−

−

− −0.22 −

−

− −0.24 −

−

−

−

−

− −0.35 − −0.35

MEbrown

−

−

−

−0.2

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

MEgreen

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

MEblue

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

−

0.22

−

−

−

−

−

−

−0.2

−

−

−

−

−

−

−

−

−

−

−

−

0.35

−

0.35

−

−

−

0.2

−

−

− −0.19 −

−

−

−

− −0.21−0.2 60.26

−

−

−

0.25

−

−

0.25

−

−

−

−

−

−

MEgreenyellow

MEred

0.36 0.18 0.27

−0.18−0.25 −

−0.28 − −0.26 −

−

0.28

−

−

0.2

0

−0.2

−0.4

− −0.27−0.21−0.28

O r lv igin ic Bo w i dy d t h G Ga dep ill r a pe t h w k G e r n idth ill ra u m ke b e rl r en gt h an S y e an Ce x s a n yNe tod yD m e ip ato lo d an any sto e yP Tr mu r o em m te o c ato e p de an ha Bu yE lus no xt de ra an ern o r yIn a l Bu ter n o na de l Pa rin ra an a s Sh ite y G N a D um nn iv ut be on ers r o Di i t y f M ve H rsit y C al le m les hc L m D1 hc m LD2 hc D C A

MEturquoise

0.4

Pe

F I G U R E 2 WGCNA reveals correlations between modules of coexpressed genes and traits in wild lake and stream fish. Cell values are Pearson correlation coefficients. Only correlations with p values <.1 are presented. Modules shown are the same as Figure 5. Pelvic width is the width of the pelvic girdle, body depth is the distance from the base of the first dorsal spine and the anterior point of the pelvic girdle, gape width is the distance between mouth corners. Columns labelled “Any” refer to any species of a broad group

−

new habitat and that lake fish exhibit stronger plasticity. The latter

adaptive in several ways. First, generic stress responses could be used

finding matches the greater plasticity of lake fish in our gene-by-

to protect immigrants’ from unfamiliar environmental conditions, par-

gene analysis with DESeq2 (above; Figure 7).

asites, low energy income, etc. Such genes might be upregulated for

We also considered convergence at the individual gene level.

all migrants regardless of origin or destination. Second, organisms

Using the DESeq2 linear model estimates, we found that destination

may adjust their expression to better match the local environment,

effects were positively correlated with origin effects (r = .67, Fig-

converging on residents of the migrants’ new habitat. Such genes

ure 7). That is, transcripts that were more abundant in lake natives

would exhibit transcriptional convergence, a main effect of destina-

were also more abundant in fish placed in lake cages, and vice versa

tion habitat. Of course, we also must acknowledge that transcriptional

for stream-biased transcripts (Figure 7). This implies that for many

plasticity can be maladaptive: a signal of stress or poor condition, or a

genes, expression differences between the native populations are

misguided response to an unfamiliar environmental cue.

recapitulated by plastic responses to animals’ recent environment.

To look for static and plastic responses in gene expression asso-

The observed destination origin relationship has a slope less than 1

ciated with emigration, we tested for differences in gene expression

(0.42, p  .001) indicating that the plasticity is not, however, com-

among stickleback reciprocally transplanted between two adjoining

plete, which fits with the fact that the major LDA axis still separates

habitats containing genetically divergent populations. We found

lake vs. stream natives and explains more variation than the second

expression differences between these populations, and changes in

LDA axis that measured plasticity.

response to emigration, at the level of individual genes, gene coexpression and the whole transcriptome [see Jones et al., (2012) for similar discussion contrasting observations of wild marine and fresh-

4 | DISCUSSION

water stickleback gene expression and (Ishikawa et al., 2017) for an eQTL study testing for the relative contribution of cis and trans reg-

Organisms’ adaptation to their native habitats means that migrants

ulatory elements and their connection with genomic islands of adap-

will often be maladapted to novel environments. One way that

tation in adaptation to freshwater].

migrants may be able to ameliorate stressors of new habitats is by modulating gene expression. Prior studies have used reciprocal transplants to uncover plasticity in select candidate genes, but this approach could miss a myriad of responses to the environment

4.1 | There are constitutive differences in gene expression between lake and stream stickleback

(although see (Ghalambor et al., 2015; Kenkel & Matz, 2016) for a

Although Roberts Lake and Stream are adjoining habitats that permit

transcriptome-wide approach). Transcriptional plasticity could be

easy movement of stickleback between sites, the resident stickleback

|

0.0

−0.10

0

3

4

5

6 7 8 N um ber of alleles

9

10

Infection by any Nematodes

0

0.0

0.1

0.2

0.3

6

MEpurple expression

5

(c)

4

0.1

–0.05

–0.1

3

0.05

1.0

0.00

0.5

S hannon diversity

−0.05

Intensity of MEgreenyellow expression

0.05

2

0.10

−0.1

(b)

1

1.5

(a)

ET AL.

0

Expression of greenyellow is correlated with MHC alle # and Shannon diversity

1 2 3 4 5 6 7 8 9 10 11 12

LOHMAN

Infection by any proteocephallus

4664

–0.1

0.0

0.1

0.2

0.3

MEred expression

F I G U R E 3 Module trait correlations from WGCNA. (a) Heat map of MEgreenyellow expression shows negative correlation with Shannon diversity of parasite infection and positive correlation with MHC allele number. (b) Increased expression of MEpurple is correlated with decreased infection by nematodes. (c) Increased expression of MEred is correlated with increased infection by proteocephalus and stream fish. Macrophages contribute to initiation of immune

g6pca.6

4.5

6 5

the tapeworm Schistocephalus solidus (Kurtz et al., 2006), whose

4.0

4

defences against a variety of parasites including but not limited to infectious procercoids are deposited by loons and mergansers (which prefer lakes over streams) and carried by zooplankton (which are

3

3.5

4.0

5.0

6.0

5.0

dhx58

3.0

3.0

Variance stabilized count (Expression)

cyp24a1

Lake

Stream

Lake

Destination

Stream Destination

Lake

Stream Destination

more abundant in lakes than streams). MHC class II is another parasite defence-related GO category which is different between lake

F I G U R E 4 Reaction norm plots of genes significant for interaction between origin and destination. X-axis is destination, teal line is lake and magenta line is stream origin/genotype. Vertical line indicates standard error

and stream. Prior work in the Roberts Lake stickleback has revealed that MHC II allele frequencies differ between this particular lake and stream (Stutz & Bolnick, 2014), as well as many other lake–stream pairs (Eizaguirre et al., 2010; Kurtz et al., 2006; Wegner et al., 2006). Furthermore, individuals who carry local MHC alleles are

populations are genetically distinct. Fish from this lake and stream

more heavily infected with parasites than individuals carrying foreign

differ in a range of morphological and parasitological traits (Berner

MHC alleles (Bolnick & Stutz, 2017). Our WGCNA results suggest

et al., 2009; Bolnick & Stutz, 2017; Oke et al., 2016; Stutz & Bol-

that MHC allele diversity and parasite diversity are negatively corre-

nick, 2017; Weber, Bradburd et al., 2017), as is true for many such

lated with each other and jointly associated with a set of coex-

lake–stream pairs (Stuart et al., 2017). Given these genetic and phe-

pressed genes. Specifically, the greenyellow module has a negative

notypic differences, we expected to find differences in gene expres-

correlation with parasite diversity and a positive correlation with the

sion between these populations. Approximately 7% of the 9,748

number of MHC alleles (Figures 2 and 3a). While this result stands

genes in our transcriptome data set exhibited between-population

out as support for a large body of theory (Eizaguirre & Lenz, 2010;

differences in relative abundance.

Spurgin & Richardson, 2010) and agrees with prior empirical evi-

Some of these differences fit well within the existing literature

dence (Piertney, Telfer, & Oliver, 2009; Wegner et al., 2003), we

of lake–stream divergence. For example, our GO enrichment results

would have expected greater correlations between modules and par-

suggest that macrophage differentiation is different between lake

asite infection. However, this lack of correlation is likely due to

LOHMAN

|

ET AL.

Module−trait relationships Transplant

4665

plasticity. Genes more highly expressed in lake (stream) natives were also upregulated in all fish placed in lake (stream) cages (Figure 7). If

−

−

−

MEpink

−0.25

−0.21

−

−

MEmagenta

−0.3

−0.22

−0.21

−

MEyellow

−0.26

−

−

−

MEpurple

−0.19

−

0.23

−

MEbrown

−

−

−

−

MEgreen

−

−

−

−

MEblue

−

−

−0.2

−

MEturquoise

0.36

−

−

−

MEgreenyellow

−0.18

−0.22

−

−

MEred

−0.28

−

0.18

0.2

gin and destination. This is consistent with prior observations that there are no interactions effects between origin and destination for

s

0.4

the case. So, this correlation between origin and destination effects the same direction, to between-ecotype differences in expression. The fact that the origin–destination effect correlation has a slope less than 1 confirms the statement, above, that heritable (origin)

0

effects were somewhat stronger than the environmental (destination) effects. Moreover, the paucity of genes in the top-left and bot-

−0.2

tom-right quadrants of Figure 7 suggests that remarkably few genes exhibited plastic responses that opposed the heritable lake–stream

−0.4

differences. Very few genes (10) were significant for the interaction of ori-

Le

parasite load, survival, growth or condition in this experiment (Bol-

in ng

e

nick & Stutz, 2017). The few interactions that do exist follow two distinct patterns. First, some genes were downregulated after indi-

ha

ha

C

C

ment), we would expect to see no origin effect at all, which is not suggests that heritable and plastic differences jointly contribute, in

0.2

ng

as M

ng

D

e

es

in

tin

O

at

r ig

io

in

n

−

th

expression was exclusively plastic (on the time-scale of our experi-

MEblack

viduals were placed in a foreign habitat. Second, other genes were

F I G U R E 5 WGCNA reveals correlations between suites of coexpressed genes and traits in transplanted fish. Cell values are Pearson correlation coefficients. Only correlations with p values <.1 are presented. Modules shown are the same as Figure 2

more highly expressed by lake fish, but also showed stronger plastic downregulation in lake fish placed in the foreign stream habitat. The absence of genes which were more highly expressed by stream fish, regardless of habitat, is notable. Some of the interaction genes we do detect may be involved in ROS production and antiviral response, both of which may be potentially important to fitness.

sparse and overdispersed parasite infections, which make correla-

For example, ROS production was recently shown to be a heritable

tions difficult to estimate well.

response to infection by S. solidus (Weber, Steinel, Shim, & Bolnick, 2017). We observed no cases where expression was higher in foreign

4.2 | Transplantation into alternate habitats reveals static and plastic gene expression

habitats. While this could be a product of the low number of

For

origin

in future work. Intuitively, we would have expected transplanted

accounted for more expression variation (507 genes) than did desti-

fish in either direction to upregulate stress genes, but this appar-

nation (111 genes, Figure 7). The main effect of origin represents

ently did not occur. Perhaps the absence of interaction effects on

persistent between-population differences no matter which habitat

genes here is because plasticity reinforced between-ecotype differ-

the fish were caged in. Thus, we interpret the main effect of origin

ences. The paucity of interaction effects may also be a conse-

as a probable signal of static genetic differences in expression, insen-

quence

sitive to the environment. However, we also found significant desti-

expression levels converts multiplicative (interaction) effects into

nation effects for a subset of genes, indicating appreciable plasticity

additive effects, which can reduce power or even completely

our

experimentally

transplanted

fish,

individuals’

interaction genes, this pattern is surprising and worth considering

of

our

analytical

technique:

log

transformation

of

in gene expression in response to sticklebacks’ recent (cage) environ-

obscure our ability to detect significant interactions between origin

ment. That is, expression of certain genes was higher in lake-caged

and destination effects. Nevertheless, other reciprocal transplant

fish than stream-caged fish, regardless of their origin. We infer that

studies using large-scale RNAseq have found more interaction

shifts in gene expression are the result of plasticity rather than

genes and this seems to be relatively common (Lovell et al., 2016;

selection because the expression profile of immigrants falls outside

Reid et al., 2016).

that of the natives in PCA space (Figure 6). This plasticity is consis-

Our WGCNA analysis revealed only weak correlations between

tent with prior evidence that morphological plasticity contributes to

origin and phenotypes unique to transplanted fish. For example,

phenotypic differences between the Roberts Lake and Stream stick-

change in mass and length. However, it is interesting to note that

leback (Oke et al., 2016).

changes in mass and length are most correlated with different mod-

Notably, there was a positive correlation between origin effect

ules. This may suggest a change in condition (loss of mass

and destination effect (r = .67). We therefore infer that the heritable

but increase in length due to growth but poor foraging efficiency,

lake-to-stream differences were at least partly recapitulated by

Figure 5).

4666

|

LOHMAN

ET AL.

Convergence in LDA space

2 0 −2 −4

LDA2 (Destination)

4

Destination lake stream

Origin Lake Stream

−5

5

0

10

LDA1 (Origin) F I G U R E 6 Convergence of immigrant expression profiles towards native expression profiles in transplanted fish. Red arrows are drawn between the means of each distribution. Fish originating from the lake move farther along LD2 than stream fish (two-factor ANOVA, p  .001). LDA1 explains 86% of the total variance while LDA2 explains 10%

4.3 | On the whole-transcriptome level, lake fish are more plastic than stream fish At the whole-transcriptome level, we again observe substantial and

placed in lake cages fall short of the optimum expression in the lake (Figure 6). We therefore conclude that transcriptomic plasticity is incomplete (LD1 remains intact and explains the most variance), and differs between lake and stream ecotypes. This result implies that

persistent differences between the expression profiles of lake and

sticklebacks’ transcriptional reaction norms may be evolving as they

stream fish, captured by LD axis 1. However, along LD2, we observe

adapt to different habitats. However, because we used wild-caught

substantial plastic convergence of immigrant fish towards the

rather than laboratory-raised fish for this experiment, we cannot rule

expression profile of their new population (Figure 6). Interestingly,

out effects of early rearing environment, and hence cannot defini-

we also observed convergence in parasite community composition in

tively ascribe a genetic cause to the different reaction norms of lake

this same experiment. Lake and stream natives carried distinct para-

and stream fish.

site communities, and individuals transplanted to the neighbouring

Our results lend additional support to an emerging insight that

habitat exhibited an intermediate parasite community (Bolnick &

transcriptomic plasticity may play a substantial role in migrants’

Stutz, 2017).

adaptation to novel environments. This has been very extensively

Our analysis suggests that fish from the lake exhibit a more plas-

explored in experimental settings in the laboratory, where organisms

tic response to being transplanted into the stream, compared to

may be exposed to alternative environmental conditions (often a sin-

stream fishes’ more limited plasticity when placed in the lake. This

gle variable such as salinity, temperature or a toxin). Many studies

transcriptome-wide analysis is consistent with our single-gene analy-

find plastic responses in candidate genes, or a subset of the tran-

ses which also found that lake natives tended to show greater plas-

scriptome, in response to such experimental treatments (Morris

ticity in response to transplantation. Assuming fish caged in their

et al., 2014; Reid et al., 2016; Velotta et al., 2017; Whitehead,

native habitat adopt a locally optimal expression profile, we infer

Roach, Zhang, & Galvez, 2011). Often, these plastic responses are

that lake sticklebacks’ strong plastic response is actually excessive,

genotype dependent, with one population exhibiting a stronger

overshooting the stream profile along LD2. In contrast, stream fish

response than another (e.g., PCB-tolerant killifish are less plastic than

LOHMAN

|

ET AL.

4667

Gene−by−gene convergence Transplanted Fish Only

6 4 2 0 −2 −6

−4

Destination (Lake−Stream)

8

Origin effect (473) Destination effect (77) Both (34) Not significant (9194)

−25

−20

−15

−10

−5

0

5

Origin (Lake−Stream) F I G U R E 7 Gene-by-gene convergence among transplanted fish. We included only transplanted fish in a linear model in DESeq2: with expression of each focal gene as a function of origin + destination + origin:destination. X- and Y-axis are Log2 fold changes between lake and stream fish by origin, and destination, respectively. Points are coloured when q value <.05, and colour-coded based on which effect(s) were significant. The red dashed line is 1:1, helping to visualize that the main trend has a slope <1, indicating that plasticity (destination) effects are weaker than origin effects PCB-susceptible populations (Reid et al., 2016)). Fewer studies have

Van Buskirk, 2002). Our result is thus somewhat puzzling, in that we

examined transcriptomic plasticity of migrants in natural settings.

observe greater transcriptomic plasticity in lake fish, which inhabit

Kenkel and Matz (2016) subjected corals to a reciprocal transplant

the more temporally stable habitat. While stream habitats are gener-

experiment across a temperature gradient, and also found transcrip-

ally very diverse (flow regime, overhead foliage, substrate, spatial dis-

tomic convergence of migrants towards residents, as we do. They

tribution of prey), lake habitats generally have large and smooth

also found that one genotype was more transcriptionally plastic than

transitions between any variation in environmental variables (and in

the other, as we do. In a similar design to ours (Ghalambor et al.,

most cases very little variation (Ahmed, Thompson, Bolnick, & Stuart,

2015) transplanted fish between two locations in streams and mea-

2017; Stuart et al., 2017). However, lakes may be less predictable in

sured subsequent changes in gene expression. Their results high-

other ways. For instance, lake stickleback consistently harbour more

lighted that adaptive plasticity decreases the impact of directional

diverse parasite communities (Bolnick & Stutz, 2017; Stutz & Bol-

selection, and thus the pace of evolution. Under this paradigm, we

nick, 2017), and so may have evolved greater immunological plastic-

might expect that expression profiles of stream fish might evolve

ity to handle an unpredictable suite of pathogens and helminths.

more rapidly than lake fish because of their reduced plasticity.

In conclusion, we see extensive gene expression differences

A large body of existing empirical and theoretical studies suggest

between genetically divergent stickleback populations inhabiting

that increased plasticity should evolve in more temporally or spatially

adjoining habitats but connected by gene flow (Weber, Bradburd

heterogeneous habitats (Auld & Relyea, 2011; Baythavong, 2011;

et al., 2017). But, for many genes, transcript abundance is highly

Davidson et al., 2011; Dudley & Schmitt, 1996; Murren et al., 2015;

plastic. Fish that disperse into the adjoining foreign habitat will

4668

|

partially converge on the gene expression profile typical of their new habitat. This suggests that expression plasticity can soften the impact of immigration into an unfamiliar habitat. Nonetheless, lake and stream fish differed in survival rates and parasite infection rates in our study, implying that this expression plasticity is not fast or extensive enough to fully homogenize the lake and stream fish performance.

ACKNOWLEDGEMENTS Collection and animal handling were approved by the University of Texas Institutional Animal Use and Care Committee (Protocol # 07032201) and a Scientific Fish Collection Permit from the Ministry of the Environment of British Columbia (NA07-32612). We wish to thank Chad Brock and Kelsey Jiang for their assistance with collecting stickleback head kidneys during the reciprocal transplant experiment. The staff of the Genome Sequencing and Analysis Facility at the University of Texas at Austin provided technical support. This work was supported by the Howard Hughes Medical Institute (DIB) and the David and Lucille Packard Foundation (DIB).

DATA ACCESSIBILITY Meta data, parasite data, code for processing raw reads, code for statistical analysis and plotting are located in DRYAD entry https://doi. org/10.5061/dryad.mk8ns. Raw reads are available for download via “wget http://web.corral.tacc.utexas.edu/Lohman_et_al_2017_Molecu larEcology/” from the terminal. The iRNAseq pipeline is available here: https://github.com/z0on/tag-based_RNAseq. The GO_MWU software is available here: https://github.com/z0on/GO_MWU.

AUTHOR CONTRIBUTIONS BKL built and sequenced TagSeq libraries, and performed gene expression data analysis. WES performed the experiment which generated the samples. BKL and DIB wrote the manuscript. All authors approved the final version. The authors declare no conflict of interests.

REFERENCES Ahmed, N. I., Thompson, C., Bolnick, D. I., & Stuart, Y. E. (2017). Brain morphology of the threespine stickleback (Gasterosteus aculeatus) varies inconsistently with respect to habitat complexity: A test of the clever foraging hypothesis. Ecology and Evolution, 7, 3372–3380. Annalora, A. J., Goodin, D. B., Hong, W.-X., Zhang, Q., Johnson, E. F., & Stout, C. D. (2010). Crystal structure of CYP24A1, a mitochondrial cytochrome P450 involved in vitamin D metabolism. Journal of Molecular Biology, 396, 441–451. Auld, J. R., & Relyea, R. A. (2011). Adaptive plasticity in predator-induced defenses in a common freshwater snail: Altered selection and mode of predation due to prey phenotype. Evolutionary Ecology, 25, 189–202. Baker, J., & Foster, S. (2002). Phenotypic plasticity for life history traits in a stream population of the threespine stickleback, Gasterosteus aculeatus L. Ecology of Freshwater Fish, 11, 20–29.

LOHMAN

ET AL.

Baythavong, B. S. (2011). Linking the spatial scale of environmental variation and the evolution of phenotypic plasticity: Selection favors adaptive plasticity in fine-grained environments. The American Naturalist, 178, 75–87. Berner, D., Grandchamp, A.-C., & Hendry, A. P. (2009). Variable progress toward ecological speciation in parapatry: Stickleback across eight lake-stream transitions. Evolution, 63, 1740–1753. € ANen, € Berner, D., Kaeuffer, R., Grandchamp, A. C., RAS K., Hendry, A. P., & Raeymaekers, J. A. M. (2011). Quantitative genetic inheritance of morphological divergence in a lake-stream stickleback ecotype pair: Implications for reproductive isolation. Journal of Evolutionary Biology, 9, 1975–1983. Berner, D., & Thibert-Plante, X. (2015). How mechanisms of habitat preference evolve and promote divergence with gene flow. Journal of Evolutionary Biology, 28, 1641–1655. Bolnick, D. I. (2004). Can intraspecific competition drive disruptive selection? An experimental test in natural populations of sticklebacks. Evolution, 87, 608–618. Bolnick, D. I., Snowberg, L. K., Patenia, C., Stutz, W. E., Ingram, T., & Lau, O. L. (2009). Phenotype-dependent native habitat preference facilitates divergence between parapatric lake and stream stickleback. Evolution, 63, 2004–2016. Bolnick, D. I., & Stutz, W. E. (2017). Frequency dependence limits divergent evolution by favouring rare immigrants over residents. Nature, 546, 285–288. Conover, D. O., & Present, T. M. (1990). Countergradient variation in growth rate: Compensation for length of the growing season among Atlantic silversides from different latitudes. Oecologia, 83, 316–324. Davidson, A. M., Jennions, M., & Nicotra, A. B. (2011). Do invasive species show higher phenotypic plasticity than native species and if so, is it adaptive? A meta-analysis. Ecology Letters, 14, 419–431. Davis, J. M., & Stamps, J. A. (2004). The effect of natal experience on habitat preferences. Trends in Ecology & Evolution, 19, 411–416. Day, T., & McPhail, J. (1996). The effect of behavioural and morphological plasticity on foraging efficiency in the threespine stickleback (Gasterosteus sp.). Oecologia, 108, 380–388. Dixon, G. B., Davies, S. W., Aglyamova, G. A., Meyer, E., Bay, L. K., & Matz, M. V. (2015). Genomic determinants of coral heat tolerance across latitudes. Science, 348, 1460–1462. Dudley, S. A., & Schmitt, J. (1996). Testing the adaptive plasticity hypothesis: Density-dependent selection on manipulated stem length in Impatiens capensis. The American Naturalist, 147, 445–465. Edelaar, P., & Bolnick, D. I. (2012). Non-random gene flow: An underappreciated force in evolution and ecology. Trends in Ecology & Evolution, 27, 659–665. Edelaar, P., Siepielski, A. M., & Clobert, J. (2008). Matching habitat choice causes directed gene flow: A neglected dimension in evolution and ecology. Evolution, 62, 2462–2472. Eizaguirre, C., & Lenz, T. L. (2010). Major histocompatibility complex polymorphism: Dynamics and consequences of parasite-mediated local adaptation in fishes. Journal of Fish Biology, 77, 2023–2047. Eizaguirre, C., Lenz, T. L., Kalbe, M., & Milinski, M. (2012). Divergent selection on locally adapted major histocompatibility complex immune genes experimentally proven in the field. Ecology Letters, 15, 723–731. Eizaguirre, C., Lenz, T. L., Sommerfeld, R. D., Harrod, C., Kalbe, M., & Milinski, M. (2010). Parasite diversity, patterns of MHC II variation and olfactory based mate choice in diverging three-spined stickleback ecotypes. Evolutionary Ecology, 25, 605–622. Eizaguirre, C., Lenz, T. L., Sommerfeld, R. D., Harrod, C., Kalbe, M., & Milinski, M. (2011). Parasite diversity, patterns of MHC II variation and olfactory based mate choice in diverging three-spined stickleback ecotypes. Evolutionary Ecology, 25, 605–622. Feulner, P. G., Chain, F. J., Panchal, M., Huang, Y., Eizaguirre, C., Kalbe, M., . . . Reusch, T. B. (2015). Genomics of divergence along a

LOHMAN

ET AL.

continuum of parapatric population differentiation. PLOS Genetics, 11, e1004966. Fischer, U., Koppang, E. O., & Nakanishi, T. (2013). Teleost T and NK cell immunity. Fish & Shellfish Immunology, 35, 197–206. Fischer, U., Utke, K., Somamoto, T., Kollner, B., Ototake, M., & Nakanishi, T. (2006). Cytotoxic activities of fish leucocytes. Fish & Shellfish Immunology, 20, 209–226. Ghalambor, C. K., Hoke, K. L., Ruell, E. W., Fischer, E. K., Reznick, D. N., & Hughes, K. A. (2015). Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature. Nature, 525, 372– 375. Ghalambor, C. K., McKay, J. K., Carroll, S. P., & Reznick, D. N. (2007). Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. Functional Ecology, 21, 394–407. Gibbons, T. C., Metzger, D. C., Healy, T. M., & Schulte, P. M. (2017). Gene expression plasticity in response to salinity acclimation in threespine stickleback ecotypes from different salinity habitats. Molecular Ecology, 26, 2711–2725. Hanson, D., Moore, J. S., Taylor, E., Barrett, R., & Hendry, A. (2016). Assessing reproductive isolation using a contact zone between parapatric lake-stream stickleback ecotypes. Journal of Evolutionary Biology, 29, 2491–2501. Hendry, A. P., Taylor, E. B., & McPhail, J. D. (2002). Adaptive divergence and the balance between selection and gene flow: lake and stream stickleback in the Misty system. Evolution, 56, 1199–1216. Hereford, J. (2009). A quantitative survey of local adaptation and fitness trade-offs. The American Naturalist, 173, 579–588. Ishikawa, A., Kusakabe, M., Yoshida, K., Ravinet, M., Makino, T., Toyoda, A., . . . Kitano, J. (2017). Different contributions of local-and distantregulatory changes to transcriptome divergence between stickleback ecotypes. Evolution, 71, 565–581. Jiang, Y., Peichel, C. L., Ling, F., & Bolnick, D. I. (2017). Sensory trait variation contributes to biased dispersal of threespine stickleback in flowing water. Journal of Evolutionary Biology, 30, 681–685. Jiang, Y., Torrance, L., Peichel, C. L., & Bolnick, D. I. (2015). Differences in rheotactic responses contribute to divergent habitat use between parapatric lake and stream threespine stickleback. Evolution, 69, 2517–2524. Jombart, T., Devillard, S., & Balloux, F. (2010). Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genetics, 11, 94. Jones, F. C., Grabherr, M. G., Chan, Y. F., Russell, P., Mauceli, E., Johnson, J., . . . Birney, E. (2012). The genomic basis of adaptive evolution in threespine sticklebacks. Nature, 484, 55–61. Kaeuffer, R., Bolnick, D., Hendry, A. P., & Peichel, C. L. (2012). Convergence and non-convergence in ecological, phenotypic, and genetic divergence across replicate population pairs of lake and stream stickleback. Evolution, 66, 402–418. Kalbe, M., & Kurtz, J. (2006). Local differences in immunocompetence reflect resistance of sticklebacks against the eye fluke Diplostomum pseudospathaceum. Parasitology, 132, 105–116. Kaufmann, J., Lenz, T. L., Kalbe, M., Milinski, M., & Eizaguirre, C. (2017). A field reciprocal transplant experiment reveals asymmetric costs of migration between lake and river ecotypes of three-spined sticklebacks (Gasterosteus aculeatus). Journal of Evolutionary Biology, 30, 938–950. Kawecki, T. J., & Ebert, D. (2004). Conceptual issues in local adaptation. Ecology Letters, 7, 1225–1241. Kenkel, C. D., & Matz, M. V. (2016). Gene expression plasticity as a mechanism of coral adaptation to a variable environment. Nature Ecology & Evolution, 1, 0014. Kurtz, J., Wegner, K. M., Kalbe, M., Reusch, T. B., Schaschl, H., Hasselquist, D., & Milinski, M. (2006). MHC genes and oxidative stress in

|

4669

sticklebacks: An immuno-ecological approach. Proceedings of the Biological Sciences/The Royal Society, 273, 1407–1414. Langfelder, P., & Horvath, S. (2008). WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics, 9, 559. Leavy, O. (2012). Antiviral immunity: LGP2 rigs CD8+ T cells for survival. Nature Reviews Immunology, 12, 616–617. Leder, E. H., McCairns, R. S., Leinonen, T., Cano, J. M., Viitaniemi, H. M., a, J. (2014). The evolution and adaptive potential Nikinmaa, M., . . . Meril€ of transcriptional variation in sticklebacks—Signatures of selection and widespread heritability. Molecular Biology and Evolution, 32, 674–689. Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends in Ecology & Evolution, 17, 183–189. Lenz, T. L., Eizaguirre, C., Rotter, B., Kalbe, M., & Milinski, M. (2013). Exploring local immunological adaptation of two stickleback ecotypes by experimental infection and transcriptome-wide digital gene expression analysis. Molecular Ecology, 22, 774–786. Lohman, B. K., Steinel, N., Weber, J. N., & Bolnick, D. (In review). The role of gene expression in the recent evolution of host resistance in a model host parasite system. Evolution. Lohman, B. K., Weber, J. N., & Bolnick, D. I. (2016). Evaluation of TagSeq, a reliable low-cost alternative for RNAseq. Molecular Ecology Resources, 16, 1315–1321.
pez-Maury, L., Marguerat, S., & B€ Lo ahler, J. (2008). Tuning gene expression to changing environments: From rapid responses to evolutionary adaptation. Nature Reviews Genetics, 9, 583–593. Love, M., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. Lovell, J. T., Schwartz, S., Lowry, D. B., Shakirov, E. V., Bonnette, J. E., Weng, X., . . . Jenkins, J. (2016). Drought responsive gene expression regulatory divergence between upland and lowland ecotypes of a perennial C4 grass. Genome Research, 26, 510–518. McCairns, R. J., & Bernatchez, L. (2010). Adaptive divergence between freshwater and marine sticklebacks: Insights into the role of phenotypic plasticity from an integrated analysis of candidate gene expression. Evolution; International Journal of Organic Evolution, 64, 1029–1047. Meyer, E., Aglyamova, G. V., & Matz, M. V. (2011). Profiling gene expression responses of coral larvae (Acropora millepora) to elevated temperature and settlement inducers using a novel RNA-Seq procedure. Molecular Ecology, 20, 3599–3616. Meyer, E., Davies, S., Wang, S., Willis, B. L., Abrego, D., Juenger, T. E., & Matz, M. V. (2009). Genetic variation in responses to a settlement cue and elevated temperature in the reef-building coral Acropora millepora. Marine Ecology Progress Series, 392, 81–92. Morris, M. R., Richard, R., Leder, E. H., Barrett, R. D., Aubin-Horth, N., & Rogers, S. M. (2014). Gene expression plasticity evolves in response to colonization of freshwater lakes in threespine stickleback. Molecular Ecology, 23, 3226–3240. Moser, D., Frey, A., & Berner, D. (2016). Fitness differences between parapatric lake and stream stickleback revealed by a field transplant. Journal of Evolutionary Biology, 29, 711–719. Murphy, K. M. (2011). Janeway’s immunobiology. New York, NY: Garland Science. Murren, C. J., Auld, J. R., Callahan, H., Ghalambor, C. K., Handelsman, C. A., Heskel, M. A., . . . Pfennig, D. W. (2015). Constraints on the evolution of phenotypic plasticity: Limits and costs of phenotype and plasticity. Heredity, 115, 293–301. Nosil, P., Vines, T. H., & Funk, D. J. (2005). Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution, 59, 705–719. Oke, K. B., Bukhari, M., Kaeuffer, R., Rolshausen, G., R€ as€ anen, K., Bolnick, D. I., . . . Hendry, A. P. (2016). Does plasticity enhance or dampen phenotypic parallelism? A test with three lake–stream stickleback pairs. Journal of Evolutionary Biology, 29, 126–143.

4670

|

Pfennig, D. W. (1992). Polyphenism in spadefoot toad tadpoles as a logically adjusted evolutionarily stable strategy. Evolution, 46, 1408– 1420. Piertney, S. B., Telfer, S., & Oliver, M. K. (2009). Major histocompatibility complex (MHC) heterozygote superiority to natural multi-parasite infections in the water vole (Arvicola terrestris). Proceedings of the Royal Society B: Biological Sciences, 276, 1119–1128. Queitsch, C., Sangster, T. A., & Lindquist, S. (2002). Hsp90 as a capacitor of phenotypic variation. Nature, 417, 618–624. R Development Core Team. (2007). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Reid, N. M., Proestou, D. A., Clark, B. W., Warren, W. C., Colbourne, J. K., Shaw, J. R., . . . Crawford, D. L. (2016). The genomic landscape of rapid repeated evolutionary adaptation to toxic pollution in wild fish. Science, 354, 1305–1308. Reusch, T. B. H., Wegner, K. M., & Kalbe, M. (2001). Rapid genetic divergence in postglacial populations of threespine stickleback (Gasterosteus aculeatus): The role of habitat type, drainage and geographical proximity. Molecular Ecology, 10, 2435–2445. Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., & Smyth, G. K. (2015). limma powers differential expression analyses for RNAsequencing and microarray studies. Nucleic Acids, 47, e47–e47. Robertson, S., Bradley, J. E., & MacColl, A. D. (2016). Measuring the immune system of the three-spined stickleback–investigating natural variation by quantifying immune expression in the laboratory and the wild. Molecular Ecology Resources, 16, 701–713. Roesti, M., Hendry, A. P., Salzburger, W., & Berner, D. (2012). Genome divergence during evolutionary diversification as revealed in replicate lake-stream stickleback population pairs. Molecular Ecology, 21, 2852– 2862. Rutherford, S. L., & Lindquist, S. (1998). Hsp90 as a capacitor for morphological evolution. Nature, 396, 336–342. Scharsack, J. P., Kalbe, M., Derner, R., & Millinski, M. (2004). Modulation of granulocyte responses in three-spined sticlebacks Gasterosteus aculeatus infected with the tapeworm Schistocephalus solidus. Diseases of Aquatic Organisms, 59, 141–150. Scharsack, J. P., Kalbe, M., Harrod, C., & Rauch, G. (2007). Habitat-specific adaptation of immune responses of stickleback (Gasterosteus aculeatus) lake and river ecotypes. Proceedings of the Royal Society B: Biological Sciences, 274, 1523–1532. Scharsack, J. P., Koch, K., & Hammerschmidt, K. (2007). Who is in control of the stickleback immune system: Interactions between Schistocephalus solidus and its specific vertebrate host. Proceedings of the Royal Society B: Biological Sciences, 274, 3151–3158. Schlichting, C. D., & Pigliucci, M. (1998). Phenotypic evolution: A reaction norm perspective. Sunderland, MA: Sinauer Associates Incorporated. Spurgin, L. G., & Richardson, D. S. (2010). How pathogens drive genetic diversity: MHC, mechanisms and misunderstandings. Proceedings of the Royal Society B: Biological Sciences, 277, 979–988. Stamps, J., & Davis, J. M. (2006). Adaptive effects of natal experience on habitat selection by dispersers. Animal Behaviour, 72, 1279–1289. Stuart, Y. E., Veen, T., Weber, J. N., Hanson, D., Ravinet, M., Lohman, B. K., . . . Ahmed, N. (2017). Contrasting effects of environment and genetics generate a continuum of parallel evolution. Nature Ecology & Evolution, 1, 0158. Stutz, W. E., & Bolnick, D. I. (2014). A stepwise threshold clustering (STC) method to infer genotypes from error-prone next-generation sequencing of multi-allele genes such as the major histocompatibility complex (MHC). PLoS ONE, 9, e100587. Stutz, W. E., & Bolnick, D. I. (2017). Natural selection on MHC IIb in parapatric lake and stream stickleback: Balancing, divergent, both or neither? Molecular Ecology, 1–15. https://doi.org/10.1111/mec. 14158.

LOHMAN

ET AL.

Stutz, W. E., Schmerer, M., Coates, J. L., & Bolnick, D. I. (2015). Amonglake reciprocal transplants induce convergent expression of immune genes in threespine stickleback. Molecular Ecology, 24, 4629–4646. Svanb€ack, R., & Schluter, D. (2012). Niche specialization influences adaptive phenotypic plasticity in the threespine stickleback. The American Naturalist, 180, 50–59. Thompson, C. E., Taylor, E. B., & McPhail, J. D. (1997). Parallel evolution of lake-stream pairs of threespine sticklebacks (Gasterosteus) inferred from mitochondrial DNA variation. Evolution, 51, 1955–1965. Todd, E. V., Black, M. A., & Gemmell, N. J. (2016). The power and promise of RNA-seq in ecology and evolution. Molecular Ecology, 25, 1224–1241. Van Buskirk, J. (2002). A comparative test of the adaptive plasticity hypothesis: Relationships between habitat and phenotype in anuran larvae. The American Naturalist, 160, 87–102. Velotta, J. P., Wegrzyn, J. L., Ginzburg, S., Kang, L., Czesny, S., O’Neill, R. J., . . . Schultz, E. T. (2017). Transcriptomic imprints of adaptation to fresh water: Parallel evolution of osmoregulatory gene expression in the Alewife. Molecular Ecology, 26, 831–848. Wang, G., Yang, E., Smith, K. J., Zeng, Y., Ji, G., Connon, R., . . . Cai, J. J. (2014). Gene expression responses of threespine stickleback to salinity: Implications for salt-sensitive hypertension. Frontiers in Genetics, 5, 312. Weber, J. N., Bradburd, G. S., Stuart, Y. E., Stutz, W. E., & Bolnick, D. I. (2017). Partitioning the effects of isolation by distance, environment, and physical barriers on genomic divergence between parapatric threespine stickleback. Evolution, 71, 342–356. Weber, J. N., Steinel, N. C., Shim, K. C., & Bolnick, D. I. (2017). Recent evolution of cestode growth suppression by threespine stickleback. PNAS, 114, 6575–6580. Wegner, K. M., Kalbe, M., Kurtz, J., Reusch, T. B. H., & Milinski, M. (2003). Parasite selection for immunogenetic optimality. Science, 301, 1343. Wegner, K. M., Kalbe, M., Rauch, G., Kurtz, J., Schaschl, H., & Reusch, T. B. (2006). Genetic variation in MHC class II expression and interactions with MHC sequence polymorphism in three-spined sticklebacks. Molecular Ecology, 15, 1153–1164. West-Eberhard, M. J. (2003). Developmental plasticity and evolution. New York, NY: Oxford University Press. Whitehead, A., Roach, J. L., Zhang, S., & Galvez, F. (2011). Genomic mechanisms of evolved physiological plasticity in killifish distributed along an environmental salinity gradient. Proceedings of the National Academy of Sciences, 108, 6193–6198. Wright, R. M., Aglyamova, G. V., Meyer, E., & Matz, M. V. (2015). Gene expression associated with white syndromes in a reef building coral, Acropora hyacinthus. BMC Genomics, 16, 371.

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

How to cite this article: Lohman BK, Stutz WE, Bolnick DI. Gene expression stasis and plasticity following migration into a foreign environment. Mol Ecol. 2017;26:4657–4670. https://doi.org/10.1111/mec.14234