In absence of local adaptation, plasticity and spatially varying selection rule: a view from genomic reaction norms in a panmictic species (Anguilla rostrata)
BMC Genomics volume 15, Article number: 403 (2014)
American eel (Anguilla rostrata) is one of the few species for which panmixia has been demonstrated at the scale of the entire species. As such, the development of long term local adaptation is impossible. However, both plasticity and spatially varying selection have been invoked in explaining how American eel may cope with an unusual broad scope of environmental conditions. Here, we address this question through transcriptomic analyses and genomic reaction norms of eels from two geographic origins reared in controlled environments.
The null hypothesis of no difference in gene expression between eels from the two origins was rejected. Many unique transcripts and two out of seven gene clusters showed significant difference in expression, both at time of capture and after three months of common rearing. Differences in expression were observed at numerous genes representing many functional groups when comparing eels from a same origin reared under different salinity conditions. Plastic response to different rearing conditions varied among gene clusters with three clusters showing significant origin-environment interactions translating into differential genomic norms of reaction. Most genes and functional categories showing differences between origins were previously shown to be differentially expressed in a study comparing transcription profiles between adult European eels acclimated to different salinities.
These results emphasize that while plasticity in expression may be important, there is also a role for local genetic (and/or epigenetic) differences in explaining differences in gene expression between eels from different geographic origins. Such differences match those reported in genetically distinct populations in other fishes, both in terms of the proportion of genes that are differentially expressed and the diversity of biological functions involved. We thus propose that genetic differences between glass eels of different origins caused by spatially varying selection due to local environmental conditions translates into transcriptomic differences (including different genomic norms of reaction) which may in turn explain part of the phenotypic variance observed between different habitats colonized by eels.
Similar to other members of the freshwater eels of the genus Anguilla, the American eel (Anguilla rostrata, Lesueur 1817) represents one of the most scientifically intriguing cases of life history among fishes. Our recent findings unambiguously confirmed that this temperate euryhaline semelparous species comprises a single genetically panmictic population . As for the European eel, American eel exhibits extreme inter-individual phenotypic variance among teleost fishes during their so-called yellow eel (juvenile) phase in terms of growth rate, length at maturity, morphology, sex ratio, fecundity, and other life history traits as well as pronounced regional variation in recruitment across the highly heterogeneous environments occupied by the postlarval up to adult stages [2–12].
It has been hypothesized that phenotypic variation observed across the species’ distribution range is associated with differential mortality caused by spatially varying selection associated to environmental heterogeneity [12, 13]. Differential mortality is thus causing local genetic differences during a single generation and/or genetically-based intra-specific variation in migratory behavior and alternative habitat use . However, this does not exclude a role for plasticity in coping with heterogeneous environments [15, 16].
Understanding the relative importance of genetic versus plasticity factors represents a major challenge, and a first attempt to decipher both, we previously compared growth patterns in controlled conditions between young eels from two geographic origins . This study revealed both an origin (genetic) and environmental (plastic) effect on growth, and suggested a possible origin-environment interaction. The relative roles of these factors on patterns of phenotypic variation in eels and the mechanistic nature of such variation can be investigated by means of gene expression analyses. Thus, although gene transcripts themselves represent discrete phenotypes, they may also result in externally visible (e.g. size, color) or more cryptic (e.g. physiological tolerance to environmental conditions, age at reproduction) phenotypes [7, 17, 18]. Phenotypic variance is ultimately controlled by gene expression, which is modulated by the interactions between the genotype and the environment. In absence of any local genetic differences between locations, and given panmixia, one expects that eels from different origins reared in similar environments should present similar gene expression profiles. By contrast, observing differences of transcription profiles between eels from different origins or different genomic norms of reaction across different rearing environments would suggest a possible role for genetic differences between eels from different origins, possibly resulting from spatially varying selection. Admittedly, epigenetic effects could also be involved in such a case. However, these can hardly be investigated in non-model species given the current state of the knowledge and tools at hands, although the role of epigenetics in controlling phenotypic variance in eel would obviously be of major interest in subsequent steps . Finally, a better understanding of the underlying causes of phenotypic variation in the American eel across different environments is of outmost applied importance as this may impact on management decisions (e.g. translocations, defining management units) made towards improving its conservation [20, 21].
The goal of this study was to perform a transcriptomic study by means of cDNA microarrays in order to assess the relative role of geographic origin, environment (salinity) and their interaction on transcription profiles observed both in natural and controlled conditions. Our results revealed that all three factors play a role in explaining variation of gene expression in American eel, thus indicating that phenotypic variation observed in the species is mainly but not purely plastic and therefore also involves genetically (and/or epigenetic) based differences between eels from different locations.
This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the Canadian Council on Animal Care in science and all efforts were made to optimize environmental enrichment and avoid suffering. The protocol was approved by the Committee on the Good Care of Animal Experiments of Université Laval (Permit Number: 2007-165-1).
Sampling and rearing conditions
We previously described the details of sampling and rearing conditions in Côté et al. . In brief, glass eels were captured at two sampling locations in spring 2007: Mira river, 45°56’N 60°07’W, hereafter MR and Grande-Rivière-Blanche, 48°78’N 67°69’W, thereafter GRB. Glass eels correspond to the young life history stage at which they colonize coastal and continental waters following a 12 to 18 month oceanic migration from their spawning ground in the Sargasso Sea toward coastal or freshwater habitats .
Eels of both origins are characterized by contrasted life history in their natural settings. GRB represents eels from the St. Lawrence R. freshwater system which produces exclusively very large females  while MR represents a brackish environment were both males and small females occur in variable proportions and variable adult size (Dr. Martha Jones, Cape Breton University unpublished data). Glass eels of similar sizes were collected at the mouth of each river as soon as they approached the coast.
A sample of glass eels for each site was fixed for 24 h at 4°C in a saturated salt buffer immediately after collection in the field (T0) to preserve RNA and DNA integrity. The buffer was replaced within 24 h to a week and samples were stored at room temperature. The salt saturated buffer consisted of the following for one liter saturated salt solution; 40 mL of 0.5 M EDTA, 25 mL 1 M Sodium citrate, 700 g Ammonium sulfate and 935 mL ultrapure water, stirred on low heat combine until the Ammonium sulfate dissolves completely. Following cooling, pH was adjusted to 5.2 with H2SO4. Glass eel from both sites and captured at different time points were treated exactly in the same way, that is they were all randomly distributed into individual tanks after a same acclimation period of 12 days. This consisted in maintaining the freshwater condition at 3-5% and for the brackish water treatment, salinity was raised gradually from 3% to 22% with increment of approximately 2% per day. For each origin, glass eels were distributed and reared in 25 L aquaria at the Laboratoire de Recherche en Sciences Aquatiques (LARSA) at Université Laval for three months in two different salinity conditions: fresh (FW: 2–3 ppt) or brackish (BW: 20–22 ppt) water. At the end of this period, three samples per aquarium, for a total of 12 eels (hereafter named elvers) per origin and salinity conditions (total of 48 eels) were randomly collected. Individuals were sedated with 10% eugenol, measured for both total body length (LT ± 1 mm) and wet mass (W ± 0.1 mg), and frozen in liquid nitrogen.
High quality RNA preparation
Extraction, conservation and retro-transcription of mRNA were processed as follows. Each sample, consisting of whole body were ground in extraction buffer (TRIzol® reagent, Invitrogen) supplemented with β-mercaptoethanol (SIGMA-ALDRICH) and RNAse inhibitors (Ambion, Invitrogen) in 2 mL RNase free Eppendorfs (Qiagen). Total RNA was extracted using PureLink™ Micro-to-Midi Total RNA purification kit (Life Technologies). Residual genomic DNA was eliminated using DNAse for more accurate measurement of total RNA and more efficient cDNA synthesis. The DNAse (Promega) was inactivated properly prior to retro-transcription (RT) through phenol:chloroform:isoamyl purification. Total RNA was re-suspended in ultrapure water and measured with a NanoDrop 2000c spectrophotometer (Thermo Scientific, DE) and quality was checked using the Experion automated electrophoresis system (Biorad, CA, USA). An amount of 15 μg total RNA was used per RT reactions and fluorophore tagging was completed prior to the hybridisation on the microarray slide (Genisphere Array 350 3DNA™ and SuperScript II retro-transcriptase Invitrogen’ kit).
cDNA microarray description, hybridisation and analysis
The cDNA microarray used in this study was developed by Kalujnaia et al.  for European eel (Anguilla anguilla) and was successfully used in a previous study on A. rostrata. This array comprises 6144 expressed sequenced tags (ESTs) spotted in triplicate on glass silicate matrix lamella (GAPS II: Corning Inc, NY) by the service for microarray printing of the Biological Sciences School at Liverpool University. Initial ESTs were extracted from various tissues including brain, gills, guts and kidney.
Experimental design, data acquisition and analyses
The experimental design consisted in paired comparisons where two individuals were hybridized simultaneously on each array using two fluorophores (cyanine3 and cyanine5). From the initial 12 individuals randomly sampled per group and later hybridized, 10 were retained for the analysis (total of 60 individuals) for standardization in terms of both numbers of samples and the quality of the hybridization signal. Glass eels correspond to samples analyzed at T0 and 10 samples from each location were compared for the final analysis, half with each different fluorophore (Cy3, Cy5). For the analysis after three months of common rearing (T3), each sample (now named elvers) was used for two comparisons in order to assess the relative role of origin, salinity conditions, and their interactions. Ten samples (five for each fluorophore) for each salinity condition and for each origin were compared.
Hybridized microarrays were scanned using the ScanArray Express equipment (Packard Biosciences). Images were scanned at a 10 μm resolution and saved as *TIF files. Signal intensity (532 nm for Cy3 and 635 nm for Cy5) was extracted using QuantArray® (Perkin Elmer Life Sciences). Outlier spots and background noise were subtracted using the «Row Average Imputer» function and missing data were calculated through the «K-nearest neighbours» algorithm, both included in the R SAM software . Once the background noise was subtracted, the average intensity of each transcript (gene, cDNA, in triplicate) was divided by the average intensity of its correspondent fluorophore (channel). A given gene was retained for further analyses if its intensity exceeded twice the average intensity plus twice the standard deviation of the negative control (empty spots) found on the same array with the same fluorophore. Base-2 logarithm was used to normalise the distribution of the average values retained, which were then normalized with the regional lowess method of the R/MANOVA package.
To detect gene expression differences among the tested groups, data were analysed using a mixed model ANOVA  and the R/MAANOVA package . We tested for origin effects with the following model at T0 (Yijkl = μ + Ai + Dj + Ok + ϵijkl) and the following model for the experiments at T3 (Yijkl = μ + Ai + Dj + Ok + Sl + ϵijkl) where A: array; D: Dye; O: Origin; S: Salinity. The "Array" and "Sample" terms are treated as random whereas "Dye" and "Origin" are fixed terms. A permutation-based F-test (Fs statistic, with 1000 sample ID permutations) was performed and maximum likelihood was used to solve the mixed model equations . We tested for origin type effects with the ANOVA model and used the p-values to determine the significance of differential expression. In order to limit type I error, a False Discovery Rate correction (FDR = 0.05 and 0.10) was applied within the R/MAANOVA package to account for multiple testing. F tests with 1000 permutation and maximum-likelihood restriction were used to solve the mixed model.
In order to reduce type II errors, a hierarchical clustering analysis using the complete-linkage method was performed on the entire raw fold-changes dataset. This method creates compact classes of very similar observations in terms of co-variation in expression. The number of classes to be considered was established using the cubic clustering criterion  implemented in SAS 9.2. This revealed an optimized grouping into seven gene clusters (see Results). Then, the GLM procedure in SAS was applied to test for overall differences among clusters followed by a K-means analysis to identify which clusters differed among themselves (XLSTAT, 1995). The reaction norms for the different clusters were based on gene expression mean fold change difference between any two groups compared. Here, the group which was the least expressed relative to the other was set on the zero base line, for the purpose of the comparison in relative terms.
Transcript sequences were annotated using the software «Blast2GO PRO» . Sequences were uploaded in Blast2GO and analyses were run as recommended by the developer. First, the loaded sequences were run through blastx program on non-redundant database (nr Blast DB from QBlast-NCBI), 1.0E-3 Blast ExpectValue, with the maximum of 20 Blast Hits, word size of 3, low complexity filter, HSP length cutoff of 33. Secondly, the GO-Mapping sequences step was performed to link BLAST Hits to the functional information stored in the Gene Ontology database. Blast2GO uses different public resources provided by the NCBI, PIR, and GO to link the different protein IDs to the information stored in the GO database. All annotations are associated to Evidence Code which provides information about quality of this functional assignment (Blast2GO support information). The third step consisted in the Annotation default parameters which were set as followed: 1.0 E-6 E-Value-Hit-Filter, 55 of Annotation CutOff, five of GO Weight and no Hsp-Hit Coverage CutOff. Annotation Expander (ANNEX) was run to use an additional Gene Ontology structure which suggested new biological processes and cellular component annotations, based on the genes’ existing ‘’Molecular-Function” annotations (Ref.: ANNEX by GOAT, the Gene Ontology Annotation Toolbox, http://goat.man.ac.uk/ Simen Myhre and Henrik Tveit). The fourth step was the Enzyme Code Mapping and Kegg Pathway information, followed by InterProScan. Once all annotation steps were completed, GO and InterProScan informations were merged before running GO-slim.
Gene annotations and clustering
Genes that significantly differed in expression in the various comparisons (see below) were identified through homology in translated proteins in non-redundant NCBI, non-redundant (nr) and GO mapping. Since these putative proteins have been mainly characterized in a few model organisms, and only represents partial information in cDNA microarrays, we considered them as part of “gene families” rather than precise unigenes and avoided over-interpretation about their biological function as this must await further functional characterization in terms of ecological annotation . The proportion of annotated gene families varied between 28% and 68% (mean = 47%) depending on comparisons (see Results below). As hypothesized by Kalujnaia et al. , the high number of unidentified sequences obtained probably reflects the large number of 3’-untranslated region fragments that were spotted on the microarray. These gene families were then classified into eight functional categories previously described in Kalujnaia et al.  in their transcriptomic study of osmoregulation and development regulation in European eel (A. anguilla). These functional categories are: Cell protection/immunity, Detoxification, Energy metabolism/respiration, Growth/differentiation/development, Membrane transporters/carrier proteins, Signal transduction, Structural/junctional complex, and Transcription/translation. Table 1 presents the list of differentiated gene families within each of these functional categories, as well as their representation across different comparisons based on single transcript analysis and in each cluster based for the hierarchical clustering analysis (details below). Gene families that have previously been reported in other relevant studies on the transcriptomics of osmoregulation in fish are also indicated, with an emphasis on the study of Kalujnaia et al.  on European eel.
The hierarchical clustering analysis, performed on all genes detected that were common to all eel groups compared, was resolved by seven distinct clusters of genes that co-varied in expression (Table 2). The total number of transcripts per cluster varied between 183 (cluster 1) and 8 (cluster 7) whereas the number of annotated transcripts varied between 93 (cluster 3) and 4 (cluster 7). Annotated transcripts were assigned to different gene families that represented between one to eight functional categories depending on clusters. Thus, most gene families were shared between two clusters or more but one private gene family (detoxification) was observed for cluster 3. Details on the gene families belonging to different functional categories and represented in each cluster are presented in Table 1. These seven clusters were used for subsequent analyses, as detailed below.
Differentiation between glass eels at time of capture
The number of uniquely expressed transcripts varied between 1400 and 2179 depending on comparisons (Table 3). The null hypothesis of no difference of gene expression between glass eels from the two origins at the time of capture in natural conditions was rejected whereby 883 unique transcripts showed significant difference in expression between individuals from GRB and MR. Of these, 57% were over-expressed in MR relative to GRB. The number of transcripts passing both FDR thresholds were comparable to that observed based on p-value. A total of 154 unique transcripts were annotated and these were assigned to different gene families which represented all eight functional categories (Tables 1 and 4). The three functional categories that were the most represented in terms of absolute number of gene families were: Cell protection/immunity, Energy metabolism/respiration, and Growth/Differentiation/Development.
As for differentiation at individual transcripts, the null hypothesis of no difference in gene clusters that co-varied in expression was rejected. Figure 1a illustrates mean fold changes (in Log2) for transcripts included in each of the seven clusters and provides a general picture of how these differed between origins. The GLM procedure revealed a highly significant difference between mean expression level among the seven clusters (F = 123.69, F < 0.0001). K-means test indicated that the seven clusters formed three distinct groups (A:1,2; B:3,7: C4, D:5–6). The most salient result of this comparison was the pronounced down-regulation of clusters 3 and 7 in MR relative to GRB glass eels by approximately 60%. All eight functional categories are represented in cluster 3, among which three are shared with cluster 7 (Energy metabolism/respiration; Growth/differentiation/development; Membrane transporters/carrier proteins) (Table 2). Two annotated gene families were common to both clusters (Keratin and Apolipoproteins) (Table 1), and both were previously reported to show differences between European silver eels acclimated to different salinity conditions . Other gene families belonging to cluster 3 that were both represented by several transcripts and that were also reported by Kalujnaia et al.  were Cytochrome c oxydase, Maltase-glucoamylase, Claudin, Cathepsin, Ictacalcin, and Ribosomal protein (Table 1).
Differentiation between eels from distinct origins reared in fresh and brackish water
Significant differences between the two origins were maintained after three months of rearing in common environment but were less pronounced than at time of sampling in natural conditions. Thus, 259 and 422 unique transcripts showed differences in expression between origins when reared in BW and FW, respectively (Table 3). Of these, 55% were significantly over-expressed in GRB relative to MR in BW whereas this was the case for 54% in FW. However, the number of significant transcripts passing both FDR thresholds was substantially smaller, particularly in the BW comparison. Overall then, there was a trend towards less pronounced gene expression between origins in BW relative to freshwater. In the BW comparison, 67 of the differentiated transcripts were annotated and these represented seven of the eight functional categories (Table 4). In the FW comparison, 96 of the differentiated transcripts were annotated and these also represented seven of the eight functional categories. As in T0, the three functional categories that were the most represented in both BW and FW comparisons in terms of absolute number of gene families were: Cell protection/immunity, energy metabolism/respiration, and Growth/Differentiation/Development. The list of gene families that differed between origins in both environments and that were also reported by Kalujnaia et al.  was comparable to those observed at T0, comprising: Aminopeptidase, Cytochrome c oxydase, Maltase-glucoamylase, Claudin, Keratin, and Apolipoproteins, Ribosomal protein (Table 1).
We then tested for differentiation between origins depending on the rearing environment after three months (T3). Figure 1b illustrates mean fold change differences between MR and GRB for each of the seven clusters when reared in BW whereas Figure 1c represents the same comparison in FW. The GLM procedure revealed a highly significant difference between mean expression level among the seven clusters, both in BW and FW environment (respectively F = 88.73 and F = 85.16, p < 0.0001). The K-means test indicated that the behavior of some of the clusters varied between rearing conditions whereas others were similar. Thus, clusters 1 and 3 showed no clear difference in expression between origins both in BW and FW environments whereas clusters 4 and 7 were up-regulated in MR for both environments but the fold change for both clusters tended to be higher in FW than BW. Thus, there was a limited plastic response in expression for these genes. In contrary, the other three clusters showed a contrasted pattern of differentiation depending on the environment. For BW, cluster 2 was down-regulated in MR whereas no clear difference between origins was observed in FW. Conversely, cluster 6 was down-regulated for FW in MR whereas no clear difference in expression was observed for BW. Cluster 5 was down-regulated for BW and up-regulated in FW for MR. In addition to functional categories comprised in cluster 7 and defined above, genes that showed up-regulation for both environments in MR belonged to the Cell protection/immunity category (Table 2). This category was represented only by Aminopeptidase which was also reported by Kalujnaia et al.  (Table 1). Genes of cluster 5 presented the most pronounced pattern of plasticity, and all belonged to the Cell protection/immunity functional category represented by three gene families: two reported by Kalujnaia et al.  Aminopeptidase, Mucin, as well as Nattectin (Table 1). Genes belonging to cluster 6 were also plastic and represented four functional groups: Growth/differentiation/development, Membrane transporter/carrier proteins, Signal transduction and Transcription/translation (Table 2). These comprised four gene families that were all reported by Kalujnaia et al. : S100 calcium binding protein, Apolipoproteins, Interferon-induced gene 2, and Ribosomal protein (Table 1).
Another way to look at environmental effect on the plasticity of gene expression is to compare the extent of differentiation between FW and BW within a same origin. Thus, 32% and 28% of detected transcripts showed differences in expression between salinity conditions for MR and GRB, respectively (Table 3). Of these, 50% were over-expressed in BW relative to FW for MR whereas this was the case for 53% for GRB. The number of significant transcripts passing both FDR thresholds was important although lower to those observed with p-values. Thus, at T3, there were more differences in expression between environments than between origins based on absolute number of differentially expressed transcripts. The number of differentiated transcripts between environments was comparable for both origins and these represented all eight functional categories (Table 4). There were numerous gene families that differed in expression between environments for both origins and that were also reported to show differences between European silver eels acclimated to different salinity conditions . Those represented by the largest number of transcripts were: Aminopeptidase, Lectin (β-galactoside-binding), Cytochrome c oxydase, Maltase-glucoamylase, Claudin, Keratin, Apolipoproteins, and Ribosomal protein (Table 1). Na.K-ATPase (one transcript) also differed in expression for both origins. It is also noteworthy that seven transcripts of the KDELR1 (not reported in previous studies) gene family differed in expression for both origins. Finally, a few gene families showed contrasted patterns between origins. In particular, many more transcripts of Cytochrome c oxydase were differentially expressed between environments in GRB vs. MR (17 vs. 3). Also, five transcripts of Solute carrier 15 and 22 transcripts of Chemotaxin differed in expression in GRB vs. none in MR. Finally; there were twice as many transcripts of Apolipoproteins differentially expressed in MR vs. GRB (12 vs. 6).
We then compared the plastic response of the seven gene clusters between environments for eels from the two origins. Figure 1d illustrates mean fold change differences between environments for MR and Figure 1e represents the same comparison for GRB. The GLM procedure revealed a highly significant difference between mean expression level among the seven clusters, both for MR (F = 155.10, F < 0.0001) and GRB (F = 196.70, F < 0.0001). Again, the K-means test indicated that the plastic response to different rearing conditions varied among clusters, with four different cluster groups observed in MR and all seven gene clusters showing a different plastic response in GRB. Variation between rearing environments for clusters 1 to 6 were marginally correlated between MR and GRB (r2 = 0.646, F = 7.293 and p = 0.054), suggesting similar directionality in plastic response between origins to environmental conditions for these genes. In contrast, cluster 7 showed a pronounced gene-by-environment interaction whereby genes of this cluster were highly overexpressed in FW relative to BW water for GRB (fold change = 1.49) whereas they slightly tended to be overexpressed in brackish water in MR (mean fold change = 1.07). These differences translated into contrasted genomic reaction norms between origins for cluster 7 (Figure 2). Variable patterns of genomic reaction norms between origins were observed among the other six clusters. Cluster 1 showed no plastic response for both origins, cluster 2 showed no plastic response in MR but was overexpressed in BW relative to FW in GRB, translating origin-environment interaction. Clusters 3 and 6 were overexpressed in FW for both origins, whereas clusters 4 and 5 also showed a pronounced plastic response, both being overexpressed in BW for both origins. For cluster 4, however, the plastic response was more pronounced in MR than GRB, thus also translating origin-environment interactions. Thus, gene families that showed the most pronounced differences in genomic reaction norms between origins belonged to clusters 2, 4 and 7, the most important being: Lectin, Pentraxin, Maltase-glucoamylase, Elongation factor 1-alpha, Profilin, Chemotaxin, Aminopeptidase, Enolase, Keratin, Apolipoproteins and Ribosomal protein, all reported by Kalujnaia et al.  except for Chemotaxin (Table 1).
Since panmixia precludes any possibility for long-term local adaptations, a next step towards understanding how American eel can cope with highly heterogeneous environments is to identify factors explaining the pronounced local variations in life history as well as physiological and ecological traits observed in this species throughout its range of distribution . In order to address this issue from a gene expression perspective, we assessed the relative role of geographic origin, environment (salinity) and their interaction on transcription profiles observed both in natural and controlled conditions. Our results revealed that all three factors play a role in explaining variation of gene expression in American eel.
Transcriptomic divergence in natural conditions
Our results rejected the null hypothesis of no differences in gene expression between glass eels collected at same sizes and developmental stage, at the same time of arrival at different river mouths just prior to entering freshwater, thus avoiding differential effect of time spent in freshwater before experiments . Nearly 40% of the 2179 unique transcripts that were detected showed significant difference in expression between individuals from GRB and MR. We cannot rule out the possibility that transcriptomic differences between sites are due to differences in local environmental conditions at time of sampling, namely temperature. Thus, mean sea surface temperatures are higher along the coast of Nova Scotia than in the St. Lawrence estuary (http://www.class.ncdc.noaa.gov data base SST14NA). However, gene expression generally increases with temperature in fishes [32, 33] whereas we observed a trend towards overexpression in GRB relative to MR. Also, both samples were taken in similar brackish water just at the mouth of both rivers. Thus, different temperatures and salinities encountered by glass eels at time of capture may not be the most parsimonious explanation for the observation of differential expression between origins. Alternatively, differential selective pressures associated with variable environmental conditions being encountered during marine dispersal can result in genetic differences between glass eels from different locations through the process of spatially varying selection as demonstrated by Gagnaire et al. (2012) . Notwithstanding a role for environmental plastic response in explaining differential expression between glass eels of different geographic origins (see below), this first result along with evidence for local genetic differences at coding genes points towards a role for some genetic differences between glass eels from different origins. Admittedly, epigenetic effects could also be invoked in explaining the observed differential pattern of expression between origins and cannot be refuted , although there is at present very little empirical evidence for this in nature . Also, given the impossibility to date to routinely perform artificial reproduction and rear American eel larvae beyond the age of about 10 days , it was not possible to perform crosses and rear F2 progeny to rigorously rule out possible environmental and maternal carry-over effects that could impact on the observed patterns of gene expression.
Transcriptomic divergence in controlled conditions: comparisons with other studies
The experiment conducted over a three-month period revealed that transcriptomic differences were maintained in controlled conditions. Further evidence for differential response between young eels from GRB and MR was provided by the origin-environment interactions observed for some of the gene clusters. Overall, general patterns of gene expression between distinct origins in controlled conditions corroborate those we observed previously for growth between glass eels from these same two locations . Moreover, growth differences between origins associated with sex-specific effects were also maintained after three years of common rearing and were similar to those observed in natural conditions (Côté et al. in revision).
The scale of differences observed here between eels of both origins is comparable and even more pronounced to what has been reported in other fish studies comparing the transcriptome of genetically distinct populations within a same species. In a comparison between genetically distinct resident and anadromous brook charr (Salvelinus fontinalis) reared in common conditions, Boulet et al. (2012)  observed that gene expression in the gills differed for 4% of the genes, a value comparable to what we observed after 3 months of common rearing. In a study of divergence in gene regulation at young life history stages of whitefish (Coregonus sp.), Nolte et al. (2009)  found that 8.6% of genes differed in expressions between dwarf and normal whitefish populations. In a study comparing the ecological transcriptomics of lake-type and riverine sockeye salmon (Oncorhynchus nerka), Pavey et al. (2011)  found that about 1% transcripts differed between these locally adapted populations to distinct habitat types. In a comparative analysis of gene transcription between genetically distinct populations of the Atlantic killifish (Fundulus heteroclitus), Whitehead et al. (2011)  found that 8.2% of the 6800 genes analyzed differed in expression between populations.
In summary, results to date suggest that local genetic differences between glass eels of different origins are at least partly due to spatially varying selection translating into transcriptomic differences (including different genomic norms of reaction), which could partly explain the phenotypic variance observed between eels from different habitats. This offers a plausible working hypothesis towards solving the apparent paradox between observations of important regional variation in life-history traits, demography, and ecology on the one hand and on the other hand, evidence for the existence of a single, randomly mating population in Atlantic eel.
Since growth is a predominant phenotypic trait that differs between eels from GRB and MR, one should expect that genes potentially involved in growth functions differ between them as well. Indeed, this was the case for some genes involved. Thus, some gene clusters (clusters 4, 7) showed either a net difference between origins when reared either in FW or BW or showed an origin effect on plastic response to environmental conditions (clusters 5, 6). Some of the genes super-families found in these clusters could potentially influence growth, for instance Cytochrome c oxidase (COX) which is involved in the respiratory chain , Keratin which is involved in critical cells functions such as migration, size, growth, transport and proliferation , and Profilin-2 which enhances actin growth . However, growth is a complex phenotypic trait with a highly polygenic basis such that it would be risky at best to attempt any further causal links between specific functions and variation in growth. Besides these super-families that could play a role in growth, many other functions differed in expression between eels from different origins, as reported also in other fish studies comparing genetically distinct locally adapted populations [30, 37–39].
Many of the gene families that showed significant differences in expression or expressed different genomic norms of reactions between geographic origins overlapped with genes that showed significant differences between adult silver European eels from a same location acclimated to different salinity conditions . Many of those genes do not have any obvious direct link to osmoregulatory functions, indicating that salinity conditions may impact on the expression of many biological functions, including cell protection, immunity, detoxification, energy metabolism, growth, differentiation and development, membrane transporters, signal transduction, structural and junctional complex as well as transcription and translation. Here, while we observed a plastic response to salinity conditions for many of the gene families represented (see discussion below), we also found that some of those genes showed a non-plastic response that differed between origins or differential plastic responses between origins.
Genomic reaction norms and possible evidence for adaptive plasticity
These results rule out the general belief that all variation in phenotypes and habitat use observed in eel is purely plastic [43, 44]. On the other hand, they also support an important role of phenotypic plasticity in explaining eel evolutionary success in coping with a very broad array of ecological conditions. . Thus, we observed more differences in expression between environments within origin than between origins for both rearing environments. However, the extent of plasticity to different salinity conditions was highly variable among genes. First, it is noteworthy that the overall pattern of plasticity in expression was relatively similar between origins whereby variation between clusters 1 to 6 was marginally correlated (Figure 1d and Figure 1e). However, different clusters showed different patterns of plasticity. The 183 genes belonging to cluster 1 did not show any variation in any of the comparison and as such, represent the group of genes that were the least plastic in expression. Clusters 3, 5 and 6 showed similar plastic responses in both origins. Finally, clusters 2, 4 and 7 revealed distinct reaction norms between GRB and MR. Differential plastic responses to environmental conditions among genes have commonly been reported and used to identify genes that may be the most important in coping with different environments [23, 32, 33, 45–47]. These studies typically identified a broad array of genes that vary in expression despite no a priori expectation based on functions. This was also the case here whereby we observed a plastic response in expression to salinity condition for genes belonging to many gene super-families and functional groups without any known a priori role in osmoregulation. Nevertheless, several individual genes potentially involved in osmoregulation were also identified (Table 3), such has prolactin, Na-K ATPase, and kinase C [23, 39, 48–50].
The different genomic reactions observed between eels from different origins could reflect an adaptive plastic response to the environment. This could be the case if they were causally linked with differential plastic responses in growth to these same salinity conditions we observed previously . This, however, needs to be rigorously investigated in future studies. Determining how plastic developmental changes that occur in response to environmental conditions (e.g. growth) are coordinated at the molecular (e.g. transcriptome level) represents a major challenge to be addressed in future studies by combining ecology with developmental biology .
Conclusions and perspectives
In summary, our results revealed that geographic origin, environment and the interaction between them play a role in explaining variation of gene expression in American eel. This indicates that phenotypic variation observed in American eel is not purely plastic but also involves genetically (and/or epigenetic) based differences between locations. Whether these differences are adaptive remains to be demonstrated. Yet, the extent and nature of differential gene expression between eels of different geographic origins compares to differences reported between genetically distinct and locally adapted populations in other species. Given the increasing availability of genomic resources for studying Anguilla species, including genome sequences [52, 53] large annotated EST libraries  and dense microarrays , it will be feasible to investigate into more details how plastic developmental changes that occur in response to environmental conditions are controlled at the molecular level, including the detection of epigenetic regulation mechanism , in link with genome-wide spatially varying selection caused by natural environmental heterogeneity , habitat degradation [22, 57], selective exploitation  or predation , as well investigating the molecular basis of behavioral traits that may be involved in habitat selection . Such knowledge is essential towards improving conservation strategies for this threatened species in Canada .
The data set supporting the results of this article is available from the Dryad Digital Repository: http://doi.org/10.5061/dryad.nq38s.
We are grateful to Guy Verreault (MRN) for his comments and field tips, to Eric Normandeau for laboratory and statistical assistance, Lucie Papillon, Laure Devine, Nathalie Brodeur and Hacène Tamdrari for field assistance, to Yvonne Carey from Atlantic Elver Fishery for gracefully providing dip nets, and to Mike Campbell from South Shore Trading Co. for providing samples and advice on elver rearing. The assistance of Serge Higgins and Jean-Christophe Therrien from LARSA was invaluable throughout the rearing experiment. We are also grateful to Laura Castonguay for proof-reading and improving the quality of language. This research was funded through the Department of Fisheries and Oceans Canada as well as a Discovery grant from the Natural Sciences and Engineering Research Council of Canada to LB. This is a contribution to the research program of Québec-Océan.
Grande Rivière Blanche
Time of collection on the field
Time of collection after three months of common rearing
Laboratoire de Recherche en Science Aquatique
Fresh water 2–3 ppt
Brackish water 20–22 ppt
Total body length
Expressed sequenced tags
- Cy3 and Cy5:
Cyanine 3 and cyanine 5 fluorophores
False discovery rate correction
Cytochrome c oxidase.
Côté CL, Gagnaire P-A, Bourret V, Verreault G, Castonguay M, Bernatchez L: Population genetics of the American eel (Anguilla rostrata): F(ST) = 0 and North Atlantic Oscillation effects on demographic fluctuations of a panmictic species. Mol Ecol. 2013, 22 (7): 1763-1776.
Pratt TC, Threader RW: Preliminary evaluation of a large-scale American eel conservation stocking experiment. North Am J Fisheries Manage. 2011, 31 (4): 619-627.
Tremblay V: Reproductive strategy of female American eels among five subpopulations in the St. Lawrence River Watershed. Eels at the edge: science, status, and conservation concerns. Edited by: Casselman JM, Cairns DK. 2009, Bethesda: American Fisheries Society, 85-102. vol. Symposium 58,
Jessop BM, Shiao JC, Iizuka Y, Tzeng WN: Migratory behaviour and habitat use by American eels Anguilla rostrata as revealed by otolith microchemistry. Mar Ecol Prog Ser. 2002, 233 (2): 217-229.
Jessop BM: Geographic effects on American eel (Anguilla rostrata) life history characteristics and strategies. Can J Fish Aquat Sci. 2010, 67 (2): 326-346.
Côté CL, Castonguay M, Verreault G, Bernatchez L: Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes. J Fish Biol. 2009, 74 (9): 1934-1948.
Wang CH, Tzeng WN: Interpretation of geographic variation in size of American eel Anguilla rostrata elvers on the Atlantic coast of North America using their life history and otolith ageing. Mar Ecol Prog Ser. 1998, 168 (1): 35-43.
Vladykov VD, Liew PKL: Sex of adult American eels (Anguilla rostrata) collected as elvers in two different streams along the eastern shore of Canada, and raised in the same freshwater pond in Ontario. Proceedings of the 1980 North American eel Conference. Ont Fishries Tech Rep Ser. 1982, 4: 88-93.
Ide C, De Schepper N, Christiaens J, Van Liefferinge C, Herrel A, Goemans G, Meire P, Belpaire C, Geeraerts C, Adriaens D: Bimodality in head shape in European eel. J Zool. 2011, 285 (2): 230-238.
Verreault G, Dumont P: An estimation of American eel escapement from the upper St. Lawrence River and Lake Ontario in 1996 and 1997. Biology, management, and protection of catadromous eels. Edited by: Dixon DA. 2003, Bethesda, Maryland: American Fisheries Society, Symposium, 243-251. 33,
Laflamme S, Côté C, Gagnaire P-A, Castonguay M, Bernatchez L: RNA/DNA ratios in American glass eels (Anguilla rostrata): evidence for latitudinal variation in physiological status and constraints to oceanic migration?. Ecol Evol. 2012, 2 (5): 875-884.
Gagnaire P-A, Normandeau E, Côté C, Hansen MM, Bernatchez L: The genetic consequences of spatially varying selection in the panmictic American eel (Anguilla rostrata). Genetics. 2012, 190 (2): 725-735.
Hedrick PW: Genetic polymorphism in heterogeneous environments: the age of genomics. Annu Rev Ecol Evol Syst. 2006, 37 (1): 67-93.
Edeline E: Adaptive phenotypic plasticity of eel diadromy. Mar Ecol Prog Ser. 2007, 341 (2): 229-232.
West-Eberhard MJ: Phenotypic accommodation: Adaptive innovation due to developmental plasticity, with or without genetic change. Integr Comp Biol. 2003, 43 (6): 970-970.
West-Eberhard MJ: Phenotypic accommodation: Adaptive innovation due to developmental plasticity. J Exp Zool Part B-Molecular Dev Evol. 2005, 304B (6): 610-618.
Wang CH, Tzeng WN: The timing of metamorphosis and growth rates of American and European eel leptocephali: a mechanism of larval segregative migration. Fish Res. 2000, 46 (2): 191-205.
Tzeng WN, Severin KP, Wang CH, Wickstrom H: Elemental composition of otoliths as a discriminator of life stage and growth habitat of the European eel. Anguilla anguilla Mar Freshwater Res. 2005, 56 (5): 629-635.
Lloyd VK, Brisson JA, Fitzpatrick KA, McEachern LA, Verhulst EC: The Epigenetics of Emerging and Nonmodel Organisms. Genet Res Int. 2012, 2012: 491204-
Allendorf FW, Hohenlohe PA, Luikart G: Genomics and the future of conservation genetics. Nat Rev Genet. 2010, 11 (10): 697-709.
Maes GE, Volckaert FAM: Challenges for genetic research in European eel management. Ices J Mar Sci. 2007, 64 (7): 1463-1471.
Castonguay M, Hodson PV, Couillard CM, Eckersley MJ, Dutil JD, Verreault G: Why Is Recruitment of the American Eel, Anguilla-Rostrata, Declining in the St-Lawrence-River and Gulf. Can J Fish Aquat Sci. 1994, 51 (2): 479-488.
Kalujnaia S, McWilliam IS, Zaguinaiko VA, Feilen AL, Nicholson J, Hazon N, Cutler CP, Cramb G: Transcriptomic approach to the study of osmoregulation in the European eel Anguilla anguilla. Physiol Genomics. 2007, 31 (3): 385-401.
Bernatchez L, St-Cyr J, Normandeau E, Maes GE, Als TD, Kalujnaia S, Cramb G, Castonguay M, Hansen MM: Differential timing of gene expression regulation between leptocephali of the two Anguilla eel species in the Sargasso Sea. Ecol Evol. 2011, 1 (4): 459-468.
Conesa A, Götz S, Miguel J, García-Gómez , Terol J, Talón M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005, 21 (18): 3674-3676.
Wolfinger RD, Gibson G, Wolfinger ED, Bennett L, Hamadeh H, Bushel P, Afshari C, Paules RS: Assessing gene significance from cDNA microarray expression data via mixed models. J Comput Biol. 2001, 8 (6): 625-637.
Kerr MK, Leiter EH, Picard L, CG A: Sources of Variation in Microarray Experiments. Computational and Statistical: Approaches to Genomics. Edited by: Zhang W, Shmulevich I. 2002, US: Springer, 37-47.
Cui X, Hwang JT, Qiu J, Blades NJ, Churchill GA: Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatistics. 2005, 6 (1): 59-65.
Milligan GW, Cooper MC: An examination of procedures for determining the number of clusters in a data set. Psychometrika. 1985, 50 (2): 159-178.
Pavey SA, Sutherland BJ, Leong J, Robb A, Schalburg KV, Hamon TR, Koop BF, Nielsen JL: Ecological transcriptomics of lake-type and riverine sockeye salmon (Oncorhynchus nerka). BMC Ecol. 2011, 2011 (1): 11-18.
Castonguay M, Hodson PV, Moriarty C, Drinkwater KF, Jessop BM: Is there a role of ocean environment in American and European eel decline?. Fish Oceanogr. 1994, 3 (3): 197-203.
Jeffries KM, Hinch SG, Sierocinski T, Clark TD, Eliason EJ, Donaldson MR, Li S, Pavlidis P, Miller KM: Consequences of high temperatures and premature mortality on the transcriptome and blood physiology of wild adult sockeye salmon (Oncorhynchus nerka). Ecol Evol. 2012, 2 (7): 1747-1764.
Smith S, Bernatchez L, Beheregaray LB: RNA-seq analysis reveals extensive transcriptional plasticity to temperature stress in a freshwater fish species. BMC Genomics. 2013, 14: 375-
Aubin-Horth N, Renn SC: Genomic reaction norms: using integrative biology to understand molecular mechanisms of phenotypic plasticity. Mol Ecol. 2009, 18 (18): 3763-3780.
Angers B, Castonguay E, Massicotte R: Environmentally induced phenotypes and DNA methylation: how to deal with unpredictable conditions until the next generation and after. Mol Ecol. 2010, 19 (6): 1283-1294.
Oliveira K, Hable WE: Artificial maturation, fertilization, and early development of the American eel (Anguilla rostrata). Can J Zool. 2010, 88 (11): 1121-1128.
Boulet M, Normandeau E, Bougas B, Audet C, Bernatchez L: Comparative transcriptomics of anadromous and resident brook charr (Salvelinus fontinalis) before their first salt water transition. Curr Zool. 2012, 58 (2): 155-167.
Nolte AW, Renaut S, Bernatchez L: Divergence in gene regulation at young life history stages of whitefish (Coregonus sp.) and the emergence of genomic isolation. BMC Evol Biol. 2009, 9 (1): 59-70.
Whitehead A, Roach J, Zhang S, Galvez F: Genomic mechanisms of evolved physiological plasticity in killifish distributed along an environmental salinity gradient. Proc Natl Acad Sci U S A. 2011, 108 (15): 6193-6198.
Horn D, Barrientos A: Mitochondrial copper metabolism and delivery to cytochrome c oxidase. Int Union Biochem Mol Biol Life. 2008, 60 (7): 421-430.
Bose A, Teh M-T, Hutchison IL, Wan H, Leigh IM, Waseem A: Two Mechanisms Regulate Keratin K15 Expression In Keratinocytes: Role of PKC/AP-1 and FOXM1 Mediated Signalling. Plos One. 2012, 7 (6): e38599-
Di Nardo A, Gareus R, Kwiatkowski D, Witke W: Alternative splicing of the mouse profilin II gene generates functionally different profilin isoforms. J Cell Sci. 2000, 113 (21): 3795-3802.
Daverat F, Limburg KE, Thibault I, Shiao JC, Dodson JJ, Caron FO, Tzeng WN, Iizuka Y, Wickstrom H: Phenotypic plasticity of habitat use by three temperate eel species, Anguilla anguilla. A-japonica A-rostrata Mar Ecol-Progr Ser. 2006, 308 (2): 231-241.
Capoccioni F, Lin DY, Iizuka Y, Tzeng W-N, Ciccotti E: Phenotypic plasticity in habitat use and growth of the European eel (Anguilla anguilla) in transitional waters in the Mediterranean area. Ecol Freshw Fish. 2014, 23 (1): 65-76.
Tzeng WN, Mail Y-HT, Han Y-S, Hsu C-C, Chang C-W, Lorenzo ED, Hsieh CH: Evaluation of Multi-Scale Climate Effects on Annual Recruitment Levels of the Japanese Eel, Anguilla japonica, to Taiwan. Plos One. 2012, 7 (2): e30805-
Evans ML, Bernatchez L: Oxidative phosphorylation gene transcription in whitefish species pairs reveals patterns of parallel and non-parallel physiological divergence. J Evol Biol. 2012, 25 (9): 1823-1834.
Richards RC, Short CE, Driedzic WR, Ewart KV: Seasonal changes in hepatic gene expression reveal modulation of multiple processes in rainbow smelt (Osmerus mordax). Marine Biotechnol. 2010, 12 (6): 650-663.
Hoffman EA, Goodisman MAD: Gene expression and the evolution of phenotypic diversity in social wasps. BMC Biol. 2007, 5 (1): 23-31.
Wilson JM, Leitao A, Goncalves A, Ferreira C, Pinto C, Novo D, Pereira-Wilson C, Fonseca AV, dos Santos PNR, da Silva JC, Coimbra J: The effect of salinity on branchial ion transport protein expression in the euryhaline glass eel (Anguilla anguilla). Comp Biochem Physiol a-Mol Integr Physiol. 2005, 141 (3): S211-S211.
Wilson JM, Antunes JC, Bouça PD, Coimbra J: Osmoregulatory plasticity of the glass eel of Anguilla anguilla: freshwater entry and changes in branchial ion-transport protein expression. Can J Fisheries Aquat Sci. 2004, 61 (3): 432-442.
Gilbert SF: Mechanisms for the environmental regulation of gene expression: Ecological aspects of animal development. J Biosci. 2005, 30 (1): 65-74.
Henkel CV, Dirks RP, Wijze DL, Minegishi Y, Aoyama J, Jansen HJ, Turner B, Knudsen B, Bundgaard M, Hvam KL, Boetzer M, Pirovano W, Weltzien FA, Dufour S, Tsukamoto K, Spaink HP, van den Thillart GE: First draft genome sequence of the Japanese eel, Anguilla japonica. Gene. 2012, 511 (2): 195-201.
Henkel CV, Burgerhout E, de Wijze DL, Dirks RP, Minegishi Y, Jansen HJ, Spaink HP, Dufour S, Weltzien FA, Tsukamoto K, van den Thillart GE: Primitive Duplicate Hox Clusters in the European Eel's Genome. PlosOne. 2012, 7 (2): e32231-
Coppe A, Pujolar J, Maes G, Larsen P, Hansen M, Bernatchez L, Zane L, Bortoluzzi S: Sequencing, de novo annotation and analysis of the first Anguilla anguilla transcriptome: EeelBase opens new perspectives for the study of the critically endangered European eel. BMC Genomics. 2010, 11 (11): 635-647.
Pujolar J, Marino I, Milan M, Coppe A, Maes G, Capoccioni F, Ciccotti E, Bervoets L, Covaci A, Belpaire C, Cramb G, Patarnello T, Bargelloni L, Bortoluzzi S, Zane L: Surviving in a toxic world: transcriptomics and gene expression profiling in response to environmental pollution in the critically endangered European eel. BMC Genomics. 2012, 13 (11): 507-527.
Bossdorf O, Richards CL, Pigliucci M: Epigenetics for ecologists. Ecol Lett. 2008, 11 (2): 106-115.
Couillard CM, Hodson PV, Castonguay M: Correlations between pathological changes and chemical contamination in American eels, Anguilla rostrata, from the St. Lawrence River. Can J Fisheries Aquat Sci. 1997, 54 (8): 1916-1927.
Bevacqua D, Capoccioni F, Melià P, Vincenzi S, Pujolar JM, Giulio ADL, Eleonora C: Fishery-Induced Selection for Slow Somatic Growth in European Eel. PLoS One. 2012, 7 (5): e37622-
Béguer-Pon M, Benchetrit J, Castonguay M, Aarestrup K, Campana SE, Michael JWS, Julian JD: Shark Predation on Migrating Adult American Eels (Anguilla rostrata) in the Gulf of St. Lawrence. Plos One. 2012, 7 (10): e46830-
Edeline E, Dufour S, Elie P: Role of glass eel salinity preference in the control of habitat selection and growth plasticity in Anguilla anguilla. Mar Ecol Prog Ser. 2005, 304: 191-199.
COSEWIC: COSEWIC assessment and status report on the American Eel Anguilla rostrata in Canada. Committee on the Status of Endangered Wildlife in Canada. Edited by: Tremblay V. 2012, Ottawa: Her Majestu the Queen in Right of Canada, 109-Catalogue No. CW69-14/458-2012E-PDF. xii + (http://www.registrelep-sararegistry.gc.ca/default_e.cfm) ISBN ISBN 978-1-100-20712-4
The data of this paper are original and no part of this manuscript has been published or submitted for publication elsewhere. The authors have no competing interests in this study.
CLC conceived and coordinated of the study, captured, reared and processed the fish specimens used in this study, carried out the molecular study, statistics and draft the manuscript. SKM and GC designed, annotated, and provided the microarray. MC and LB participated in all the steps of the conception the study, coordination and helped to draft the manuscript. All authors read and approved the final manuscript.
About this article
Cite this article
Côté, C.L., Castonguay, M., Kalujnaia, M.S. et al. In absence of local adaptation, plasticity and spatially varying selection rule: a view from genomic reaction norms in a panmictic species (Anguilla rostrata). BMC Genomics 15, 403 (2014). https://doi.org/10.1186/1471-2164-15-403