An integrated transcriptomic and comparative genomic analysis of differential gene expression in Arctic charr (Salvelinus alpinus) following seawater exposure

Migration between freshwater and seawater habitats is an important stage in the ontogeny of anadromous salmonids (Hoar, 1988). Freshwater and seawater impose adverse and contrasting types of osmoregulatory stress, such that maintenance of internal ion homeostasis is paramount for survival. The transition from freshwater to seawater requires osmoregulatory mechanisms to switch from a state of ion absorption (i.e. hyper-osmoregulation) to one of ion excretion (i.e. hypo-osmoregulation). Gill tissue epithelial cells play an integral role in this process (Evans et al., 2005; Marshall and Grosell, 2006).

Mitochondrion-rich cells (MRCs) form a complex with accessory cells and pavement cells in the epithelial layer of gill tissue (Evans et al., 2005; Marshall and Grosell, 2006). In seawater, this complex facilitates excretion of excess Cl⁻ and Na⁺ from blood to seawater. Cl⁻ excretion is achieved via the collective function of three interdependent membrane-bound ion transporters found in MRCs: Na⁺/K⁺/2Cl⁻, Na⁺/K⁺-ATPase and cystic fibrosis transmembrane conductance regulator. Na⁺ excretion occurs through leaky paracellular pores located between MRCs and accessory cells (Silva et al., 1977; Marshall and Grosell, 2006), where permeability seems to be dependent upon the claudin isoform located at the tight junction (Furuse et al., 2001; Van Itallie et al., 2006; Tipsmark et al., 2008). Elevated expression of these transporters occurs concomitantly with an increase in the size and quantity of MRCs, which is controlled by cortisol and growth hormone (McCormick, 2001; Marshall and Grosell, 2006). Despite the prominent role that MRCs have in hypo-osmoregulation, they make up less than 10% of gill tissue epithelium, whereas pavement cells cover over 90% (Wilson and Laurent, 2002; Evans et al., 2005).

Variable functions for the genes regulating ion transport in teleost gill tissues can be achieved as a result of the whole genome duplications that teleost fishes underwent in their evolutionary origins (Spring, 1997; Jaillon et al., 2004). Furthermore, as a result of tandem duplications and an additional whole-genome duplication in salmonids (Allendorf and Thorgaard, 1984), these species possess a repertoire of functionally distinct isoforms of genes that are components of the seawater MRC complex. Na⁺/K⁺-ATPase isoforms ATP1α1a and ATP1α1b exhibit reciprocal expression in freshwater and seawater. In seawater, ATP1α1b is up-regulated and ATP1α1a is down-regulated, while the opposite occurs in freshwater (Richards et al., 2003; Mackie et al., 2005; Bystriansky et al., 2006; Nilsen et al., 2007; Larsen et al., 2008; Madsen et al., 2009; McCormick et al., 2009). ATP1α1a and ATP1α1b are also expressed in separate MRCs (McCormick et al., 2009). There are two NKCC1 isoforms (i.e. NKCC1a and NKCC1b), but evidence suggests that only NKCC1a is significantly expressed in the gills of fishes in seawater (Cutler and Cramb, 2002; Tipsmark et al., 2002; Scott et al., 2004; Lorin-Nebel and Charmantier, 2007; Tse et al., 2006; Hiroi et al., 2008; Kang et al., 2010). In Atlantic salmon (Salmo salar), two CFTR isoforms are expressed in gill (i.e. CFTRI and CFTRII) (Chen et al., 2001; Mackie et al., 2005). Although both are up-regulated in seawater, temporal expression patterns are variable (Singer et al., 2002; Mackie et al., 2007; Aykanat et al., 2011). Five claudin isoforms have been detected in Atlantic salmon gill tissue: CLDN10E, CLDN27A, CLDN28A, CLDN28B and CLDN30. In seawater, CLDN10E mRNA is significantly up-regulated, while CLDN27A and CLDN30 mRNA are down-regulated (Tipsmark et al., 2008; Engelund et al., 2012). Overall, the pattern emerging seems to be an evolutionary divergence of osmoregulatory function towards distinct freshwater and seawater isoforms.

Large-scale transcriptomics analyses in a variety of teleosts have broadened our understanding of hypo-osmoregulation in fishes. Studies that assess differential gene expression in gill between fish in freshwater and seawater (Boutet et al., 2006; Kalujnaia et al., 2007; Evans and Somero, 2008; Tine et al., 2008; Evans et al., 2011) demonstrate that hypo-osmoregulation is influenced by diverse biological processes that extend beyond the mechanisms described by the seawater MRC complex. As the number of studies has increased, patterns have emerged suggesting that hypo-osmoregulation in different taxa is influenced by some of the same biological processes. These include cell-cycle regulation, cell metabolism, protein synthesis and regulation, cytoskeleton and structural proteins, and immune system response (Boutet et al., 2006; Kalujnaia et al., 2007; Evans and Somero, 2008; Tine et al., 2008). Notable in their absence are controlled laboratory experiments on salmonids, species for which salinity tolerance has important implications from both life-history and economic perspectives. Salmonids have been studied in their natural state during different life-history stages in their anadromous lifecycle, although such studies are far from complete. For example, Evans et al. (Evans et al., 2011) examined the hyper-osmoregulatory aspects of gill function by sampling sockeye salmon (Oncorhynchus nerka) on their return to freshwater natal spawning grounds. Because their work was based on fish migrating from seawater to freshwater (instead of freshwater to seawater), it may only reflect part of the story about differential expression (DE) in genes involved in seawater acclimation.

Large-scale gene expression analysis can also be used to gain insight into the genetic basis of salinity tolerance. Salmonids exhibit diverse salinity tolerance capacities, both within and among species (Hoar, 1988; Delabbio et al., 1990; Schmitz, 1995; Singer et al., 2002; Mackie et al., 2005; Shrimpton et al., 2005; Bystriansky et al., 2006; Hiroi and McCormick, 2007; Nilsen et al., 2007). Within species, quantitative trait locus (i.e. QTL) studies indicate that variation has a genetic basis (Le Bras et al., 2011; Norman et al., 2011; Norman et al., 2012). Genes differentially regulated in Arctic charr (Salvelinus alpinus) with contrasting salinity tolerance capacities have been identified using large-scale gene expression analysis on gill tissue. Evidence suggests that the genes affecting variation in salinity tolerance capacity might be different from the genes that are required for ionoregulation in hyperosmotic conditions (Norman et al., 2014). Transcriptome-wide experiments assessing DE in Arctic charr between freshwater and seawater environments would permit further investigation of these differences. As the number of large-scale expression studies on salinity tolerance increases, comparisons between species may provide insight into the nature of the diverse salinity tolerance capacities that characterize salmonids. Such work would also facilitate further research on the propensity of new candidate genes to confer genetically derived phenotypic variation.

Differences in the genomic organization of salinity tolerance candidate genes may also influence inter- and intraspecific variation in salinity tolerance capacity. Clustered genes in vertebrates are more likely to be co-expressed (Hurst et al., 2004; Ng et al., 2009; Woo et al., 2010) and co-methylated (McGowan et al., 2011) and thus are more likely to be involved in the same or closely related biological pathways. For instance, candidate genes in salmonids show evidence of clustering along linkage groups (Norman et al., 2012) and similar patterns are evident on model teleost chromosomes (e.g. medaka, stickleback, zebrafish) (Norman et al., 2012) and on reconstructed proto-actinopterygian ancestral chromosomes (Kasahara et al., 2007). Current data indicate that interchromosomal rearrangements in the Atlantic salmon genome have produced unique combinations of salinity tolerance candidate genes, which may be related to the high salinity tolerance capacity that is characteristic of this species (Norman et al., 2012). However, these observations are based on relatively few candidate genes. The production of candidate gene lists produced by large-scale gene expression analysis coupled with in silico comparative genomics approaches will facilitate a more in-depth investigation of these patterns, providing additional insight into the genetic architecture of salinity tolerance in salmonids.

In this study we collected large-scale gene expression data from Arctic charr to gain insight into the genes and biological processes that affect salinity tolerance in this species. Specifically, we used high-throughput mRNA sequencing technology to compare expression profiles between fish exposed to freshwater and seawater. We then investigated the genome organization of differentially expressed genes (DEGs) to gain insight into putative non-random patterns in the genetic architecture of salinity tolerance. Our specific objectives were to identify genes differentially expressed in Arctic charr gill tissue after transfer to seawater from freshwater, and to characterize which biological processes affect hypo-osmoregulation in this species. Based upon comparisons with other teleost species, inferences on the locations of the DEGs were made (i.e. putative retention of synteny blocks) and used to test if DEGs are clustered, and whether or not genes or gene clusters are randomly distributed among Arctic charr linkage groups.

Analyses conducted using EBSeq (Leng et al., 2013) revealed that 1075 and 1544 contigs were DE between seawater and freshwater treatments for family 10 (F10) and family 12 (F12), respectively (Fig. 1). For F10, 568 (52.8%) contigs were up-regulated and 507 were down-regulated, while for F12, 812 (53.9%) contigs were up-regulated and 732 were down-regulated. A comparison between families revealed that 257 contigs were DE in both families (supplementary material Table S1). Interestingly, a striking congruence in both the direction and magnitude of expression was observed for all 257 contigs (Fig. 2; F=3118.10; regression d.f.=1; residual d.f.=255; R²=0.9244; P<0.000001). The majority (84.4%) of these exhibited an inter-familial difference in expression ratio that was less than 0.5 (i.e. a 1.4-fold increase) (Fig. 3), and slightly more contigs were up-regulated (51.7%; Table 1) than down-regulated (Table 2).

Alignments against the non-redundant protein database from NCBI were successful for 60,565 of 108,645 assembled contigs (55.7%). To test if they were accurate, we aligned all unannotated contigs (i.e. 40,080) against expressed sequence tags (EST) databases for rainbow trout (Oncorhynchus mykiss) (release 8) and Atlantic salmon (release 6) from the Gene Index Project at the Dana-Farber Cancer Institute (Quackenbush et al., 2001), using criteria described in the Materials and methods. Significant alignments were observed for 36,879 contigs (76.7%). Thus, in total, 97,444 of 108,645 assembled contigs (89.6%) were verified by significant alignments to external databases. The number of annotated DE contigs for F10 and F12 was 832 (supplementary material Table S2) and 1177 (supplementary material Table S3), respectively.

F10 and F12 shared 25 gene ontology (GO) categories significantly enriched for up-regulated genes (Table 3; supplementary material Tables S4–S9). Collectively, these categories represented biological processes related to cell division, cell proliferation and responses to stress. In contrast, F10 and F12 shared only six categories enriched for down-regulated genes, which were related to developmental processes, such as organ development. Tests for disproportionate representation of genes in biological process categories revealed that seven categories for up-regulated genes were different between families (Table 4). Together they represented two general processes, mRNA processing and splicing, and cellular metabolic processes involving nucleic acids. No differences were observed for down-regulated gene categories.

Comparisons of mean intergenic distance (IGD) between stickleback genes and genes DE in seawater gill tissue revealed that DEGs homologous with certain stickleback chromosomes were significantly closer to one another than the respective chromosomal-wide means (Fig. 4; supplementary material Table S10). In F10 these patterns were found on Ga-II, Ga-V, Ga-XV and Ga-XVII, all of which exhibited P<0.001. DEGs that aligned to Ga-II and Ga-V were over 1 Mbp closer to one another on average, while those aligned to Ga-XV and Ga-XVII were 660 kbp and 704 kbp closer to one another, respectively. For eight stickleback chromosomes the chromosomal-wide mean IGD was significantly lower than the mean IGD of the corresponding set of homologous DEGs from Arctic charr F10, which were 460 kbp to 1.13 Mbp further apart. In F12, DEGs that aligned to Ga-I, Ga-IV, Ga-V and Ga-X were significantly closer to one another than the corresponding chromosomal-wide means by at least 480 kbp. With the exception of Ga-X (P=0.0281), these differences were also highly significant (P≤0.0013). Only four stickleback chromosomes (Ga-III, Ga-VI, Ga-VIII and Ga-XIV) had a mean IGD that was significantly lower than the corresponding set of homologous DEGs from F12. Consensus between families was evident for Ga-III, Ga-V and Ga-VIII; however, only for Ga-V were the mean IGDs from DEGs significantly lower than the chromosomal-wide mean.

Translated blast alignments (i.e. BLASTX) of DE contigs against the stickleback cDNA database yielded 507 hits for F10, of which 54.8% (i.e. 278 contigs) were up-regulated. For F12, 695 significant alignments were observed and of these 51.7% (i.e. 359 contigs) were up-regulated. In terms of DEGs located within 40 kbp of one another, 107 DEGs formed 50 clusters of two or more in F10 (supplementary material Table S11), which contained five three-gene clusters and one four-gene cluster. In F12, 252 DEGs formed 105 clusters of two or more genes (supplementary material Table S12). Most genes occurred in clusters of two (72.3%); however, clusters of three (19.0%), four (5.7%), five (1.9%) and six (1.0%) were also evident. In F10, 36 clusters (50% were up-regulated) showed congruency in the direction of expression (i.e. genes in these clusters were either all up-regulated or all down-regulated), while in F12, 48 clusters were congruent (56% were up-regulated). Notably, the proportion of IGDs that were ≤40 kbp was highly similar between all protein-coding genes from the stickleback genome (63918/9180573=0.0069) and the DEGs from F10 (53/8457=0.0063) and F12 (126/18095=0.0069), which suggests that the proportion of DEGs within 40 kbp of one another was the same as would be expected by chance.

Three synteny blocks (i.e. chromosomal segments where a contingent of genes is conserved between two species) contained a greater number of DEGs than expected by chance (Table 5). All were observed in F10, while only the AC-1/Ga-XVII synteny block passed the false discovery rate (FDR) corrected threshold. No synteny blocks deviated from expectations for F12.

This is the first study to use high-throughput RNA sequencing in a controlled laboratory experiment to compare gene expression profiles in the gills of a salmonid species reared in freshwater and seawater. Using a de novo transcriptome assembly we conducted a large-scale analysis to gain insight into the genes and biological processes important for hypo-osmoregulation in Arctic charr. To make inferences about the conservation of candidate genes and biological processes in different taxa, we studied two families and compared our findings with those from similar studies on other salmonids. We examined gene cluster patterns and the distribution of DEGs to test hypotheses about the genetic architecture of salinity tolerance. Our large-scale gene expression analysis has provided a broad perspective of hypo-osmoregulation in Arctic charr, which will be a foundation for future studies to pursue answers to new questions about the physiological and genetic processes that influence hypo-osmoregulation.

Multiple salinity tolerance candidate genes from the seawater MRC complex (Evans et al., 2005; Marshall and Grosell, 2006) were DE in the gill of Arctic charr after transfer to seawater. Among the up-regulated genes were Na⁺/K⁺-ATPase subunit alpha 1b (ATP1α1b), Na⁺/K⁺/2Cl⁻ co-transporter 1a (NKCC1a), cystic fibrosis transmembrane conductance regulator II (CFTRII) and CLDN10E (supplementary material Table S3). These findings confirm previous research in other salmonid species (Singer et al., 2002; Tipsmark et al., 2002; Richards et al., 2003; Mackie et al., 2005; Bystriansky et al., 2006; Nilsen et al., 2007; Larsen et al., 2008; Tipsmark et al., 2008; Madsen et al., 2009; McCormick et al., 2009). Furthermore, ATP1α1a was down-regulated in seawater, as expected given that ATP1α1a and ATP1α1b are reciprocally expressed (Richards et al., 2003; Mackie et al., 2005; Bystriansky et al., 2006; Nilsen et al., 2007; Larsen et al., 2008; Madsen et al., 2009; McCormick et al., 2009) and regulated (McCormick et al., 2009) between freshwater and seawater environments. This consistency with the general literature on hypo-osmoregulation indicates that our experimental design and analysis produced reliable data.

The two families exhibited different gene expression patterns for major transporters from the seawater MRC complex. ATP1α1b, NKCC1a, CFTRI/II and CLDN10E were DE in F12 but not in F10, while ATP1α1a exhibited DE in both families. This difference is consistent with previous research showing that only F12 QTL co-localize with the locations of seawater MRC complex candidate genes (Norman et al., 2014). Genome scans for salinity tolerance QTL using the same fish revealed that only 17.6% (i.e. 3/17) were detected in both families. Moreover, QTL co-localize with ATP1α1b and NKCC1 on Arctic charr linkage groups in F12 but not F10 (Norman et al., 2012). Thus differences in salinity tolerance between F10 and F12 are evident at both nucleotide and transcriptional levels.

Our findings suggest that hypo-osmoregulation in Arctic charr may be regulated by a ‘core’ gene subset that has limited capacity for plasticity in expression. In total, 257 DEGs from both families had expression ratios showing a remarkable consistency in magnitude and direction (Fig. 1). Predicting whether this pattern is conserved in other populations or in other salmonids is difficult, as current evidence is equivocal. For instance, the consistency in expression of this core subset might be related to the shared genetic background of the families, as both were derived from the same strain. Nevertheless, both families were generated by outbred crosses from unrelated parents and results from QTL experiments do support the concept of genetic heterogeneity between the families (Norman et al., 2011).

Comparing expression profiles among genes from the core subset provides clues about the molecular pathways involved, permitting inferences about biological significance to be made. For instance, changes in the direction of transcription of genes involved in methylation suggest that methylation states were altered in seawater-exposed fish. Adenosylhomocysteinase was up-regulated, and as part of the methionine cycle (Fig. 5), this enzyme catalyses a reversible reaction that is thermodynamically favoured towards the production of S-adenosylhomocysteine from homocysteine and adenosine. S-adenosylhomocysteine is a potent inhibitor of methyltransferases (Hoffman et al., 1980), leading to increased hypomethylation levels of proteins, lipids, polyamines and chromatin, which in turn have been linked to altered gene expression and cellular differentiation states (reviewed in James et al., 2002). Adenosine removal by adenosine kinase (up-regulated in F10) may reduce the availability of adenosine to adenosylhomocysteinase, thereby promoting conditions favourable to the formation of homocysteine and adenosine. Gene transcription profiles for other components of the methionine cycle provide further support for this proposition: methyltransferases were up-regulated (both families); 5′-nucleotidase, which converts adenosine monophosphate to adenosine, was downregulated (F12); and methionine adenosyltransferases, which catalyse the formation of S-adenosylmethionine from methionine, were up-regulated (F10). Because the gene sets that affect hyper- and hypo-osmoregulation are different, changes in methylation states should be anticipated when transcription is compared between freshwater- and seawater-exposed fish. Observations from other species support these findings: seawater-induced changes in methylation states have been observed in brown trout (Morán et al., 2013), and transcription of the betain-homocysteine S-methyltransferase gene in ayu (Plecoglossus altivelis) gill is reduced after transfer from freshwater to brackish water (Lu et al., 2010). Thus, changes in methylation states are likely to be part of the physiological transition to seawater in anadromous fishes.

Comparisons between transcriptional profiles for genes that are part of the same molecular pathways as inducible nitric oxide synthase (iNOS) provide clues about why immune response gene regulation coincides with elevated salinity and salinity tolerance capacity (Norman et al., 2014). Inducible nitric oxide synthase, highly up-regulated in the core subset, produces nitric oxide (NO), a bio-signalling molecule that elicits many physiological responses (Pacher et al., 2007). NO produced by iNOS in particular, triggers an immune system response to pathogenic bacteria and viruses (Lowenstein and Snyder, 1992). NO also increases the expression of cytochrome P450 (Zamora et al., 2001), which in turn produces hydrogen peroxide that is catalysed with chloride by myeloperoxidase (MPO) to produce hypochlorous acid (HOCl), a potent antimicrobial agent (Klebanoff, 2005). Interestingly, chloride availability may limit MPO activity (Pacher et al., 2007). If true, elevated plasma chloride levels may provide conditions that are conducive towards strengthening antimicrobial defence, which could have implications for ionoregulation, as hypochlorous acid production would necessitate chloride removal from blood plasma.

Elevated iNOS transcription probably fostered increased peroxynitrite formation, suggesting that further research about the putative impact and role of peroxynitrite in hypo-osmoregulation may be warranted. Peroxynitrite is a highly potent oxidizing compound formed from the spontaneous reaction of NO with superoxide (Huie and Padmaja, 1993). Although superoxide can be cytotoxic, it is detoxified by scavenging enzymes known as superoxide dismutases (SOD) (Pacher et al., 2007). In F12, SOD was down-regulated in seawater, suggesting that superoxide detoxification was reduced during hypo-osmoregulation. Moreover, superoxide detoxification by SOD is three times slower than peroxynitrite formation (Pohanka, 2013). Therefore, increased superoxide levels combined with elevated NO production would probably result in increased peroxynitrite levels in seawater gill. Reduced transcription of SOD in high salinity tolerance capacity Arctic charr (Norman et al., 2014) also suggests that variation in salinity tolerance capacity may correlate with variation in peroxynitrite levels. Although peroxynitrite can protect against invading pathogens, it is also cytotoxic and can induce apoptosis and necrotic cell death (Pacher et al., 2007). Protection from pathogens would probably ease the transition to a hypo-osmoregulatory state, which can be metabolically demanding for salmonids (Rao, 1968; McCormick et al., 1989; Morgan and Iwama, 1999). Additionally, increased rates of cell death may accelerate cellular turnover in gill epithelia, where freshwater-type MRCs are replaced with seawater-type MRCs (McCormick et al., 2009). Hypomethylated DNA is also more susceptible to oxidative stress, leading to increased rates of apoptosis and cellular turnover (James et al., 2002). Elevated cellular turnover may, in turn, more rapidly increase hypo-osmoregulatory capacity, which would conserve energy and mitigate physiological stress conferred by a suboptimal osmoregulatory condition. Notably, peroxynitrite also modulates Na⁺/K⁺-ATPase activity via cysteine oxidation (Muriel and Sandoval, 2000; Varela et al., 2004).

While certain genes from the core subset have been implicated in hypo-osmoregulation by other studies, the majority are DE only in Arctic charr, providing limited support for the hypothesis that it is functionally conserved in other salmonids. We compared the results of this study with similar research on Atlantic salmon (Seear et al., 2010; Robertson and McCormick, 2012) and sockeye salmon (Evans et al., 2011). Only nine genes in sockeye salmon (Evans et al., 2011) and seven genes in Atlantic salmon (Seear et al., 2010; Robertson and McCormick, 2012) showed patterns of DE that were consistent with our observations (Table 6). It must be emphasized that differences in experimental design and in the technical platforms used by others (Seear et al., 2010; Evans et al., 2011; Robertson and McCormick, 2012) may account for these differences.

The over-representation of up-regulated cell division genes suggests that MRC proliferation was occurring in the gill of seawater-acclimating fish. Though short-term (e.g. 5 h) acclimation to a hyperosmotic environment is characterized by the transformation of freshwater- to seawater-type MRCs (Hiroi et al., 1999), long-term (i.e. days to weeks) acclimation involves an increase in the number of MRCs (McCormick, 2001; Evans et al., 2005; Hwang and Lee, 2007). In the Atlantic salmon gill, the seawater isoform ATP1α1b occurs in separate MRCs than its freshwater counterpart, ATP1α1a (McCormick et al., 2009). Considering that seawater acclimation in salmonids entails elevated ATP1α1b expression (Richards et al., 2003; Mackie et al., 2005; Bystriansky et al., 2006; Nilsen et al., 2007; Larsen et al., 2008; Madsen et al., 2009; McCormick et al., 2009), and that tissue samples were collected 10 days after transfer to seawater, it is plausible that the enrichment of up-regulated cell division genes could reflect the proliferation of ATP1α1b-type seawater MRCs. However, gill epithelium is composed of many different cell types, and pavement cells, which have a limited role in ionoregulation, constitute most of the epithelial surface area (>90%), followed by MRCs (<10%) (Wilson and Laurent, 2002; Evans et al., 2005). Therefore, cell division gene enrichment is also probably evidential of the proliferation of other cell types (see Wilson and Laurent, 2002). For instance, iNOS, which may play a role in hypo-osmoregulation, localizes to small cells buried deep within gill filaments (Ebbesson et al., 2005). Enhanced expression of the methionine cycle, specifically methionine adenosyltransferase, is also regarded as essential in regulating cell cycle turnover rates (Lin et al., 2014), which complements the results from this study (Fig. 5). Enrichment of up-regulated genes related to cell division GO categories has also been described for other teleosts: cell-cycle regulation in sea bass (Boutet et al., 2006), cellular processes in tilapia (Tine et al., 2008), and cell-cycle in goby (Evans and Somero, 2008).

GO analyses suggest that the repression of developmental gene expression may be related to increased osmoregulatory stress. Over-representation of down-regulated genes was confined to a GO branch related to development (i.e. multicellular organismal development, system development, developmental processes, anatomical structure development, multicellular organismal process and organ development). This appears paradoxical because cell division is integral to organ development and thus these processes would be expected to exhibit some congruency in expression. However, cellular turnover driven by senescence and apoptosis can occur independently of development (Pellettieri and Sánchez Alvarado, 2007), while development necessitates cell division. Given that stress response genes were also up-regulated and over-represented, we hypothesize that the stress induced by hyperosmotic conditions may cause a repression of developmental gene expression. This would have to be temporary, however, perhaps until gill tissue and other ionoregulatory organs have transitioned to an optimal ion-excreting state, because the seawater phase of salmonid ontogeny is typically characterized by heightened growth rates. Repression of developmental gene expression may liberate energy for the transition of osmoregulatory tissues towards an ion-excreting state, as the metabolic demand of hypo-osmoregulation can be high, ranging from <4% in cut-throat trout (Oncorhynchus clarkii clarkii) (Morgan and Iwama, 1999) to 30% in rainbow trout (Rao, 1968), and 20% of gill metabolism in Atlantic salmon (McCormick et al., 1989). To our knowledge, there are no data on the metabolic demands of hypo-osmoregulation in Arctic charr, although our hypothesis suggests that the cost would be high. The current study may have captured the transition period in gill physiology when hyper-osmoregulatory cells are being ablated prior to the metabolic machinery of hypo-osmoregulatory cells being fully activated.

Comparative genomics analyses provide support for the hypothesis that the relative location of salinity tolerance candidate genes along a chromosome is a feature of the genetic architecture of hypo-osmoregulation. The physical clustering within linkage groups of DEGs in seawater gill tissue is consistent with previous observations in teleosts (Norman et al., 2012). The majority of genes that met our clustering criteria (i.e. ≤40 kbp) were also identified as paired by previous research conducted on the same fish, but on a smaller scale (Norman et al., 2012). Furthermore, DEGs that align to certain stickleback chromosomes tend to have lower mean IGDs than the corresponding chromosomal-wide means, and several of these chromosomes are homologous with salinity tolerance QTL. For instance, Ga-I, which contains DEGs from F12 that are on average 890 kbp closer together, is homologous with AC-4 and AC-12 (Norman et al., 2012), where salinity tolerance QTL have been localized in the same family (Norman et al., 2011). In addition, Ga-I not only contains two key genes from the seawater MRC model (i.e. ATP1α1b and CLDN10E), but an inversion on this chromosome is correlated with disparity between marine and freshwater stickleback ecotypes (Jones et al., 2012). Interestingly, ATP1α1b is positioned within this inversion while CLDN10E is positioned outside it. Consequently, ATP1α1b and CLDN10E are 1.007 Mbp apart on Ga-I in the freshwater ecotype and 1.012 Mbp apart in the marine ecotype, which suggests that ecotype development could be influenced by changes in the relative proximity of these two genes along the chromosome. Salinity tolerance QTL localized to Atlantic salmon linkage groups AS-5 and AS-22 are also homologous with Ga-I (Norman et al., 2012). Norman et al. (Norman et al., 2012) hypothesized that the superior hypo-osmoregulatory capacity that is characteristic of Atlantic salmon (Hoar, 1988; Hiroi and McCormick, 2007) could be related to the formation of novel gene clusters by interchromosomal rearrangements. Thus evidence suggests that a contingent of genes on Ga-I and their homologous counterparts in other species (e.g. Arctic charr and Atlantic salmon) may have implications for hypo-osmoregulation. Other stickleback chromosomes where DEGs are significantly closer than expected and that share homologous affinities with Arctic charr salinity tolerance QTL are AC-18 and Ga-V (F12), AC-4 and Ga-I (F10), and AC-31 and Ga-X (F10) (Norman et al., 2011).

The genetic architecture of hypo-osmoregulation in Arctic charr shows some evidence of large-scale structure. We found that stickleback synteny blocks on Arctic charr linkage groups AC-1, AC-11 [henceforth referred to as AC-1/11, due to the strong homeologous affinities of these linkage groups (Danzmann et al., 2005; Timusk et al., 2011; Norman et al., 2012)], and AC-28 contain a greater number of DEGs than expected by chance. These linkage groups also contain QTL for the salinity tolerance performance traits Na⁺/K⁺-ATPase activity and growth in seawater (Norman et al., 2011), suggesting that these QTL could be the product of multiple genes containing cis-regulatory polymorphisms. This also suggests that some of the phenotypic variation attributed to QTL could stem from disparity in expression rate and transcript availability. In addition, stickleback chromosomes Ga-XVII and Ga-XIV share syntenic blocks with AC-1/11 and AC-28, respectively. Together with their respective homeologues, Ga-XVII&XII and Ga-XIII&XIV primarily originated from the ancestral L and I chromosome lineages of actinopterygian fishes, respectively (Kasahara et al., 2007; Nakatani et al., 2007; Norman et al., 2012), and were duplicated with the whole-genome duplication event in the teleost ancestor ~350 MYA (Christoffels et al., 2004). Comparisons with the teleost ancestor reveal that Ga-XVII&XII and Ga-XIII&XIV have been affected by very few interchromosomal rearrangements since the stickleback and teleost ancestral lineages diverged. Furthermore, the interchromosomal rearrangements that have occurred reflect translocations mostly between these homeologous pairs. Ga-XVII (orthologous with ancestral chromosome L), contains a 575 kbp segment that is homologous with Ga-XIII (orthologous with ancestral chromosome I), which localizes to the AC-1/11 and Ga-XVII synteny block (Norman et al., 2012). Close inspection of this region reveals that a single gene, tubulin polyglutamylase complex subunit 2 (TPGS2), was DE in seawater. These findings suggest that interchromosomal rearrangements between orthologues of L and I in the most recent common ancestor of stickleback and Arctic charr could have contributed to the apparent large-scale structure of AC-1/11. Furthermore, this suggests that interchromosomal rearrangements may have played a role in the evolution of hypo-osmoregulation in Arctic charr, and is consistent with the findings of Norman et al. (Norman et al., 2012), who found that interchromosomal rearrangements produced a coalescence of salinity tolerance candidate genes on certain chromosomes in Atlantic salmon, which they hypothesized could contribute to the high salinity tolerance capacity of that species. Further research is necessary to test these hypotheses.

This study characterized changes in gene expression profiles and biological processes that occurred in Arctic charr gill tissue during seawater acclimation. Gene expression profiles of major transporters from the seawater MRC complex were consistent with previous research for one family, and inconsistent for the other family. Despite these differences, a subset of DEGs in both families showed a remarkable consistency in both magnitude and direction of expression. We hypothesized that this subset may comprise genes that have a limited capacity for plasticity in expression, which may be integral for hypo-osmoregulation in Arctic charr. Examination of these genes suggested that in seawater-exposed fish, methylation states were changing and that peroxynitrite formation and antimicrobial defence may have implications for hypo-osmoregulation. Based on patterns from GO, we postulated that the enrichment of cell division genes could reflect the proliferation of ATP1α1b-type seawater MRCs, and that the stress induced by hyperosmotic conditions may have caused a repression of developmental gene expression. Our findings provide additional support for previous research that suggests the relative proximity among a contingent of genes on Ga-I has important implications for hypo-osmoregulation. We also found that clustering patterns of DEGs were consistent with previous observations in salmonids and other teleosts, and that AC-1/11 and AC-28 show evidence of large-scale structure for hypo-osmoregulation candidate genes. Finally, we provided evidence that interchromosomal rearrangements have played a role in the evolution of hypo-osmoregulation in Arctic charr.

Two full-sib Arctic charr (Salvelinus alpinus Linnaeus 1758) families (F10 and F12) were created from fish descended from the anadromous population resident in the Fraser River in Labrador, Canada. Crosses were made in November of 2006 at the Coastal Zones Research Institute, Shippagan, New Brunswick (Canada). Eight months after hatching, 150 progeny per family were implanted with passive integrated transponder tags and transferred to St Andrews Biological Station, St Andrews, New Brunswick (Canada). There they were held in 1 m³ freshwater tanks and reared under controlled simulated-natural conditions for photoperiod and water temperature. Tanks contained filtered and aerated freshwater (9.9–10.7°C, flow rate 18 l min⁻¹, dissolved O₂ 10.0–10.6 mg l⁻¹). For the duration of sampling (i.e. 9 June 2008 to 6 July 2008) a 16 h:8 h light:dark photoperiod regime was used to minimize the confounding effects associated with a naturally changing photoperiod (for details, see Norman et al., 2011). All experiments were carried out in accordance with the Canadian Council for Animal Care Guidelines under the University of Guelph Animal Care Committee approved protocol number 08R033.

The salinity tolerance trials were conducted during June 2008. Freshwater was replaced with 100% filtered seawater (salinity 31,000–33,000 mg l⁻¹, 10.5–11.9°C, flow rate 18 l min⁻¹, dissolved O₂ 8.1–11.4 mg l⁻¹) over 24 h. After 10 days exposure to full strength seawater non-lethal gill tissue biopsies (McCormick, 1993) were collected for 50 fish per family. Samples were also collected from freshwater fish at this time. Each sample consisted of ~3–6 primary filaments, which were stored in 1 ml of RNAlater (Life Technologies) at 4°C for 24 h, then transferred to −20°C for long-term storage. To compare expression profiles in gill tissue between fish exposed to seawater and freshwater, we selected six fish from seawater and three fish from freshwater from each family for mRNA sequencing (N=18). Notably, the same seawater-exposed fish used for this study were also used for a companion study, and thus were individuals that exhibited the greatest and least number of salinity tolerance QTL (see Norman et al., 2014).

Immediately prior to homogenization, gill tissue samples were removed from RNAlater and rinsed. Trizol reagent (Life Technologies) was used to extract total RNA according to the manufacturer's instructions. A disposable pestle grinder system (Thermo Fisher Scientific, Waltham, MA, USA) was used to homogenize tissue immersed in 1 ml Trizol. Purity of extracted total RNA was assessed with a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific). For further details regarding sample preparation, we refer the reader to Norman et al. (Norman et al., 2014). The Bio-Rad Experion system (Bio-Rad, Hercules, CA, USA) was used to assess RNA integrity. Only samples with an RNA quality index (RQI) of 8.0 or greater were submitted for mRNA sequencing. mRNA libraries were constructed separately for all 18 fish. Illumina sequencing was performed by the Clinical Genomics Centre (Mount Sinai Hospital, Toronto, Ontario, Canada). Sequencing libraries were created according to the manufacturer's instructions (Illumina). Briefly, poly-A mRNA was isolated from total RNA and purified using poly-T oligo-attached magnetic beads. mRNA was fragmented into pieces using divalent cations under elevated temperature, and copied into cDNA using reverse transcriptase and random primers. Following the ligation of adaptors, cDNA templates were purified and amplified with PCR. Sequencing was performed with a HiSeq 2000 instrument (Illumina, San Diego, CA, USA) and subjected to 100 cycles of paired-end sequencing. Image analysis, base-calling and quality value calculation were performed using Illumina's sequence analysis software, Casava (version 1.8.2). The raw reads are available from the Sequence Read Archive (SRA) at the National Centre for Biotechnology Information (NCBI), under BioProject accession no. SRP026259 and BioSample accession no. SRS452215.

The procedures used for read pre-processing, the measures taken for quality control, and the steps for de novo construction of the reference transcriptome are described in detail by Norman et al. (Norman et al., 2014). Briefly, the reads used to quantify expression were filtered and trimmed according to quality scores using FastQC and Trimmomatic software (Andrews, 2010; Lohse et al., 2012). To facilitate de novo assembly, reads with rare k-mers, low quality reads and identical reads were removed using Rnnotator (Martin et al., 2010). The Velvet-Oases software pipeline was used to generate a single merged assembly from multiple assemblies spanning a range of k-mers. Redundancy in the merged assembly was reduced using CD-HIT-EST (Li and Godzik, 2006), producing a final consensus assembly of 108,645 contigs.

Read alignment to the reference transcriptome was conducted with RNA-Seq by expectation maximization (RSEM) (Li and Dewey, 2011). RSEM performs accurate transcript quantification when a reference genome is unavailable. It also allows for the inclusion of non-uniquely mapping reads in abundance estimates through the use of a statistical model that accounts for read mapping uncertainty. Because scaffolding was not used during assembly with Velvet-Oases, paired-end reads could not be used for accurate expression quantification. Accordingly, only forward reads were mapped as they typically yielded higher quality scores than reverse reads. Differential expression between fish exposed to seawater and freshwater was analysed separately for each family, using EBSeq (Leng et al., 2013). This software uses an empirical Bayes hierarchical model that accommodates read mapping uncertainty from RSEM. To account for differences in library size, the number of reads in each sample was adjusted by median normalization (Anders and Huber, 2010; Leng et al., 2013). Differential expression was established at a threshold of P≤0.05 and a Benjamini–Hochberg (B–H) FDR of α=0.05 (Sokal and Rohlf, 1995). Linear regression was used using SPSS 16 (IBM) to compare the magnitude and direction of expression ratios for genes that were DE in both families.

Sequences were annotated using Blast2GO (Götz et al., 2008). The single best alignment for each contig was obtained using BLASTX against the non-redundant protein sequence database at NCBI at an E-value threshold of 1×10⁻⁵.

The Biological Networks GO tool (BiNGO) (Maere et al., 2005) was used to perform hypergeometric and B–H FDR tests (α=0.05) to identify over-represented GO categories under the Biological Processes umbrella. Independent tests were conducted for each family for up-regulated and down-regulated genes. Data were entered in the form of official gene symbols and tests were run against the annotated reference de novo constructed transcriptome. All annotations were obtained from Homo sapiens. Fisher's exact tests followed by B–H FDR tests (α=0.05) were performed to determine whether any GO categories contained different quantities of genes between families. These tests were done for up-regulated and down-regulated gene groups independently.

BLAST+ (version 2.2.26) (Camacho et al., 2009) was used to align all contigs showing DE between seawater and freshwater groups to a database of cDNA sequences for protein-coding genes in stickleback using the same alignment criteria noted previously (Ensembl version 68) (Jones et al., 2012; Flicek et al., 2013). Based on a synteny map of Arctic charr and stickleback (Norman et al., 2012), patterns between gene location and DE were assessed to ascertain the extent that DEGs cluster along Arctic charr chromosomes. Genes were considered clustered if their coding regions were within 40 kbp of one another on a stickleback chromosome (Ng et al., 2009), as Ng et al. (Ng et al., 2009) detected a significant negative correlation (R²=0.50) between gene co-expression and IGD in zebrafish, which was strongest (R²=0.79) for genes 10 to 40 kbp apart.

To test whether DEGs were located in closer proximity to one another along a chromosome, permutation analysis was performed using the R statistical language (R Core Team, 2013). The analysis included all of the DEGs that aligned to the stickleback genome. Tests were conducted independently for each chromosome. The IGDs among all genes on a chromosome were pooled with those from the DEGs that aligned to that chromosome. The difference between means was estimated between a random sample of the pooled data and the remaining non-sampled data over 10,000 iterations. The size of each random sample was equal to the number of DEGs that aligned to the chromosome in question. The observed mean difference was then compared with the distribution of mean differences from the 10,000 random samples. Significance was denoted for P≤0.05 following a B–H FDR test to correct for multiple tests (Sokal and Rohlf, 1995). This analysis assumes that the relative proximity and distribution of genes located on stickleback chromosomes are conserved on Arctic charr chromosomes. We acknowledge that this may not be realistic across entire chromosomes, although conserved synteny blocks between stickleback and Arctic charr chromosomes have been described (Norman et al., 2012).

To determine if the quantity of DEGs in each of the Arctic charr–stickleback synteny blocks was greater than expected, two-tailed χ² tests of association were performed using 2×2 contingency tables. This test is an expansion of the analysis originally performed and described in detail by Norman et al. (Norman et al., 2013). Briefly, for each synteny block and family, tests compared the quantity of genes showing significant DE in a synteny block and the total quantity of genes contained within that synteny block, with the quantity of genes in these two classes that were not contained within the synteny block. These tests assumed that the composition of genes within synteny blocks was not different between Arctic charr and stickleback. An adjustment to correct for multiple tests was done using a B–H FDR test (α=0.05) (Sokal and Rohlf, 1995).

The sample size was not large enough to warrant statistical analyses of co-variation between gene expression and salinity tolerance phenotypes. However, in a previous study with the same fish (Norman et al., 2011), significant associations were found between variation in gill filament Na⁺/K⁺-ATPase activity and blood plasma osmolality with genetic variation on several Arctic charr linkage groups (i.e. chromosomes). Given that the fish chosen for this study were those individuals with the greatest and least number of salinity tolerance QTL, and considering that a QTL represents a statistical association between phenotypic variation and genetic variation, changes in gene expression are related to physiological responses.

The authors acknowledge Claude Pelletier of the Coastal Zones Research Institute for generating the family crosses, and Brian Glebe, Steven Leadbeater and Wilfred Young-Lai for overseeing the rearing and maintenance of fish at St Andrews Biological Station. We also thank Xuejiang (Roger) Shi of the Clinical Genomics Centre at Mount Sinai Hospital for overseeing the mRNA sequencing. We give special thanks to Xia Yue, Andrea Kocmarek and Jing Zhang for help and advice in the laboratory. We also thank two anonymous reviewers for their comments, which improved the quality of the manuscript.