Arctic charr brain transcriptome strongly affected by summer seasonal growth but only subtly by feed deprivation

Abstract

Background

The Arctic charr (Salvelinus alpinus) has a highly seasonal feeding cycle that comprises long periods of voluntary fasting and a short but intense feeding period during summer. Therefore, the charr represents an interesting species for studying appetite-regulating mechanisms in fish.

Results

In this study, we compared the brain transcriptomes of fed and feed deprived charr over a 4 weeks trial during their summer feeding season. Despite prominent differences in body condition between fed and feed deprived charr at the end of the trial, feed deprivation affected the brain transcriptome only slightly. In contrast, the transcriptome differed markedly over time in both fed and feed deprived charr, indicating strong shifts in basic cell metabolic processes possibly due to season, growth, temperature, or combinations thereof. The GO enrichment analysis revealed that many biological processes appeared to change in the same direction in both fed and feed deprived fish. In the feed deprived charr processes linked to oxygen transport and apoptosis were down- and up-regulated, respectively. Known genes encoding for appetite regulators did not respond to feed deprivation. Gene expression of Deiodinase 2 (DIO2), an enzyme implicated in the regulation of seasonal processes in mammals, was lower in response to season and feed deprivation. We further found a higher expression of VGF (non-acronymic) in the feed deprived than in the fed fish. This gene encodes for a neuropeptide associated with the control of energy metabolism in mammals, and has not been studied in relation to regulation of appetite and energy homeostasis in fish.

Conclusions

In the Arctic charr, external and endogenous seasonal factors for example the increase in temperature and their circannual growth cycle, respectively, evoke much stronger responses in the brain than 4 weeks feed deprivation. The absence of a central hunger response in feed deprived charr give support for a strong resilience to the lack of food in this high Arctic species. DIO2 and VGF may play a role in the regulation of energy homeostasis and need to be further studied in seasonal fish.

Background

Feeding is pivotal for animals in order to sustain their energy and substrate needs to live, grow and reproduce. In mammals, energy intake and expenditure are tightly regulated by a crosstalk of peripheral and central signalling actors and pathways [1]. Peripherally derived hunger (orexigenic) and satiety (anorexigenic) signals as well as long-term signals reporting energy status are perceived and processed in a number of brain nuclei in order to control short-term (meal-to-meal) appetite and long-term energy homeostasis [2, 3]. Of these, the arcuate nucleus (ARC) in the hypothalamus represents the pivot for controlling food intake and energy balance [4]. The ARC contains two populations of “first order” neurons, one expressing the anorexigenic proopiomelanocortin (POMC) and cocaine-and amphetamine regulated transcript (CART), the other the orexigenic agouti-related peptide (AgRP) and neuropeptide Y (NPY) [3]. These project to “second order” neurons that transduce orexigenic and anorexigenic signals via NPY and melanocortin receptors (MCR) [1]. While NPY signalling through its receptors causes an orexigenic response, signalling through MCRs results in either an anorexigenic or an orexigenic response. POMC-derived α-melanocyte-stimulating hormone (α-MSH) is a melanocortin 4 receptor (MC4R) agonist and a potent appetite suppressor in mammals [5]. AgRP, on the other hand, is an inverse agonist to the constitutively active MCRs and increases food intake [6]. These appetite and energy signalling neuropeptides have been shown to be evolutionary conserved [7, 8] and to be involved in the control of food intake in fish [9, 10]. However, responses of central appetite regulators to energy perturbation vary across species and even within species depending on the experimental design [11]. For example, NPY expression is higher after 7 days of feed deprivation in the hypothalamus of zebrafish (Danio rerio) [12] and in the preoptic area of chinook salmon (Oncorhynchus tshawytscha) and coho salmon (Oncorhynchus kisutch) [13], whereas in cunner (Tautogolabrus adsperus), hypothalamic NPY expression remains unaffected after 7 days feed deprivation [14]. Such differences in the response to feed deprivation are not unexpected as fish represent the most diverse group of vertebrates with a myriad of adaptations to spatially different and temporal changing environments. For example, in the high-latitude inhabiting anadromous (sea-migrating) Arctic charr (Salvelinus alpinus), food intake varies dramatically from little or no feeding while residing in fresh water during winter to voracious feeding during their short summer residence in the sea [15, 16]. This behaviour, which presumably developed as a response to seasonal and spatial differences in water temperature and food availability at high latitudes, now appears to be regulated independently of these factors; captive offspring of anadromous Arctic charr exhibit pronounced seasonal changes in appetite and growth when held at constant temperature and given food in excess [17]. Despite such seasonal changes in food intake, the expressions of orexigenic and anorexigenic appetite regulators in different brain regions is similar between anorexic winter charr with hyperphagic summer charr [18]. Furthermore, the expressions of orexigenic and anorexigenic neuropeptides in the hypothalamus are unaffected by short- and long-term feed deprivation in this species [19]. However, in these and in most other studies investigating appetite regulation in fish, expression levels of known appetite regulators have been measured by RT-qPCR, a method with the drawback of restricting the focus to a limited number of genes. The fact that novel actors in the complex control of food intake are still being discovered calls for a more global approach when investigating appetite regulation in fish. Today, high-throughput RNA sequencing is a powerful tool in experimental biology. Transcriptomic approaches have been applied in various contexts to improve knowledge of the biology of the seasonal Arctic charr [20, 21], albeit not with focus on appetite regulation. Consequently, we sequenced the brain transcriptome of fed and feed deprived charr during their natural summertime growth in an attempt to advance the knowledge on global responses to feed deprivation, assess alterations of known central appetite regulators, and to identify possible new actors involved in the control of appetite and energy metabolism in fish.

Results

Morphometrics

The experiment was performed on two-year old immature anadromous Arctic charr held under natural photoperiod (24 h light at that time of the year) and ambient water temperature. Growth development during the 4 weeks experimental period was assessed on subsamples of 15 tagged fish per treatment group. Average body mass and condition factor were 98.0 ± 4.2 g and 1.04 ± 0.02 on June 25, 2014, the start of the experiment (T₀) (Fig. 1). On July 23, the end of the experiment (T₁), body mass and condition factor were markedly higher in fed (Fed) than in feed deprived (FDP) charr with 161.4 ± 4.7 g and 1.27 ± 0.01, and 84.4 ± 5.2 g and 0.86 ± 0.02 in the fed and FDP charr, respectively. Average body mass and condition factor of the 5 fish sampled for transcriptomics fell in the same range as for the tagged fish and were 108.3 ± 10.2 g and 1.05 ± 0.03 at T₀ and 157.2 ± 8.3 g and 1.25 ± 0.03 (Fed) and 80.6 ± 6.3 g and 0.92 ± 0.04 (FDP) at T₁.

Fig. 1

Fulton’s condition factor (a) and body mass (b) of tattooed fed (black dots) and feed deprived Arctic charr (white dots) during the experiment. Charr sampled at T₀ were measured before distribution to Tank 1 and Tank 2. Dashed line: water temperature. n = 21 per treatment group. Values are shown as mean ± SEM. Different capital and lowercase letters denote differences within treatment group at different time points and differences between treatment groups at the given time point, respectively

De novo transcriptome assembly and annotation

The de novo brain transcriptome assembly produced 49829 contigs with a FPKM greater than 1 for at least one of the 15 samples. Their total length equalled 84028148 base pairs. The N50 (i.e., the contig length that produces half the bases of the assembly) reached 2663 bp (Table 1).

Table 1 General statistics of contigs generated by RNAseq technology for brain gene expression characterization. Only contigs possessing a FPKM greater than 1 for at least one library were considered for annotation and further expression analysis

The annotation rate reached 85.07% of the 49829 contigs. The Atlantic salmon (Salmo salar) contributed most to the annotation of the Arctic charr brain transcriptome with 57.6% of the contigs annotated on this species (Fig. 2). All other species contributing to the annotation, except Homo sapiens, were exclusively teleosts including another salmonid, the rainbow trout (Oncorhynchus mykiss). One or more GO identifier could be assigned to 13231 out of the 49829 contigs. Different approaches were used to verify the quality of the contigs. First, the assembly quality and annotation completeness of the transcriptome were assessed by BUSCO analysis. Out of 4584 single-copy ortholog genes common to Actinopterygii, the assembly was 67.2% complete (2560 complete single-copy BUSCOs and 524 complete duplicated BUSCOs), while only 2.8% of contigs were fragmented (130 BUSCOs) and 30.0% were missing (1370 BUSCOs). Second, a comparison to a phylogenetically close species reference proteome was achieved using the Atlantic salmon database from the NCBI. Out of 97555 Salmo salar proteins, 47419 were aligned with at least 80% identity over 80% of their length on our assembly, which corresponded to 12238 Arctic charr contigs. Third, initial reads were mapped to the contigs in order to validate to what extent the contigs reflected the initial information. The contigs had high realignment rates that ranged between 88 and 89% depending on the sample (15 libraries). Furthermore, the construction of the RNAseq data set was verified by mapping the reads from an Arctic charr gill transcriptome [22] on our set of contigs. On average, 79.35% of reads from this gill transcriptome (SRA accession: SRX314607) were aligned to the 15 Arctic charr libraries. Finally, both proteins and RNAs from the S. alpinus reference genome (NCBI ID: 12179) [23] were aligned to our set of contigs. Out of 59926 proteins, 49212 (82%) had a hit on the contigs, 31316 of these proteins shared over 80% similarity and 80% of coverage with the built contigs. Regarding RNAs from the reference genome, 76% (51217 out of the 67196) had a hit on our contigs, with 24432 RNAs having over 80% of similarity and 80% of coverage with the contigs.

Fig. 2

Top 20 species contributing the most to best-hit annotations. Best hits were based on all databases involved in the annotation process. Yellow: teleosts, orange: salmonids, blue: mammals

Gene expression patterns and differential gene expression

In order to inspect the overall expression patterns of the 15 samples, a correlation heatmap based on the raw counts of the 49829 contigs possessing a FPKM greater than 1 for at least one sample was drawn (Fig. 3). This analysis highlighted that all samples from T₀ clustered together and were markedly different from samples taken at T₁. Samples of the fed and feed deprived charr at T₁ did not cluster in accordance with the treatment group (Fig. 3).

Fig. 3

Correlation heatmap based on raw counts of the 49829 contigs possessing a FPKM greater than 1 for at least one sample. X01-X15: number of sampled fish during the experiment. T0: fish sampled at the start of the experiment (June 25). FED: fed fish sampled at the end of the experiment (July 23). FDP: feed deprived fish sampled at the end of experiment (July 23)

In fed Arctic charr, 2819 contigs differed over time (Table 2). This corresponded to 5.7% of all sequenced contigs. Among these, 1534 contigs were up-regulated, while 1285 were down-regulated.

Table 2 Number of up- and down-regulated transcripts in the different comparisons returned by EdgeR analysis. Cut-off at FDR < 0.05 and at LogFC 0.5/−0.5

In feed deprived charr, 4570 (9.2%) contigs differed over time (Table 2). Of these 2616 were up-regulated while 1954 were down-regulated from start to the end of experiment. At T₁, only 175 contigs (0.4% of all contigs) were found to be differentially expressed between FDP and Fed charr (Table 2), with 68 contigs being up-regulated and 107 being down-regulated. Matching the six lists of up- and down regulated contigs from the three comparisons in Up-set graphs enabled us to break down the lists of differentially expressed genes and find the intersection points of all three comparisons (Fig. 4). Over time, 1449 contigs were up-regulated and 825 down-regulated in both FDP and fed charr. Interestingly, at the same time, 1421 contigs were uniquely up-regulated in the FDP versus T₀ comparison.

Fig. 4

Number of up- and down-regulated transcripts in the different pairwise comparisons. The total number of differentially expressed transcripts for the different comparisons are presented as horizontal bars. Vertical bars represent lists of differentially expressed genes that were only found in one specific comparison or shared by two or more comparisons, indicated by dots and dots connected by lines, respectively, below

In contrast, only 825 contigs were uniquely up-regulated in fed charr over time. Similarly, 1058 contigs, were uniquely down-regulated in FDP charr over time, whereas only 337 contigs were uniquely down-regulated contigs in the Fed versus T₀ comparison. Chi-square tests of the 2 × 2 contingency tables of up- and downregulated genes indicated that there was a correlation of the proportion between uniquely and shared differentially expressed contigs and feeding regime over time (up-regulated: X-squared = 387.74, p-value < 0.0001, downregulated: X-squared = 153.5, p-value < 0.0001) (Fig. 4). Only five contigs were up-regulated in the Fed vs T₀ comparison whilst down regulated in the FDP versus T₀ comparison. No contigs were down-regulated in the Fed vs T₀ and at the same time up-regulated in the FDP vs T₀ comparison (Fig. 4).

GO enrichment analysis for biological processes (BP) was employed to study the biological functions of the differentially expressed contigs. This would depict global effects of season and feed deprivation during summer on the charr’s brain transcriptome.

Contigs that were only up-regulated in fed fish over time were associated with biological processes that included oxygen transport and protein related biological processes (Table 3). Down-regulated contigs in fed fish contributed amongst others to terms like ion transport, protein complex assembly (Table 4), and “feeding behaviour” (Additional file 4: Table S2). In both fed and feed deprived charr, we found that up-regulated contigs over time were involved in basic cell metabolism processes such as DNA replication and RNA metabolism (Table 5). Contigs that were down-regulated over time regardless of the feeding regime were amongst others related to ion transport, protein related processes and wnt signaling (Table 6, Additional file 6: Table S4). Feed deprivation over time led to an up-regulation of contigs relating to processes such as catabolism, apoptosis and immune system (Table 7) and a down-regulation of e.g. oxygen transport (Table 8). An up-regulation of contigs that take part in catabolic and apoptotic processes and down-regulation of oxygen transport was also found in the endpoint comparison between FDP and fed charr (Additional file 9: Table S7, Table 9, Additional file 10: Table S8).

Table 3 Biological processes enriched by up-regulated contigs only found in Fed versus T₀ (see yellow fraction in Venn diagram Additional file 1: Fig. S1). Terms sorted by the number of contributing contigs

Table 4 Biological processes enriched by down-regulated contigs only found in Fed versus T₀ (see yellow fraction in Venn diagram Additional file 1: Figure S1) Terms sorted by the number of contributing contigs

Table 5 Biological processes enriched by up-regulated contigs found in Fed versus T₀ and FDP versus T₀ (see white fraction in Venn diagram Additional file 1: Figure S1). Terms sorted by the number of contributing contigs

Table 6 Biological processes enriched by down-regulated contigs found in Fed versus T₀ and FDP versus T₀ (see white fraction in Venn diagram Additional file 1: Figure S1). Terms sorted by the number of contributing contigs

Table 7 Biological processes enriched by up-regulated contigs only found in FDP versus T₀ (see blue fraction in Venn diagram Additional file 1: Figure S1). Terms sorted by the number of contributing contigs

Table 8 Biological processes enriched by down-regulated contigs only found FDP versus T₀ (see blue fraction in Venn diagram Additional file 1: Figure S1). Terms sorted by the number of contributing contigs

Table 9 Biological processes enriched by down-regulated contigs comparing feed deprived versus fed charr at end of experiment. Terms sorted by the number of contributing contigs

In an attempt to unravel whether central appetite signalling pathways in the charr brain transcriptome were modulated by feed-deprivation, we screened the lists of differentially expressed contigs systematically for candidate genes that have previously been demonstrated to be involved in the regulation of appetite and energy homeostasis in fish [9, 10]. Further, due to a strong effect of season and temperature seen on the brain transcriptome (Fig. 3), genes involved in seasonal rhythms were included in the search. Lastly, we searched the dataset for possible new actors involved in the regulation of energy homeostasis and food intake, not previously described in fish.

Several candidate genes possibly involved in seasonality and food intake control in fish were found to be differentially expressed between one or several comparisons and are displayed in Table 10. Differences in expression were mainly found between the T₁ and T₀ for FDP and/or fed fish. Deiodinase 2b (Dio2b) was also two-fold lower expressed in FDP charr than in fed charr at the end of the experiment (Fig. 5). The lists of top annotated genes that were found differentially expressed in the different comparisons (Additional file 11: Table S9, Additional file 12: Table S10, Additional file 13: Table S11, Additional file 14: Table S12, Additional file 15: Table S13, and Additional file 16: Table S14) were screened for genes that, based on prior knowledge from mammals, are known to be involved in energy homeostasis and control of food intake. Here, the gene nerve growth factor inducible (VGF) was found as a new candidate that may be involved the regulation of appetite and energy metabolism in fish (Fig. 6).

Table 10 Differentially expressed candidate genes involved in food intake control and seasonality

Fig. 5

Log counts per million reads (logCPM) for contigs corresponding to Dio2b in the three treatment groups. Data are presented as box and whisker plots with median, 25th and 75th percentiles and 1.5 * interquartile range. In addition, individual data points are indicated within the plot

Fig. 6

Log counts per million reads (logCPM) for the contig vgf3.3 in the three treatment groups. Data are presented as box and whisker plots with median, 25th and 75th percentiles and 1.5 * interquartile range. In addition, individual data points are indicated within the plot

Data release

Our raw RNA-seq data are available in the SRA database [24] under accession number SRP151983. In addition, we specifically developed an interface for the present de novo transcriptome. The entire dataset can be browsed at the following URL: http://ngspipelines2.toulouse.inra.fr:9002/ the main menu enables the user to enter the contigs, SNPs and download sections. The contig section main page displays three blocks providing (1) statistics upon the assembly and the annotations, (2) information and access to statistical analysis tools (digital differential display, Venn diagram) on the library and (3) a table that displays data of the contigs that are of the user’s strongest interest. Once a contig has been selected, the user can view its general annotation including Gene Ontology, check its sequence for possible open reading frames, visualize the annotation location in the “jbrowse view” and scan the alignment coverage along the contig for each library. All these features are accessible through the menu located at the top of the page. The SNP section has also a main page presenting some statistics and the favourite table. Once a variable locus is selected, the user can access general information such as the list of alleles and the flanking sequences, allele information with the allele count graph and table for each library and the corresponding functional annotation.

Discussion

Feed deprived and fed charr showed a diverging development of weight and condition

The strong increase in K and body mass observed over the 4 weeks period in the fed charr (Fig. 1) was consistent with the high food intake and growth of anadromous Arctic charr during summer [17, 25]. The FDP charr underwent a strong mobilization of energy reserves during the experiment resulting in a markedly lower K and body mass in FDP charr than in fed charr at T₁ (Fig. 1). This effect of treatment was also evident in the 5 fish from each treatment sampled for transcriptome analyses.

A reliable de novo brain transcriptome

Out of the 49829 generated contigs, more than 85% were annotated, mainly based on Atlantic salmon (Fig. 2). This was expected as Atlantic salmon was phylogenetically the closest related species with a sequenced genome at the time the charr de novo transcriptome was built (spring 2016) [26]. Besides other teleost species contributing to the annotation, Homo sapiens also gave hits: the human genome is well characterized and might have brought annotations of genes that, so far, have only been described in humans.

Even if it is known for its complex transcriptomic signature, the brain alone cannot be representative of the entire diversity of the Arctic charr transcriptome. Thus, a completeness of 67% (BUSCO analysis) confirmed that the assembly produced correct contigs. Indeed other fish transcriptome de novo assemblies brought results in the same order of magnitude with 70.2% completeness for the gut tench (Tinca tinca) [27] and 64% for 4 combined tissues of the noble crayfish (Astacus astacus) [28]. In addition, (1) the comparison to the phylogenetically close Atlantic salmon reference proteome, (2) the calculation of the realignment rates and (3) the mapping of an Arctic charr public raw data using a de novo assembly from short read RNA-seq data on our dataset also confirmed a high completeness of the built transcriptome.

Global brain gene expression patterns differ over time while feed deprivation has a moderate impact

The correlation heatmap highlighted that charr brain transcriptomes from T₀ clustered together and were markedly different from the ones sampled at T₁. Samples of the fed and FDP charr at T₁ did not cluster in accordance with the treatment group. These results stand in contrast with the strong divergence in K and weight between the fed and FDP charr. However, feed deprivation may only have affected very specific processes in the brain transcriptome leaving the overall expression pattern less strongly affected.

Seasonal changes alone (including the increased water temperature), or in combination with feed deprivation, induced a high number of differing contigs from start to end of the experiment, compared to the minor differences observed between the treatment groups at the endpoint. This pattern depicts strong shifts in the charr brain gene expression over a 4 weeks period during summer, which occurred regardless of feeding regime (Fig. 4). The results underline that seasonal processes, including changes in water temperature and summer growth, have much stronger effects on the brain transcriptome in the charr than 1 month of feed deprivation during the feeding season. However, the two times higher number of differentially expressed contigs over time in the feed deprived charr than in the fed conspecifics may suggest an enhancement of seasonal differences by feed deprivation (Table 2). The large changes seen over time correspond with an earlier transcriptome study on different Arctic charr morphs, where time, in that case during early development, was the most important factor for differentially expressed transcripts [29].

In addition, the Up-set graphs (Fig.4) showed that a larger fraction of both up- and down-regulated contigs differed only in response to feed deprivation. On the other hand, most of the differences found in the fed fish were at the same time found in the FDP charr. These shared differences over time depict robust seasonal processes that remain unaffected by feed availability.

Food deprivation partially offsets increase in brain metabolism during summer growth

Feeding specific differences over time

Up-regulated contigs over time denoted foremost oxygen transport (genes encoding for several haemoglobin (Hb) subunits) and protein related biological processes (Table 3). Brain Hb mRNA has been found in rodents and humans [30] and it has been suggested that neural haemoglobin may facilitate oxygen transport in neurons [31], but the exact mechanism remains unknown. In addition, it cannot be ruled that the cDNA library may have contained RNA from non-neuronal tissue including blood vessels and veins, e.g. derived from the highly vascularised saccus vasculosus. The metabolic rate of ectotherms is directly linked to ambient temperature. Specifically, the temperature of optimal growth performance of Arctic charr from North-Norway has been shown to be 14 °C [32]. The increase in water temperature by almost 10 °C during our study may have involved an increase in metabolic rate of Arctic charr. Hence, the observed increase in expression of genes encoding oxygen transporters may have occurred in order to meet increased oxygen demands at higher temperatures.

Furthermore, the feeding related terms such as “feeding behaviour” and “negative regulation of appetite” appeared in the list of GO terms of down-regulated contigs. The contig that contributed to these GO terms was the anorexigenic neuropeptide CART, thus suggesting an increase in appetite in the fed charr over time (Additional file 4: Table S2). The decrease in CART expression is discussed in detail in the paragraph on candidate appetite regulators.

Feed deprivation specific differences over time

In the FDP charr, up-regulated contigs were associated with processes such as catabolism, apoptosis, and immune function (Table 7). In contrast, no such trends were found in brain transcriptome analyses of 21-day feed-deprived zebrafish [33]. The finding that contigs relating to apoptosis were increased is puzzling, given the consensus that the brain is well protected from starvation in both mammals [34] and fish [35]. However, in mammals, there is a debate on to what extend feed deprivation may initiate a degeneration of the central nervous system, as different studies have given indication for both absence [36] and presence [37] of autophagy in the brain of feed deprived mice. Further experiments are needed to test whether the observed up-regulation of contigs involved in apoptosis were a sign of neuronal degradation in charr.

Furthermore, we found the GO term “ketone body catabolic process”, comprising the gene encoding 3-oxoacid CoA transferase, to be up-regulated in the feed deprived charr, pointing towards an increase in ketone catabolic activity from start to the end of experiment (Additional file 7: Table S5). This is in line with previous studies on Atlantic salmon and rainbow trout, where ketone bodies were found to serve as an important energy source for the brain when food is absent [38, 39].

In contrast to the fed charr, contigs pertaining to oxygen transport were down-regulated in feed deprived charr (Table 8). This finding is in agreement with the down-regulation of transcripts related to oxygen transport in response to feed deprivation previously seen in Atlantic salmon liver transcriptome [40], rainbow trout liver transcriptome [41] and in zebrafish brain transcriptome [33]. The lower expression of Hb in feed deprived charr may be related to a metabolic suppression in FDP fish in order to save energy when feed is absent. Brain metabolic suppression, indicated by a reduction of glucose oxidation has previously been observed in feed deprived rainbow trout [39]. However, these feed deprived rainbow trout showed a decrease in hexokinase and 6-phosphofruktokinase activities in the brain. In contrast, expression of these glycolytic enzymes did not differ between fed and FDP charr in the present study.

Differences between fed and feed deprived charr at the end of the experiment

At T₁, there was an up-regulation of contigs involved in biological processes related to apoptosis by feed deprivation as well as a down-regulation of biological processes related to oxygen transport (Additional file 9: Table S7, Table 9, Additional file 12: Table S10). These findings further support the possibility that brain metabolic processes may have been partly impaired by feed deprivation, as the changes were both visible over time and between the fed and FDP fish at the end of the experiment.

Differences over time regardless of feeding regime

Over time, there was an up-regulation of contigs involved in biological processes such as DNA replication, RNA metabolism, response to steroid hormones and immune response (Table 5). This may indicate that basic cell metabolic processes such as cell proliferation and neuronal development were positively affected by seasonal growth and temperature. At the same time, there was a down-regulation of contigs affecting biological processes such as ion transport, protein related processes and wnt signalling from start to end of the experiment (Additional file 6: Table S4).

Interestingly, in adult zebrafish, activation and deactivation of wnt signalling in a sequential manner has been shown to accommodate proliferation and differentiation of progenitor cells in the hypothalamus [42]. Furthermore, the finding that wnt signalling in the ARC was stimulated by leptin in mouse [43] and by both leptin and long photoperiod in the seasonal Djungarian hamster (Phodopus sungorus) [44] have triggered a discussion for a role of hypothalamic wnt signalling in the seasonal control of energy balance [45]. Yet, we only found differences in expression for WNT4, and the question of a possible seasonal related function of wnt signalling in charr requires further study.

Taken together, the results from the GO analyses indicate that the enforced energy perturbation by feed deprivation may have affected several metabolic processes in the brain at the mRNA level. Yet, most elementary biological processes, including cell division processes and immune responses, differed similarly over time in both fed and feed deprived charr.

Effect of feed deprivation on candidate genes involved in the regulation of appetite and energy metabolism in fish

Expression patterns of known anorexigenic and orexigenic neuropeptides do not indicate an amplified hunger signalling in feed deprived charr

Despite the lower weight and K in feed deprived charr, expression of central appetite regulators did not reflect the feeding regime of the fish (Table 10). Several genes differed in one or both feeding regimes over time, but none between FDP and fish at T₁. There was a lower expression of the anorexigenic CART at T₁ in both fed (LogFC = − 0.66) and FDP (LogFC = − 0.33) charr than in fed charr at T₀ (Table 10). If CART exhibits an anorexigenic function in charr, the lower expression of CART in both fed and FDP fish over time points towards a seasonal increase in hunger signalling in both treatment groups, but no response to feed deprivation. This corresponds to the lack of responses in hypothalamic CART expression seen in previous long-term feed deprivation studies with rainbow trout (4 months) and Arctic charr (4 weeks) [19, 46]. In contrast, brain CART expression decreased in zebrafish after 3 days feed deprivation [47], in Atlantic salmon after 6 days of feed deprivation [48] and in Atlantic cod (Gadus morhua) after 7 days feed deprivation [49].

Being described as a potent satiety signal in mammals and several fish species, hypothalamic POMCA had, unexpectedly, markedly (LogFC = 5.99) higher expression levels in FDP charr at T₁ than in charr sampled at T₀. In contrast to the results seen for CART, this finding may suggest a decrease in hunger signalling in FDP charr over time. This result corresponds to an increased hypothalamic expression of POMCA1 and POMCB in rainbow trout after 4 months of feed deprivation [46]. In another experiment with rainbow trout, hypothalamic POMCA1 was down-regulated after 28 days of feed deprivation [50]. However, POMCA paralogues could not be distinguished in the present study, and a possible subfunctionalisation of these in the charr, which may be reflected in different responses to feed deprivation, cannot be excluded.

Putatively anorexigenic corticotrophin releasing factor (CRF) expression was lower in FDP charr at T₁ than at T₀ (LogFC − 0.55). Previously, reduced CRF brain expression levels were observed in goldfish (Carassius auratus) after 7 days of feed deprivation [51] whereas no changes were observed after short (7 days)- and long-term (4 weeks) feed deprivation in charr [19] or after long-term feed deprivation in rainbow trout [46].

In rat, the preprotachykinin 1 (PPT) protein, encoded by the TAC1 gene, has been shown to be negatively regulated by ghrelin and high fat diets, and hence is assumed to be involved in regulating adiposity in rodents [52]. In goldfish, the post-prandial increase of hypothalamic expression of μ-PPT has led to the suggestion that μ-PPT may signal satiety [53]. In the present study, brain TAC1 expression was lower expressed in FDP charr at T₁ than at T₀ (LogFC − 0.51).

The function of locally produced LEP in the brain is still a matter of debate, both in mammals and in fish [54,55,56], and results from other studies have so far not provided evidence for a role of central LEP in appetite regulation in fish [18, 56]. This study revealed a higher LEP expression in both fed (LogFC 1.7) and FDP fish (LogFC 1.3) at T₁ compared with T₀ fish. The lack of differences between fed and FDP charr at T₁, despite a profound difference in condition factor (Fig. 1), may indicate that central leptin expression is not linked to adiposity. On the other hand, LEPA1 was found to be more highly expressed in the hypothalamus of hyperphagic charr in July than of anorexic charr in May and January [18].

Previous studies have shown a conserved orexigenic function of NPY in fish [12, 57, 58]. NPY expression was lower in fed fish at T₁ compared to T₀ (LogFC − 0.55) while no difference was seen between FDP andT₀. As such, based on the brain transcriptome, no hunger signalling by up-regulation of NPY could be found in the FDP fish. This finding is in accordance with the lack of responses seen in other feed deprivation studies with in Atlantic cod [49], Atlantic salmon [48], charr [19] and rainbow trout [46].

Apelin is considered another potent orexigenic actor in fish [59, 60]. Our data did not reveal an effect of long-term feed deprivation on apelin expression. Brain expression of apelin receptor (APJA), was, however, higher in fed fish at T₁ than in fed fish at T₀ (LogFC = 0.99) but not different between FDP and fed fish at T₀.

In summary, the results of the present study did not show expected responses to feed deprivation in the expression of candidate genes involved in appetite regulation in fish. This result does not necessarily contradict an appetite regulatory role of these actors in this fish. In a previous study with charr, no differences in the hypothalamic expression of AgRP, MC4R, CRF, NPY, CART, POMCs and LEPR were seen between fed and 4 weeks feed deprived fish, whereas the expression of, CART, MC4R and AgRP changed in feed deprived upon re-feeding or exposure to feed flavour for 1 or 5 h [19]. This indicates that changes in gene expression are more likely to be seen during transition stages than during steady-state situations.

Such paradoxical results may be interpreted as an adaptation in high-latitude fish to save energy by reducing feed searching behaviour when feed is absent.

Genes related to energy metabolism and seasonality

Insulin-like growth factor 1 (IGF1) is key growth regulating hormone in vertebrates, and plasma levels of IGF1 usually correlate positively with growth in fish [61]. Accordingly, it has been shown that plasma IGF1 levels vary proportionally with increases and decreases in feeding rate in Arctic charr [62]. Furthermore, hypothalamic IGF1 expression was reduced by 1 month feed deprivation resulting in a positive correlation also between hypothalamic IGF1 expression and K of fed and feed deprived fish [19]. In the present study, IGF1 expression was, as expected, downregulated in feed deprived charr at T₁ compared to fed charr at T₀. However, no difference was found in central IGF1 expression between fed and feed deprived charr at T₁, despite the huge difference in K between feed deprived and fed charr at T₁ (Fig. 1). This discrepancy in results between studies on Arctic charr may relate to the fact that hypothalamic IGF1 expression was measured in the former study by Striberny and Jørgensen [19], while whole brain IGF1 expression was measured in the present study.

Dio2 converts thyroxin (T4) to the biologically active triiodothyronine (T3) which, in turn, is known as an enhancer of several biological processes and exerts pleiotropic functions in the mammalian brain [63]. In mammals and birds, the increase in day length in spring stimulates hypothalamic Dio2 expression, thereby stimulating a range of processes related to seasonal phenotype transitions, including appetite [64]. Similarly, it was recently shown that brain expression of the paralogue Dio2b was elevated in response to an increased day length in Atlantic salmon [65]. We found a significantly lower hypothalamic Dio2b expression at T₁ than at T₀ in both fed and FDP charr, and a lower expression in FDP charr than in fed charr at T₁ (LogFC = − 1.19/− 0.81) (Fig. 5). Our findings may be interpreted as a general decline in Dio2b expression during late summer, a decline that may have been enhanced by feed deprivation. In support of the latter, in the seasonal Djungarian hamster (Phodopus sungorus), hypothalamic Dio2 expression was reduced in response to fasting induced torpor during summer [66]. Further experiments are needed to characterize the function of Dio2 in seasonal processes, including feeding behaviour, in the highly seasonal Arctic charr.

Nerve growth factor inducible (VGF) - a novel candidate involved in the control of appetite and energy homeostasis in fish?

Interestingly, we found brain VGF (non-acronymic, nerve growth factor inducible) to be up-regulated (LogFC = 0.54) in FDP compared to fed charr at T₁ (Fig. 6). To the best of our knowledge, VGF has not been linked to energy metabolism and appetite control in fish. In mammals, the VGF gene encodes for a 68 kDa protein precursor that is abundantly expressed in the brain, particularly in the hypothalamus. VGF cleaves into several smaller peptides that have been shown to be involved in a multitude of processes including nerve growth upon injury, seasonality, and food intake/energy metabolism [67]. Several studies in rodents have given evidence for a role of VGF in the control of energy metabolism [68,69,70] and food intake [71]. Targeted deletion of VGF produces a lean, small, and hyperactive mouse [69]. In mammals, the function of VGF is complex and not entirely understood. For example, 48 h feed deprivation in mice caused in one study an up-regulation of hypothalamic VGF expression [69], and down-regulation in another study [70]. In mice, VGF derived neuropeptide TLQP-21 increases energy expenditure without affecting expression of POMC/CART and AgRP/NPY, suggesting that TLQP-21 exerts its effects downstream of MC4R signalling [72]. Furthermore, in Siberian hamster, ARC VGF expression was induced by a decrease in photoperiod [73] and reduced by T₃ [74], raising evidence that VGF is involved in the control of seasonal feeding in this species. We found Dio2b expression to be lower in FDP charr than in fed charr at T₁. This indicates a reduced thyroid hormone action, which based on the results from Siberian hamster, could be underlying the increased VGF expression seen in feed deprived charr. The implication of this gene in both appetite regulation and seasonality in mammals makes it of particular interest in the strongly seasonal Arctic charr.

Conclusion

In conclusion, the general gene expression patterns in brain transcriptome of fed and feed deprived charr displayed strong shifts in expression of contigs involved in basic cell metabolic processes over time, and only minor differences were seen in response to feed deprivation. However, these changes during the charr’s summer growth appeared to be enhanced by feed deprivation, indicated by a higher number of differentially expressed contigs over time in feed deprived than in fed charr. A decrease in the expression of haemoglobin subunits together with an increase in expression of genes involved in apoptosis, revealed from GO analysis, may indicate a negative effect of feed deprivation on brain metabolism. This is also supported by a substantially stronger reduction in Dio2b expression in feed deprived than in fed charr from start to end of the experiment. However, the brain is a heterogeneous tissue that consists of many different cell types, including neural and non-neural tissue; consequently, further studies are needed to get a better spatial resolution of where the observed changes occur.

Generally, a note of caution is expressed when measurements only include RNA abundances and not the corresponding protein levels. However, in a recent study in which a combined transcriptomic and proteomic approach were applied it was concluded that transcriptomic analyses indeed can be used to predict protein copy numbers [75]. Despite the marked divergence of body mass and K, no clear hunger signalling was found between fed and feed deprived charr, when searching the lists of differentially expressed contigs for anorexigenic and orexigenic candidate genes known to be involved in appetite regulation in fish. Hence, these results illustrate that even during summer the anadromous charr have a vast ability and flexibility to deal with food deprivation.

Methods

Ethics statement

Fish handling and euthanasia (see below) was performed by a competent person and in accordance with the European Union Regulations concerning the protection and welfare of experimental animals (European directive 91/492/CCE). The experiment was approved by the Norwegian Committee on Ethics in Animal Experimentation (ID 3630).

Feed deprivation experiment and sampling of fish

The charr used in the present study were two-year-old immature offspring of the anadromous Hammerfest strain, originating from wild charr caught in 1984 and since then bred at Tromsø Aquaculture Research Station, where the experiment was carried out. Until the start of the experiment they had been held on natural water temperature and light conditions (transparent roof) and fed a commercial Arctic charr feed (Skretting, Stavanger, Norway) ad libitum by automatic feeders. On June 25 (T₀), 2014, 42 fish were anesthetized in Benzocaine (60 ppm) and tattoo-tagged with Alcian Blue staining dye using a Pan Jet needleless injector (Wright Dental, Dundee, UK). Body mass and fork length were measured and the fish were distributed between two 300 L tanks supplied with fresh water. Another 12 fish were dip netted from the stock tank and euthanized by an overdose of Benzocaine (150 ppm). Body mass and fork length were measured. Subsequently, the fish were decapitated and the belly was cut open for sex determination. In order to rule out potential sex-specific differences, only male fish were sampled for RNA-seq. On a total of 5 fish, brains were dissected out and separated into telencephalon, mesencephalon and hypothalamus. Tissues were stored in 1.5 ml Eppendorf tubes containing 1 ml of RNAlater (ThermoFisher Scientific, MA, USA). Samples were kept at 4 °C for 24 h, and then frozen at − 20 °C until RNA extraction.

On the same day, 220 fish from the stock tank were distributed amongst the two tanks in which the tattooed fish had been placed (130 fish per tank). From then on, the fish in one tank were fed (Fed) with the same commercial feed as before. Fed fish were fed two main meals at 08.00 AM and 3.00 PM by automatic feeders and in addition, between the main meals to ensure excess feed availability. The fish in the other tank were feed deprived (FDP) until the end of the experiment. All fish were held at simulated natural photoperiod (69 °N), which was 24 h light at that time of the year, and both experimental tanks were supplied with fresh water provided by a flow-through system. Ambient natural water temperature was 4.5 °C at the start and 13.5 °C at the end of the experiment. On July 23 (T₁), 12 fish from each group were euthanized, from which 5 males were measured and sampled as described above. Finally, the tagged fish were anesthetized and measured for body mass and length.

The high number of fish in each treatment group compared to the number of fish sampled was justified by the need to avoid formation of social hierarchies in the fed group. The tagged fish were included in order to monitor the body mass and condition factor development of the fish in the two treatment groups. Fulton’s condition factor (K) was calculated according to Ricker (1975): K = (W × L^− 3) × 100, where W is body mass in g, and L is fork length in cm.

Sample preparation

Tissues were disrupted using TissueLyser II (QIAGEN, Hilden, Germany), and RNA was extracted using the RNeasy Plus Universal Mini Kit (QIAGEN) according to the manufacturer’s protocol. Concentration and purity of RNA were assessed using NanoDrop ND2000c (ThermoFisher Scientific, MA, USA) and when the 260/280 or 260/230 absorbance ratio was below the quality threshold (1.8), samples were further purified using ethanol precipitation. Genomic DNA was removed by treating the RNA with Ambion TURBO DNA-free™ Kit (Life Technologies, CA, USA). In order to obtain a representative view of the main brain areas that have been shown to be involved in central appetite control, 3 μg of RNA of each brain compartment (telencephalon, mesencephalon and hypothalamus) were pooled resulting in a total of 9 μg RNA per brain and individual. Finally, the overall quality of RNA samples was assessed using Bio-rad Experion Bioanalyzer (Bio-rad, CA, USA) and the RQ ranged from 8.7–9.8 indicating high quality of all samples. We did not observe any clustering in MDS plot that would reflect any potential technical variation (Additional file 2: Figure S2). Samples were then shipped on dry-ice to the GenoToul sequencing platform, Toulouse, France for RNA-seq.

cDNA library construction and paired-end RNA-seq

RNA preparation and sequencing were performed at the GenoToul sequencing platform, Toulouse, France. Fifteen RNA libraries were prepared using the TruSeq RNA sample preparation Kit (Illumina, Hayward, CA, USA), involving the following steps. Poly-As containing mRNA were isolated from 3 μg of total RNA. The mRNA was then chemically fragmented. The cleaved RNA fragments were reverse transcribed into first stranded cDNA using random primers, and second strand cDNA was then synthesized. Adaptors were ligated to the end-repaired cDNA, which contributed to fragment selection after the PCR enrichment step. Each library quality was validated measuring sample concentration and fragment size on an Agilent High Sensitivity DNA chip. Sequence hybridization to the flow cell and cluster generation was achieved using a cBot system and the cluster generation kit (Illumina, Hayward, USA). Hundred base pair fragments were sequenced in paired-end for the 15 samples. The samples were multiplexed and Sequenced By Synthesis (SBS) on an a single lane of the eight-lane flow cell of an Illumina HiSeq 2500 sequencer The sequencing lane of the flow cell was screened by a camera, driven by the HiSeq Control Software. Image correction and base calling was performed using the Real Time Analysis (RTA) software.

General statistics, data assembly and annotation

Testing for differences in body weight and K

Data for weight and K were not normally distributed and therefore, statistical testing was conducted on LOG-transformed data. Changes over time were tested using a repeated measures ANOVA. When differences were found, a pairwise comparison applying Bonferroni correction for multiple comparisons was used to determine main effects. Differences between treatment groups were tested using a 2-sample t-test. All statistical testing was done with SYSTAT 13 and figures were drawn using SigmaPlot 13 (both Systat Software, CA, USA). The significance level was set to p < 0.05.

Sequencing data

RNA-Seq data assembly, annotation and quality assessment

Read quality was checked within the ng6 environment [76] using fastQC [77] and Burrows-Wheeler Aligner BWA [78] to search for contamination. The reads were assembled with the Drap pipeline (version 1.7) [79]. The individual sample assemblies were performed with runDrap using Oases with kmers 25, 31, 37, 43, 49, 55, 61, 65, 69. The individual contig files filtered by FPKM (Fragments Per Kilobase per Million mapped reads) over one were then merged with runMeta and filtered again by FPKM over one to produce the reference contig set.

Contigs were annotated searching sequence homologies against the following Ensembl protein databases: blastx [80] Danio rerio, Gadus morhua, Oreochromis niloticus, Oryzias latipes, Takifugu rubripes, Tetraodon nigroviridis, Xiphophorus maculatus, refseq_rna blastn [81]; swissprot blastp [82]; unigene_Takifugu_rubripes.9 (blastn); unigene_Oryzias_latipes.30 (blastn); unigene_Danio_rerio.126 (blastn); NCBI Arctic char ESTs (blastn); the contigs (blastn). Repeats were identified with repeatMasker [83] (version open-4-0-3, with standard parameters) using Repbase database [84]. The GO annotations were extracted from InterproScan [85] (May 2015 version) [86]. The best SwissProt, RefSeq, or Salmo salar NCBI ESTs hit result was used to classify species by best hits contribution.

Different approaches were used to verify the quality of the built contigs. First, the contigs were processed with BUSCO V2 [87] to verify the number of actinopterygii_odb9 reference genes found and their reconstruction state (partial or complete). Then the Salmo salar protein sequences made available by the NCBI (GCF_000233375.1_ICSASG_v2_protein.faa) were aligned with BLAT (standard parameter, version 34) on the 6 frames translated contigs [88]. The alignment was filtered to retain only hits with at least 80% identity and 80% coverage giving the size of the set of well-reconstructed contigs. Finally, the assembly was validated by (1) verifying the realignment rate of the reads of each individual sample on the contigs, (2) mapping the reads of a charr gill transcriptome sequenced in 100 bp paired-end (SRA accession: SRX314607) [89] on the contigs. The Arctic charr genome [23] was released after we had completed our de novo transcriptome and related data analyses. Still, in order to further validate the data set, S. alpinus proteins and RNAs (NCBI ID: 12179) were aligned to the set of generated contigs using blat v. 35 × 1 with respectively -t = dnax -q = prot and standard parameters. The best blast hit was extracted in both cases. The alignments were filtered keeping only hits with at least 80% identity and 80% coverage.

Polymorphism: SNP and microsatellite search

Reads were aligned to the contigs with bwa mem [78]. They were deduplicated with samtools rmdup, then GATK (Version 3.0–0-g6bad1c6) base quality score recalibration was applied [90]. Indel realignment, SNP and INDEL discovery were performed with HaplotypeCaller using standard hard filtering parameters according to GATK Best Practices recommendations [91, 92]. Indels and SNP were independently filtered; 3 per window of 18b with a minimal quality of 30. The micro-satellites discovery was conducted using Tandem Repeats Finder Version 4.04 [89] using the following parameters: 2, 7, 7, 80, 10, 50, 500 -f -d –m ². These analyses will not be interpreted in the current study but were meant to be made accessible on our de novo transcriptome interface, see data mining interface section for URL and details.

Differential expression of contigs

Patterns of gene expression

Data exploration and gene expression analyses were performed using various packages implemented in R version 3.3.1 (2016-06-21). A sample correlation heatmap based on Pearson’s coefficient of correlation was drawn with pheatmap. The number of reads counts per contig was retrieved and each sample normalized, accounting for compositional differences between the libraries (calcNormFactors function EdgeR).

Paired comparisons of treatments groups were performed in EdgeR package version 3.8.6 [93] according to the users’ guide procedure. We identified differentially expressed contigs using a general linear model and a quasi-likelihood F-test, and correcting for false discovery rate (corrected Benjamini and Hochberg p < 0.05). The following contrasts were made (1) Fed fish at T₁ with T₀ fish (T₁_Fed versus T₀),feed deprived fish at T₁ with T₀ fish (T₁_FDP versus T₀) and feed deprived charr at T₁ versus fed charr at T₁ (T₁_ FDP versus T₁_Fed). Only contigs with a logFC greater than 0.5 and smaller − 0.5 were kept for a further comparison using Venn diagrams and GO enrichment analyses. Up-set graphs [94] were drawn using UpsetR to visualize unshared and shared lists of up- and down-regulated contigs between all three comparisons (total of 6 lists). Chi-square tests were run to test whether there was a relationship of feeding regime and proportions of uniquely versus shared up- and downregulated contigs (2 × 2 contingency tables) of over-time-comparisons. Lists of up- and down-regulated contigs of the T₁_Fed versus T₀ and T₁_FDP versus T₀, comparisons were further compared in the JVenn interface [95] in order to identify (1) contigs that were differentially expressed only in fed charr, (2) contigs that were only differentially expressed in FDP charr and (3) contigs that differed over time regardless of diet, i.e. a seasonal effect independent of the feeding regime (the shared portion of the Venn diagram). Resulting Venn diagrams where drawn in Venn Diagram Plotter [96] (Additional file 1: Figure S1). All lists of differentially expressed genes generated from the initial EdgeR comparison (no log fold change cut-off) were searched for candidate genes known to be involved in appetite regulation and energy metabolism.

GO enrichment

Gene ontology (GO) term enrichment was obtained using the TopGo package [97]. It consisted of the identification of terms that host more differentially expressed contigs than expected by chance in a specific comparison. Enrichment of terms by differentially expressed contigs was assessed using Fisher’s exact test (p < 0.05). This analysis focused on Biological Processes. Only terms that were enriched by 3 or more contigs were presented in the results.

Each of the three gene lists generated by Venn diagrams (Additional file 1: Figure S1) were investigated for GO enrichment, and lists of up- and down-regulated contigs were analyzed separately. Finally, GO enrichment analysis was conducted with the lists of up- and down-regulated contigs of the endpoint comparison.

Data mining interface

The assembled contigs were annotated using the RNA-seq de novo ngs-pipelines processing chain [98] and the results have been uploaded to a web-based user interface build upon biomart [99].