The transcriptome of metamorphosing flatfish

Abstract

Background

Flatfish metamorphosis denotes the extraordinary transformation of a symmetric pelagic larva into an asymmetric benthic juvenile. Metamorphosis in vertebrates is driven by thyroid hormones (THs), but how they orchestrate the cellular, morphological and functional modifications associated with maturation to juvenile/adult states in flatfish is an enigma. Since THs act via thyroid receptors that are ligand activated transcription factors, we hypothesized that the maturation of tissues during metamorphosis should be preceded by significant modifications in the transcriptome. Targeting the unique metamorphosis of flatfish and taking advantage of the large size of Atlantic halibut (Hippoglossus hippoglossus) larvae, we determined the molecular basis of TH action using RNA sequencing.

Results

De novo assembly of sequences for larval head, skin and gastrointestinal tract (GI-tract) yielded 90,676, 65,530 and 38,426 contigs, respectively. More than 57 % of the assembled sequences were successfully annotated using a multi-step Blast approach. A unique set of biological processes and candidate genes were identified specifically associated with changes in morphology and function of the head, skin and GI-tract. Transcriptome dynamics during metamorphosis were mapped with SOLiD sequencing of whole larvae and revealed greater than 8,000 differentially expressed (DE) genes significantly (p < 0.05) up- or down-regulated in comparison with the juvenile stage. Candidate transcripts quantified by SOLiD and qPCR analysis were significantly (r = 0.843; p < 0.05) correlated. The majority (98 %) of DE genes during metamorphosis were not TH-responsive. TH-responsive transcripts clustered into 6 groups based on their expression pattern during metamorphosis and the majority of the 145 DE TH-responsive genes were down-regulated.

Conclusions

A transcriptome resource has been generated for metamorphosing Atlantic halibut and over 8,000 DE transcripts per stage were identified. Unique sets of biological processes and candidate genes were associated with changes in the head, skin and GI-tract during metamorphosis. A small proportion of DE transcripts were TH-responsive, suggesting that they trigger gene networks, signalling cascades and transcription factors, leading to the overt changes in tissue occurring during metamorphosis.

Background

Metamorphosis describes the “change in form” associated with the transition between life cycle stages in a wide range of animal taxa [1–8]. This transition can be accompanied by modifications in morphology, physiology, behavior, habitat and feeding mode. The endocrine system and in particular the thyroid hormones (THs), thyroxin (T4) and triiodothyronine (T3), play a central role in vertebrate metamorphosis acting as transcription factors (TFs) when they bind to their receptors. In amphibian metamorphosis, it is well established that THs directly or indirectly stimulate apoptosis and resorption of larval tissue and also promote growth, differentiation and remodelling of tissues that are crucial for the adult life form [8]. For instance, THs are involved in the control of muscle fibre apoptosis in the amphibian tail during its regression and also promote development of the limbs [9]. In amphibians the change in feeding habit from herbivore to carnivore during the transition from tadpole to frog is associated with TH driven remodelling of the intestine that changes from a long coiled tube into a complex differentiated organ [10–12]. Similarly, THs modulate the change in the amphibian integument from a simple to a stratified structure that is better adapted to terrestrial life [13].

Teleost fish also undergo a TH driven metamorphosis that marks the larval to juvenile transition [14, 15]. However, the term is more generally applied to the profound modifications associated with the change from bilateral symmetry to asymmetry during the larval-juvenile transition of flatfish (pleuronectiformes) [7, 16–20]. In flatfish metamorphosis the external morphology is dramatically transformed and they change from symmetric pelagic larvae to asymmetric benthic juveniles with both eyes on the upper, ocular side of the head (reviewed in Power et al. [7]). The external transformation in flatfish morphology is accompanied by a plethora of changes in the structure and function of tissues and organs. Chemical disruption of the thyroid axis using thiourea or methimazole (MMI) delays or stops stomach development in the Japanese flounder (Paralichthys olivaceus) [21] and otolith mineralization in the Southern flounder (Paralichthys lethostigma) [22]. The importance of THs during metamorphosis is further emphasized in other flatfish species where they have been shown to be important for the maturation of the muscle, stomach and skin [23–29].

Flatfish have a high economic value and include species such as the Atlantic halibut (Hippoglossus hippoglossus), common sole (Solea solea), Senegalese sole (S. Senegalesis), turbot (Scophthalmus maximus) and the half-smooth tongue sole (Cynoglossus semilaevis). Overfishing and high consumer demand for flatfish has made them an interesting target for aquaculture production and a better understanding of metamorphosis is of direct relevance for their efficient and successful production. Specific problems linked to failures during metamorphosis include feeding difficulties, reduced growth rate, arrested metamorphosis, abnormal pigmentation (albinism, ambicoloration or mosaicism), failed migration of the eye and skeletal deformities (reviewed in Power et al. [7]). Control of hatchery production requires an understanding of fish biology but also a comprehension of the mediators of metamorphosis, such as the THs and potentially other endocrine factors. Although there are numerous studies of flatfish metamorphosis, the endocrine and molecular basis of the tissue-specific modifications and the timing of the cascade of events that lead to metamorphosis are still largely unknown. Moreover, experiments that have blocked the thyroid axis during metamorphosis with drugs such as MMI do not significantly modify larval viability suggesting thyroid dependent [22, 30–32] and independent processes underpin this event. A complex task now lies ahead in establishing which developmental processes during metamorphosis are fully TH dependent and which genetic pathways and endocrine systems cross-talk with THs. Furthermore, it remains to be established how profoundly different processes such as skin maturation, eye migration and craniofacial remodelling or gastrointestinal tract (GI-tract) development can be regulated by the same endocrine factors.

One of the challenges of studying metamorphosis in fish larvae is their relatively small dimension, which means pools of larvae rather than individuals or tissues have generally been used which significantly reduces the resolution of such studies. The advantage of the biggest of flatfish, the Atlantic halibut, is the large size of the larvae and their slow metamorphosis (occurring over approx. 58 days), which means it is possible to analyze individuals or individual tissues. This is advantageous as pools of larvae contain a mixture of tissues and frequently contain individuals at different developmental stages making resolution of tissue specific changes in transcripts and proteins during metamorphosis difficult or impossible.

The working hypothesis of the present study is that since THs exert their action by binding to thyroid receptors that are ligand activated TFs the overt change in flatfish during metamorphosis will be preceded by significant modifications in the transcriptome of responsive tissues. For this reason, large scale analysis of tissue-specific transcriptional changes in responsive tissue should provide insight into the underlying molecular changes of tissue specific maturation. A 454 pyrosequencing approach was used to survey the tissue specific transcriptomes in the skin, GI-tract and head of metamorphosing Atlantic halibut and to also generate a reference transcriptome. SOLiD technology was then used to map the transcriptional changes in individuals (n = 3/stage) at different stages of metamorphosis. Differentially expressed (DE) transcripts during metamorphosis were identified by comparing the transcriptome at metamorphic stages (stage 7, 8 and 9 [33]) with juvenile (benthic) stages. Subsequently, genes of the thyroid axis, TH-responsive transcripts and candidate genes that underpin remodelling and maturation of tissues during metamorphosis were identified and analyzed by quantitative PCR (qPCR).

Results

454 transcriptome sequencing

Transcriptome annotation

In spite of stringent quality control of the RNA used for 454 library construction, the number of reads resulting from the stage-specific libraries for skin, GI-tract and head of Atlantic halibut post-embryonic larvae at different metamorphic stages was highly variable. A total of 134 Mbp were produced for the tissue assemblies using MIRA V3 (http://sourceforge.net/projects/mira-assembler/files/) and they assembled into 65,530, 38,426 and 90,676, contigs for skin, GI-tract, and head, respectively. The contigs from the skin, GI-tract and head tissue assemblies were submitted to an iterative stringent four-step local Blast approach (Additional file 1). The tissue assemblies were successfully annotated and for the GI-tract library 60 % of the initial contigs had a good Blast match after the first 3 annotation steps and 57 % of the contigs were annotated for the skin and head libraries (Table 1). Most of the contigs were successfully annotated in the first step of Blastx against the zebrafish refseq protein db.

Table 1 Summary of the Blast results used for annotation of the head, skin and GI-tract assembly

Gene Ontology and KEGG analysis

The active transcriptome in each tissue analyzed was assumed to be equivalent to the number of contigs identified (Fig. 1a). Comparative analyses between the transcriptomes revealed that 2,541, 2,261 and 8,359 transcripts were unique to the skin, GI-tract and head assembly, respectively. In addition, 4,099 transcripts were common between the three tissues. The head and skin, head and GI-tract and GI-tract and skin shared a further 4,464, 955 and 506 transcripts, respectively (Fig. 1a).

Fig. 1

Atlantic halibut skin, gastrointestinal tract and head transcriptome annotation. a Venn diagram of common and unique tissue transcripts (using the transcript name); b Venn diagram representing the common and specific tissue gene ontology (GO) terms (using the unique GO terms); c Diagram representing the relative abundance of shared and tissue specific enriched GO terms by GO category (using the over/under-represented GO terms from the Fisher’s exact test)

The complete functional (GO) annotation for the skin, GI-tract and head transcriptomes (Additional files 2, 3 and 4, respectively) assigned a total of 12,577 different GO terms to the three tissue assemblies. Comparisons of the assigned GO terms for the skin, GI-tract and head transcriptome assemblies revealed 8,253 common GO terms and 293 identified only in the skin, 206 only in the GI-tract and 1,170 only in the head (Fig. 1b). The most abundant GO terms (level 2) for biological process (BP), molecular function (MF) and cellular component (CC) showed no major differences between skin, GI-tract and head (Additional file 5). Cellular process (17 %), metabolic process (16 %) and biological regulation (10 %) were the most representative GO terms in the category biological processes. Other key biological processes linked to development (7 %), localization (6 %), signaling (5 %), cell proliferation (3 %) and death (3 %) were also found (Additional file 5).

For molecular function, the most abundant GO terms (level 2) in the skin and head transcriptomes were binding, catalytic activity and structural molecule activity (Additional file 5). For the GI-tract transcriptome, binding and nucleic acid binding transcription factor activity were the most frequent MF GO terms. The exception was the GI-tract for which GO terms specific for DNA binding transcription factor activity (~20 %) were more highly represented, when compared to the other tissues (Additional file 5). Within the cellular component category, the most represented GO terms in the three tissues transcriptomes were cell, organelle, macromolecule complex and membrane (Additional file 5).

Fisher’s exact tests were applied to detect significantly over/under-represented GO terms resulting from analysis of the tissue-specific transcriptomes (FDR < 0.05). Comparison of the overall GO enrichment for the tissue specific transcriptome revealed that the highest enrichment was associated with the BP category (approx. 60–73 %), followed by the MF category (16–31 %). Of the 2,316 enriched GO terms identified for the tissue transcriptomes, 8.2 % (190), 9.2 % (214) and 82.6 % (1,912) were from skin, GI-tract, and head, respectively. Figure 1c shows the representative enriched GO terms for each tissue assembly for BP, MF, and CC gene ontology categories. The halibut head transcriptome had the highest enrichment of BP GO terms (1,395) and they corresponded to 16.7 % of the overall GO terms for this tissue, followed by the GI-tract (148 GO terms) and the skin (114 GO terms) (Fig. 1c). Similarly, MF GO terms were most enriched in the head (303 GO terms and 9.9 % of overall enriched terms) relative to the skin (58 GO terms) and the GI-tract (51 GO terms).

In the skin transcriptome, GO terms related to the muscle system, development and morphogenetic processes, including epidermis development, appendage morphogenesis, and transcripts involved in cellular response to hormone stimulus were overrepresented along with immune development and pigmentation (Additional file 6). Significant GO categories in the GI-tract that were overrepresented included digestion, proteolysis and lipid metabolism such as the cholesterol metabolic process and triglyceride mobilization (Additional file 7). In the head transcriptome, significantly overrepresented GO terms included development of the nervous system, spinal cord and otoliths, cartilage and endochondral bone. In addition, head-specific GO terms, such as pituitary gland development and thyroid hormone metabolic processes were also significantly overrepresented (Additional file 8).

REVIGO clustering of enriched Biological Process GO terms for the skin transcriptome identified phosphocreatine metabolism and response to lipopolysaccharides and biotic stimulus. In REVIGO clustering of enriched BP GO terms in the GI-tract transcriptome identified digestion, cell proliferation, rRNA transcription, lipid storage and immune system and response (e.g. foam cell differentiation, positive regulation of macrophage derived foam cell differentiation, low density lipoprotein particle remodeling, lipoproteins transport). In the head, transcriptome REVIGO clustering of enriched BP GO categories identified nervous system development (e.g. glutamate receptor and neuropeptide signaling pathways, proliferation and apoptosis of neural precursor cells, regulation of synaptic plasticity and synaptic vesicle transport), blood vessel morphogenesis, and immune and defense response (T cell activation).

More than 125 metabolic pathways were identified via KEGG mapping comprising ~1,250 different enzyme codes and more than 9,000 of the Atlantic halibut contigs matched an enzyme code (EC) (Additional file 9). Overall, there was considerable similarity in the metabolic pathway enrichment between the three tissues but this was unsurprising since many of the pathways were linked to cellular metabolism. Of note was the lower representation of sphingolipids, inositol lipid and phospholipid pathways in the GI-tract. The GI-tract, a soft tissue, had a notable reduction compared to the skin and head of metabolic pathways involved in chondrogenic matrix generation but had an increase in starch and sucrose metabolic pathways relative to the head (Additional file 9).

Tissue development/morphogenesis – identification of putative tissue-specific genes

Blast of the skin transcriptome against the in-house skin-specific database that contained genes characteristic of skin in other vertebrates identified 33 transcripts for skin development and morphogenesis and 40 for pigmentation (Table 2). Abundant transcripts included the collagens (col1a1, col1a2), genes involved in pigmentation, melanocyte differentiation and melanosome transport (e.g. apoptosis regulator bax, dedicator of cytokinesis 7, dopachrome tautomerase, lysosomal-trafficking regulator), (Table 2). Several signal transduction pathways were identified in the skin transcriptome including Notch (>70 transcripts), Wnt (>100 transcripts) and Sonic Hedgehog (Shh) (>30 transcripts) (Additional file 2).

Table 2 Genes in vertebrate skin development and pigmentation identified in the Atlantic halibut skin transcriptome

When the GI-tract transcriptome was compared with the in-house GI-tract specific database, 72 genes were identified. Identified gene transcripts included those with sequence similarity to signal transduction pathways, such as Sonic Hedgehog (Shh), Wnt and bone morphogenic protein (Bmp), as well as genes involved in gastric function (Table 3, Additional file 3).

Table 3 Genes of the vertebrate digestive system identified in the Atlantic halibut GI-tract transcriptome

The head transcriptome contained genes involved in thyroid gland development and thyroid hormone physiology, such as enzymes involved in the activation or inactivation of THs (DIO1, DIO2, DIO3), TH receptors (TRαB and TRb) and other nuclear receptors (Table 4). The pigmentation genes identified in head were coincident with those found in skin. Signalling pathways associated with development were well represented and included in the head transcriptome: i) the Wnt receptor signalling pathway (including casein kinases, low-density lipoprotein receptor related proteins, frizzled-related proteins, spontins, wnt11, wnt7, wnt8, wnt9); ii) the transforming growth factor beta (TGFβ) receptor signalling pathway (including activins, bone morphogenetic proteins, collagens, latent-transforming growth factor beta-binding proteins and tgf-beta receptors); iii) the Notch signalling pathway; iv) the Hippo signalling cascade, and v) the Hedgehog signalling pathway (Additional file 4).

Table 4 Genes involved in the TH axis identified in the Atlantic halibut head, skin and GI-tract transcriptomes

Identification of thyroid hormone responsive genes

Overall, 135 putative TH-responsive genes were identified in the skin, GI-tract and head transcriptome, which included TFs, genes involved in DNA replication, cell proliferation, cell growth and differentiation and collagen synthesis and degradation (Fig. 2). The skin transcriptome contained 113 putative TH-responsive genes that mainly corresponded to structural proteins, proteases, actins, transmembrane proteins, and several membrane transport proteins of the solute carrier group (Fig. 2). The GI-tract transcriptome contained 62 putative TH-responsive genes and included genes involved in DNA replication, cell cycle and metabolic pathways (Fig. 2). The head transcriptome was enriched with 99 putative TH-responsive genes of which 11 were specific to the head transcriptome and included TFs, DNA replication and ion binding proteins (Fig. 2).

Fig. 2

Putative thyroid hormone (TH) responsive genes identified in the Atlantic halibut transcriptomes. Heat map of putative thyroid hormone (TH) responsive genes identified in the transcriptomes of skin, GI-tract and head of metamorphosing Atlantic halibut. The acronyms of TH-responsive genes are indicated and full names are given in Additional file 18

SOLiD transcriptome comparison between metamorphic stage transitions

Identification of differentially expressed transcripts during metamorphosis

Pairwise comparisons of whole larvae transcriptomes between metamorphic stages generated a very low number of DE genes (Additional file 10f–j). For the premetamorphic stage 5, 4,155 transcripts were DE when the transcriptome was compared with the transcriptome for the juvenile stage. In contrast, more than 8,000 transcripts were DE when the transcriptome of whole larvae of each metamorphic stage was compared with the transcriptome of whole juveniles. The majority of the 8,000 DE transcripts per stage were up-regulated in the juvenile stage relative to the metamorphic stages (Fig. 3) and 3,336 of the DE transcripts were common between the metamorphic stages (7, 8, 9A, 9B and 9C). The number of DE transcripts specific to each stage was 403, 365, 1,214, 446 and 362 for stages 7, 8, 9A, 9B and 9C, respectively. The number of DE transcripts common between stages 7 and 8 was 5,999, between stages 8 and 9A was 5,638, between stages 9A and 9B was 5,272 and between 9B and 9C was 5,972.

Fig. 3

Differentially expressed transcripts between Atlantic halibut metamorphic stages and juveniles. Graphical representation of the relative number of DE transcripts (up- and down-regulated) identified when pro-metamorphic (stage 7), proclimax-metamorphic (stage 8) and metamorphic climax (9A, 9B and 9C) Atlantic halibut are compared with the juvenile post-metamorphic stage

Expression of thyroid related transcripts during metamorphosis

TH-responsive transcripts DE with SOLiD transcriptional profiling were identified by filtering all differential transcripts using the “in house” database. Overall, 145 putative TH-responsive transcripts were DE (log2 of the fold change of juvenile versus all metamorphic stages), (Fig. 4, detailed information regarding transcripts in Additional files 11 and 12). The majority of the putative TH-responsive transcripts were down-regulated in the metamorphic stages relative to the juvenile (Fig. 4a). The exception was stages 7 and 8 that had 10 and 2 up-regulated putative TH-responsive transcripts, respectively (Additional file 11). Comparison of the putative TH-responsive genes in each stage revealed 41 that were common. Stage 8 had the greatest number of putative TH-responsive transcripts (119), followed by stage 9B (101), stage 9A (98), stage 9C (96) and then stage 7 (85) (Fig. 4c). Enriched reactomes of TH-responsive transcripts during metamorphosis included cellular response to stress, the cell cycle, DNA repair, DNA replication, apoptosis, metabolism, and signal transduction (Additional file 13).

Fig. 4

Putative TH-responsive transcripts with differential expression between Atlantic halibut metamorphic stages and juveniles. Clustering of the putative thyroid hormone (TH) responsive transcripts with differential expression between metamorphic stages and juveniles of Atlantic halibut. a Heat map of the DE TH-responsive transcripts clustered by expression pattern. Transcript expression is represented as log2 of fold change for metamorphic stages versus juveniles; b Venn diagram revealing the number of DE TH-responsive transcripts that are shared between stages or that have a stage specific expression

Up-regulation of TH axis related genes during metamorphic climax

Transcripts that were not DE in SOLiD analysis (presumably due to methodological limitations) but that are involved in the thyroid axis, such as, TH production (thyroglobulin - Tg), transport (monocarboxylated transporter 8 - MCT8, monocarboxylated transporter 10 - MCT10), metabolism (deiodinase 1 - DIO1, deiodinase 2 - DIO2, deiodinase 3 - DIO3) and action (thyroid hormone receptor alpha A - TRαA, thyroid hormone receptor alpha B - TRαB, thyroid hormone receptor beta - TRβ) were analyzed by qPCR using the same samples used for SOLiD analysis (Fig. 5, Additional file 14). Tg transcript levels were lower in the pre-metamorphic stage (stage 5), significantly (p < 0.05) higher during metamorphosis, and then decreased significantly (p < 0.05) in the post-metamorphic stage (Additional file 14). TRs had a variable expression during metamorphosis and the relative transcript abundance was TRβ > TRαB > TRαA (Fig. 5a). The transcript abundance of all the TRs increased significantly (p < 0.05) during metamorphosis. TRαA transcript abundance was significantly higher at metamorphic climax (stages 9B and 9C) relative to stages 5, 6 and 7 and TRαB transcript abundance was also significantly (p < 0.05) higher at metamorphic climax (stage 9A) compared to premetamorphic stage 7. The transcript abundance of TRβ throughout metamorphosis (stages, 8, 9A, 9B and 9C) was significantly (p < 0.05) higher than during pre-metamorphosis (stages 5 and 7). The transcript abundance of the three TRs was significantly (p < 0.05) lower in juveniles relative to metamorphic climax (9C), (Fig. 5a, Additional file 14). The relative gene expression of MCT10 was higher than that of MCT8 throughout metamorphosis, although it did not change significantly at any stage (Fig. 5b, Additional file 14). MCT8 gene expression increased significantly (p < 0.05) in stage 9A and 9B relative to stage 8. In juveniles, MCT8 mRNA expression levels were significantly (p < 0.05) lower than in metamorphic stages. The gene expression profile of the three deiodinases (DIO1, DIO2 and DIO3) during metamorphosis was similar, and all were significantly (p < 0.05) up-regulated during metamorphosis (stage 9A–9C) compared to pre-metamorphic (stage 7 and 8) and juvenile stages (Fig. 5b, Additional file 14).

Fig. 5

Expression pattern of transcripts involved in the TH cascade during halibut metamorphosis. Schematic representation of the relative gene expression by qPCR of a thyroid hormones action (TRαA, TRαB, TRβ) and production (Tg); and b thyroid hormones transport (MCT8, MCT10) and regulation of the cellular availability of THs (DIO1, DIO2, DIO3) during Atlantic halibut metamorphosis. Results are presented as relative gene expression (arbitrary units). For detailed information and significance between stages for each transcript, please see Additional file 14

Confirmation of differentially expressed transcripts in SOLiD by qPCR

Six transcripts from SOLiD analysis of metamorphosing Atlantic halibut were also analyzed by qPCR (Additional file 15) and had a concordant expression pattern (Fig. 6). Transcripts with unchanged transcript abundance during metamorphosis in SOLiD: ribosomal protein L7 (RPL7) and 40S ribosomal protein S30 (FAU) were not significantly different in qPCR. Transcripts, alpha-globin 1 (Gloα1), carboxypeptidase A2 (Cpa2), apolipoprotein AI (ApoAI) and type I keratin isofom 2 (Krt1i2), significantly modified in SOLiD analysis were also significantly (p < 0.05) modified in the qPCR results during metamorphosis. A high and significant positive correlation (r = 0.843; p = 1.07 x 10⁻⁷) was obtained when results of SOLiD analysis for six genes (Gloα1, Cpa2, ApoAI, Krt1i2, RPL7, FAU) were compared with qPCR expression levels (relative to the geometric mean of 40S ribosomal protein S4 and Elongation factor 1 alpha - RPS4/EFIAI) for metamorphic stages 5 - 9C (Fig. 6a). A lower, but significant positive correlation (r = 0.576; p = 1.23 x 10⁻⁴) was obtained when the juvenile post-metamorphic stage was included in the correlation analysis between SOLiD and qPCR data (Fig. 6b).

Fig. 6

Correlation analysis between SOLiD and qPCR expression analysis. Correlation analyses between SOLiD and qPCR expression of transcripts with a constant expression and transcripts with a modified expression during Atlantic halibut metamorphosis. Comparison of normalized counts (SOLiD data) and relative gene expression profile (qPCR data) of six genes: apolipoprotein A-I (ApoAI), carboxypeptidase A2 (Cpa2), 40S ribosomal protein S30 (FAU), alpha-globin 1 (Gloα1), type I keratin isoform 2 (Krt1i2) and ribosomal protein L7 (RPL7). Different genes are represented by a specific symbol in the graph. Panel a Pearson Product Moment Correlation using metamorphic stages 5 to 9C (r = 0.843; p = 1.07 x 10⁻⁷). Panel b A Pearson Product Moment Correlation using metamorphic stages 5 to juvenile (r = 0.576; p = 1.23 x 10⁻⁴)

Discussion

Changes during flatfish metamorphosis are not limited to modifications in external morphology but include many structural and functional modifications. The role of TRs as ligand activated TFs means that a significant part of the action of THs on tissues is associated with tissue specific modifications in the transcriptome. In the present study, the large size of Atlantic halibut was utilized to establish for the first time in flatfish specific transcriptomes for larval skin, GI-tract and head during metamorphosis using 454 sequencing. SOLiD sequencing of individual larvae (n = 3/stage) and stage specific comparisons (e.g. 7 vs 8; 8 vs 9A, 9A vs 9B and 9B vs 9C) revealed a very low number of DE transcripts and no sudden or dramatic change between any particular stage. In contrast, pairwise comparisons of the juvenile transcriptome with stage specific transcriptomes revealed a high number of DE transcripts (>8,000/stage), with the majority highly up-regulated in the juvenile stage.

Comparisons of the approximately 8,000 differential transcripts per stage generated a stage specific molecular fingerprint. The large majority of DE transcripts (approx. 98 %) were not classified as TH-responsive and presumably represented transcripts underlying ontogenetic changes and belonging to gene networks that lead to the overt changes that accompany metamorphosis. The latter probably explains why blocking TH action with MMI during flatfish metamorphosis is not lethal and only modifies the development of some specific tissues [22]. In line with this observation THs maintain neoteny in only some of the tissues in salamanders [34–36]. The link between TH-responsive pathways and the numerous other gene networks that change during metamorphosis was not established in the present study, but in future studies will be explored. The majority of the putative DE TH-responsive genes clustered in specific metamorphic stages rather than over the duration of metamorphosis. The response of the majority of putative TH-responsive transcripts was none synchronous with the peak in whole body TH levels, previously reported to occur at stage 9 and 10 for Atlantic halibut [37]. The non-synchronous tissue specific response of putative TH-responsive genes during metamorphosis suggests the chronology of tissue responsiveness during metamorphosis may vary, presumably as a result of differences in cellular responsiveness to THs. Such a phenomenon was reported in a recent study of GI-tract development in the Atlantic halibut [28].

Metamorphosis-specific tissue transcriptome

Next-generation pyrosequencing 454 technology has been used to characterize the transcriptome from several flatfish species. In turbot (Scophthalmus maximus) the study focused on immune related transcripts [38], in the common sole a pooled larval and adult liver and GI-tract transcriptome was established [39] and in Senegalese sole and common sole reference transcriptomes were derived by sequencing several tissue from juveniles and adults [40]. An oligo-array study of common sole development from larva to juvenile revealed a large variety of biological processes occurred during development and that some genes of the thyroid axis were associated with the initiation of metamorphosis [39].

To our knowledge, ours is the first next generation sequencing study analyzing individual tissues and larvae of a metamorphosing flatfish. More than 400,000 sequences per specific-tissue were obtained after filtering to remove poor quality sequences and contaminating transcripts (e.g. prey in the GI-tract Artemia sp., etc). This data substantially increases available molecular resources for Atlantic halibut [41, 42]. MIRA3 assembly of the tissue transcriptomes (head: 1,186,541; skin: 830,524 and GI-tract: 418,303) generated 90,676; 65,530 and 38,426 contigs for head, skin and GI-tract, respectively, which was similar to previous studies using the same sequencing strategy [43], but higher than the gene content of the genome of model teleost species e.g. 19,388 for Takifugu rubripes and 31,953 for Danio rerio (www.ensembl.org). Technical issues, read length and the heuristic nature of the assembly methods, no doubt explain the relatively high transcript number as was previously observed in a liver transcriptome study of Zoarces viviparus [44]. Transcript annotation levels in Atlantic halibut were similar to previous 454 studies in turbot, seabream, European eel and silver carp [38, 43, 45, 46]; the contribution of alternative splicing to the high number of assembled transcripts was not established in the study.

Candidate biological processes and pathways during metamorphosis

Knowledge about the mechanisms underlying the global molecular and cellular changes during fish development, including functional gene annotation during metamorphosis, has significantly increased in the last decade due to genomics and transcriptomics technologies. Previous studies using microarrays, expressed sequence tags (ESTs) and candidate genes in Atlantic halibut identified genes involved in muscle, skin, immune system, signal transduction and transcription factor activity in adult and larvae, but information about larvae undergoing metamorphosis is limited [23, 26, 41, 47]. Studies in Solea senegalensis focusing more on metamorphosis generated 10,000 ESTs [48] that are enriched in transcripts involved in the reorganization of somatic tissues, such as, ribosomal proteins, elongation factors and cytoskeletal proteins. The global gene ontology of the individual tissue transcriptomes characterized in the present study, are far more detailed than previous EST studies, but where there is coincident sequence data the results are similar. The GO results of metamorphosing Atlantic halibut is similar to results for other developing teleosts (larval, juvenile and adult) [43, 49], but also between sexes (male and female), whole fish [50] and fish under diverse challenges (e.g. viral challenge) [38], suggesting maintenance of tissue, organ and organism function involves an overwhelming number of common genes that emerge irrespective of the experimental situation. To overcome this problem in the present study we applied a Fisher’s exact test to identify significantly over/under-represented GO terms for the tissue-specific transcriptomes of metamorphosing Atlantic halibut. The 454 transcriptome approach gave insight into tissue specific molecular changes and allied to SOLiD analysis of several individuals/stage revealed core TH-responsive genes responsible for the timing of stage specific responses of individual tissues.

The metamorphosing Atlantic halibut skin transcriptome: is enriched in GO terms related to epidermis and connective tissue development, appendage morphogenesis and pigmentation, which is concordant with the morphological modifications observed [51, 52]. Genes involved in vertebrate skin development and morphogenesis are also enriched and include components of the extracellular matrix [ECM, collagen type I (col1a1, col1a2) and type V (col5a2) [53, 54]], ECM remodeling [discoidin domain receptor 1 (ddr1) [55, 56]] and ECM degradation [stromelysin-3 (MMP11), collagenase 3 (mmp13), matrix metallopeptidase 2 (Mmp2) and metalloproteinase inhibitor 3]. Several of the enriched ECM proteins in skin are responsive to THs [13, 57–60] and in SOLiD analysis their DE pattern in stages is asynchronous presumably because of their association with opposing biological processes. The TH-responsive proteins together with the detected TH-independent growth factors, chemokines, adhesion molecules and proteoglycans have previously been identified in relation to tissue differentiation, development and morphogenesis in vertebrate skin [13, 54, 57, 61].

MMP genes identified for the first time in Atlantic halibut skin may be associated with larval-type cell apoptosis during ECM degradation as previously reported in Xenopus. For example, collagenase 3 (mmp13) has an important role in Xenopus body skin remodelling during metamorphosis [13] and this transcript is present in the halibut larval skin transcriptome and is DE in SOLiD analysis of several individuals/stage with a significant reduction post-metamorphosis. Pathways involved in “focal adhesion” and “tight junction” are also DE in pools of pre-metamorphic Solea solea [39]. The enriched transcripts identified during Atlantic halibut metamorphosis that contribute to pigmentation are of practical interest as abnormal pigmentation can have a significant impact on commercial production of flatfish [7, 20, 62, 63].

The GI-tract in metamorphosing Atlantic halibut: undergoes extensive remodeling to prepare it for the shift in habitat and diet of the juvenile [28, 64, 65]. SOLiD analysis revealed DE of digestive enzymes, such as trypsin, chymotrypsin and phospholipase A2 during Atlantic halibut metamorphosis, as previously reported in other fish species [29, 66, 67]. Our enriched GO results for GI-tract development in Atlantic halibut are similar to that reported in Xenopus, in which GO terms related to digestion are “shut down” at metamorphic climax, but increase again at the end of metamorphosis [68]. The results of the present study corroborate those of a detailed study of GI-tract development in Atlantic halibut that linked up-regulation of pepsinogen and H⁺/K⁺-ATPase α and β subunit with acquisition of a functional proteolytic stomach in early juveniles [28]. However our results diverge from those of an earlier Atlantic halibut microarray study in which genes involved in digestion are more abundant in larvae entering metamorphosis [69], and this may be a consequence of differences in staging, sample composition (pools of larvae were used in previous studies) and the more comprehensive results possible with NGS.

The head transcriptome: The enrichment in the head transcriptome of bone related genes ties in with experiments in Paralichthys lethostigma in which the development and growth of both sacculus and utricle otoliths are TH dependent during metamorphosis [22]. Alpha-tectorin, otolin and plasma membrane calcium ATPase are also enriched in the Atlantic halibut head transcriptome and a previous candidate gene study suggested they are TH-responsive during flatfish metamorphosis [70]. Several TFs specific for thyroid gland development, such as homeobox protein NK2.1, hematopoietically expressed homeobox (Hhex) and Pax8 are enriched in the Atlantic halibut head transcriptome [71], and suggests that modification of the thyroid tissue is essential for successful metamorphosis [17] and disruption of this process may explain failed metamorphosis in some cases.

Enriched pathways in metamorphosing Atlantic halibut: revealed as expected that the essential signaling pathways that trigger tissue development and cell proliferation and differentiation (e.g. Notch, Sonic hedgehog, Wnt, BMP) [72–77], are well represented in all three tissue transcriptomes. This fact lends support to the idea that in Atlantic halibut it is probably not remodeling that gives rise to juvenile tissue but rather de novo tissue development during metamorphosis, as has been demonstrated in Xenopus [78–80]. Several of the signaling pathways are regulated by THs and specific studies will be required to establish their precise mode of action and the tissue specific consequences of their up-regulation.

The majority of DE genes detected by SOLiD analysis of several individual halibut larvae per stage: during metamorphosis are not directly TH-responsive, suggesting that many of the TH effects may be indirect. Cross-referencing of putative TH-responsive genes in whole larvae with the tissue specific transcriptomes provides insight into core tissue changes during metamorphosis. The down-regulation of transcripts linked to the MAPK signalling cascade (c-raf, Kras and C-Jun) during metamorphosis suggests the coordination of TH actions may be via modulation of signalling pathways as has been suggested in Xenopus [81]. Similarly, modification in TFs may be another way in which THs bring about an indirect effect. Thus fos-related antigen-2 (fosl2 or FRA-2), sox4, TFPA2, TGFB, HMG1, CEBPD, GTF2H, NFIX and GTF2F that are all TH-responsive [82–85] peaked at metamorphic climax (stage 9A/9B). This is also reminiscent of what occurs during Xenopus laevis metamorphosis where TFs have a central role in tissue specific TH-induced programs [86–88].

The reliability of the results of SOLiD DE analysis is generally confirmed by comparison to the results of previous candidate gene studies in fish and amphibians. For example, osteonectin (Sparc), that plays an essential role in tissue morphogenesis [89–91] is strongly down-regulated in stage 7 (log2 fold change, −6.4) but its abundance increases at metamorphic climax stage 9A (log2 fold change, −1.6) and this is reminiscent of what occurs in the flatfish Scophthalmus maximus [92]. DE ECM transcripts (alpha2 Collagen type 1 and fibronectin) during Atlantic halibut metamorphosis linked with epidermal outgrowth [93, 94] are also modified in amphibian metamorphosis [86]. The stage specific fingerprint of DE TH-responsive and nonresponsive genes generated by SOLiD analysis of several individuals per stage during Atlantic halibut metamorphosis is a rich resource for future studies of the metamorphic process and its evolution [59]. Furthermore, although in general metamorphosis is comparable between fish and amphibians their divergent evolution, biology and physiology [95] makes flatfish specific data for this process a priority.

Confirmation of the TH axis role in Atlantic halibut metamorphosis

In flatfish, initiation of metamorphosis is associated with a surge in T4 and T3, which increases up until the metamorphic climax and decreases in post-climax stages [17, 37, 96–99]. The failure to detect by NGS analysis DE genes of the thyroid axis in the present and previous studies of flatfish metamorphosis may be a result of: i) their generally low tissue abundance, which is further aggravated by, ii) the dilution effect caused by using mRNA from whole larvae (or pools of larvae) rather than discrete tissue and iii) the asynchronous temporal expression pattern in different tissues. Nonetheless, Tg transcript abundance detected by qPCR in the present study mirrored the TH profiles in metamorphosing Atlantic halibut [37] and is reminiscent of what occurs in Senegalese sole [100]. Unsurprisingly, transcription of deiodinases (DIO1, DIO2, DIO3), that encode selenoproteins that activate and inactivate THs [8, 101–105] changed during metamorphosis. The results for DIO expression agreed with previous studies of metamorphosis in the Atlantic halibut [26] and Senegalese sole [106]. A limitation of the present study is the impossibility of mapping the spatial and temporal pattern of deiodinase mRNA localization, which is known to be tightly controlled in flounder metamorphosis [107]. The spatial distribution of deiodinases probably contributes to the asynchronous pattern of DE genes during metamorphosis.

In the Atlantic halibut qPCR revealed that the recently identified TH transmembrane transporters (members of the solute carrier (Slc) proteins [108–110]), MCT8 and MCT10, that regulate TH availability in peripheral tissues, are expressed in metamorphosing Atlantic halibut. However, only MCT8, a specific TH transporter [111–113], is significantly up-regulated during the metamorphic climax (stages 9A and 9B) and significantly decreases in post-metamorphic juveniles. The results in Atlantic halibut suggest that in common with metamorphosis in the amphibian (Xenopus tropicalis) the tissue distribution and abundance of Slc proteins is one of the factors explaining the differential tissue responsiveness to THs [114–118].

In the Atlantic halibut, TRs had a variable expression pattern during metamorphosis as observed in other flatfish species [18, 37, 95, 119] and this is intriguing when placed in the context of the duality model of TR actions during vertebrate development [120]. In this model, TRα is the predominant TR form during the Xenopus larval phase and is associated with repression of TH-inducible genes. Repression of TH-responsive genes occurs when TR and retinoid X receptor (RXR) bind to thyroid response elements (T₃RE) and in the absence of T3 recruit a co-repressor complex (e.g. Nuclear receptor CoRepressor (NCoR), Silencing Mediator for RAR and TR (SMRT), and other proteins). At metamorphosis the presence of T3 leads to substitution of the co-repressor complex by co-activator proteins and TH-responsive gene transcription is induced. This event is concomitant with TRβ up-regulation (reviewed by Grimaldi et al. [121] and Morvan-Dubois et al. [122]). In Atlantic halibut, our results suggest a dual TR activity model may also exist as TRαB is the main TR form expressed in premetamorphic stages, while TRβ is more abundant at the metamorphic climax. In addition, SOLiD analysis reveals DE of co-activator and repressor elements (NCoR, HDAC1, PRMT1) of TRs during metamorphosis.

In summary, although significant changes in transcripts of the thyroid axis are not detected using SOLiD transcriptome analysis, the temporal expression pattern of DIO, TH transporters and TRs varied dramatically between larvae at different stages confirming the importance of the TH axis in Atlantic halibut metamorphosis [99]. In metamorphosing frogs changes in DIOs, TH transporters and TRs are correlated with the timing of tissue specific changes during metamorphosis [10, 77, 103, 114, 121, 122]. The results of our study and those on frog highlight the importance of analysing individuals and tissues rather than pools of individuals if flatfish metamorphosis is to be understood.

Conclusions

We report for the first time the tissues specific (skin, GI-tract and head) transcriptomes during metamorphosis of a flatfish species Hippoglossus hippoglossus with a high economic value. The study contributes substantially to the molecular resources available for this species and will be an important tool for identifying new potential molecular markers for solving problems related to Atlantic halibut production during metamorphosis. The Atlantic halibut skin transcriptome is a powerful resource for studying the asymmetric pigmentation pattern, as well as the putative cross-talk with the THs axis. Questions relating to the possible asymmetric responsiveness to the THs of both ocular and abocular (blind) sides of skin during metamorphosis remain to be resolved.

The candidate TH-responsive genes identified in the transcriptomes generated will be the subject of future studies to assess tissues responsiveness, and how it is correlated with temporal changes in elements of TH signaling and metabolism during flatfish metamorphosis. Further studies will be essential to identify the tissue specific mechanisms underlying the timing and programming of the developmental events occurring during metamorphosis. The involvement of THs in a late developmental event, metamorphosis, highlights an emerging research area: the regulatory role of hormones in early development.

Methods

Sampling

The samples of Atlantic halibut larvae for sequencing were donated by a commercial producer (Fiskeldi Eyjafjarðar Ltd., Iceland) in December 2009. Samples were collected from a standard commercial production cycle [37] undergoing normal metamorphosis by a qualified member of staff and were killed humanely. The samples for analysis were collected using established husbandry procedures and were obtained in the context of routine larval sampling protocols used by the commercial producer to verify the health, welfare and quality of the larvae. The legislation and measures implemented by the commercial producer complied with Directive 98/58/EC (protection of animals kept for farming) and production and sampling conditions were optimised to avoid unnecessary pain, suffering or injury and to maximise larval survival. The study was authorized in accordance with Portuguese legislation for the use of laboratory animals under a Group-1 license from the Direcção-Geral de Veterinária, Ministério da Agricultura, do Desenvolvimento Rural e das Pescas.

After collection the Atlantic halibut larvae (n = 6 per developmental stage) were photographed and staged using myotome height (MH) and standard length (SL) [33]. A further subdivision of stage 9 (9A, 9B, 9C) was made to account for differences in eye migration. Individual larvae were collected into RNAlater (Life Technologies, Carlsbad, USA), gently agitated for 24 h at 4 °C and then transferred to −20 °C for long term storage.

Transcriptome sequencing

Larvae (stage 5 to 9A–C, n = 6 per stage) were dissected into skin, GI-tract and head. Stage 5 larvae were divided into head and body only. Total RNA was extracted from all tissues and whole individuals (n = 5 per stage) using a Maxwell®16 System (Promega, Madison, USA) following the manufacturer’s instructions. RNA integrity and concentration was verified with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA) and only samples with RIN values equal to, or above 8 were used. The 454 (GS FLX, Roche Life Sciences, Branford, USA) and SOLiD (AB 5500xl Genetic Analyser system, Applied Biosystems, USA) sequencing was performed at the Max Planck Genome Centre (Cologne, Germany).

454 sequencing

cDNA libraries of the Atlantic halibut skin, GI-tract and head were prepared from pools of 6 samples per stage to obtain 5 μg of total RNA. Ribosomal RNA was depleted using RiboMinus™ Eukaryote Kit (Life Technologies, Carlsbad, USA) and following the manufacturer’s instructions. A cDNA Rapid Library Preparation Kit (Roche Life Sciences, USA) was used to construct sixteen cDNA libraries; head from stage 5 and head, skin and GI-tract from stages 7, 8 and 9A, 9B and 9C. Each library had a unique barcode and was amplified by emulsion PCR and sequenced on a GS-FLX platform (Roche Life Sciences, USA) following the manufacturers recommendations. The sequencing assigned quality scores are available at the NCBI Short Read Archive (SRA; Accession number: SRP044664).

SOLiD mRNA sequencing

SOLiD sequencing was carried out on cDNA libraries constructed using SOLiD™ Total RNA-Seq Kit (Applied Biosystems, CA). Ribosomal RNA was depleted from total RNA of whole individual Atlantic halibut extracted as described above. Sixteen cDNA libraries were prepared from individual larvae and included three stage 7, three stage 8, three stage 9A, three stage 9B and two stage 9C libraries from whole individual Atlantic halibut. For the non-metamorphic samples, stage 5 and juvenile, a pool of RNA from 3 individuals was used to prepare the libraries due to sequencing space constraints and cost. Each library had a unique barcode and cDNA was purified to eliminate contaminants, amplified by emulsion PCR and sequenced on a SOLiD AB 5500xl Genetic Analyser platform (Applied Biosystems, USA). The sequencing assigned quality scores are available at the NCBI SRA (Accession number: SRP073364).

Assembly and annotation

454 sequence reads

Raw sequence reads (.sff format) from the sixteen libraries were extracted and quality clipped and sequencing adapters, primers and poly-A tails were removed. Only sequences above 100 bp were retained for assembly and food source contamination was removed by screening against Artemia species available in NCBI (38,287 sequences at 04.2012) and H. hippoglossus mitochondrial RNA (27 sequences at 04.2012) was removed using BLASTn (settings: score > 100; e-value <1e-25). After quality filtering and removal of potential contaminants, approximately 70 % of the initial reads from skin (1,200,186) and head (1,556,954) and 43 % from GI-tract (888,165) were kept for assembly. As the number of reads obtained from the stage-specific tissue libraries was highly variable and did not yield robust stage specific comparisons they were combined to produce 3 tissue-specific assemblies. The filtered reads were assembled into contigs using MIRA V3 (http://sourceforge.net/projects/mira-assembler/files/) with the command: mira project = xyx –job = denovo, est, accurate, 454 –DI: trt =/dev/slim, where xyx represent the file extracted from the *sff files of each library [123]. All singletons were discarded. Files containing the reads have been submitted to the National Center for Biotechnology Information Short Read Archive (Accession number: SRP044664; [124]). Validation of library assemblies was performed by Blastx sequence similarity searching (Altschul et al. 1997) against all the available ESTs for Atlantic halibut (21,018 sequences at 07.2014) and by manual annotation of 1 % of randomly chosen contigs from each individual assembly using BLASTx against the Reference Protein database (refseq_protein; NCBI) for vertebrate species.

Annotation of putative function was performed using a multi-step stringent local blast approach adapted from Yúfera et al. [85] (Additional file 1). Sequences were first compared against NCBI Danio rerio non-redundant protein database (db) (Blastx; e-value < 1e-20), then against Swissprot (Blastx; e-value < 1e-10) and finally against the non-redundant bony fish protein database (Blastx; e-value < 1e-10). All remaining contigs were then subject to a final Blast search against all the available ESTs for bonyfish (1,191,154 sequences at 2012).

Functional annotation of 454 tissue transcriptomes

Functional annotation of the assembled 454 tissue transcriptomes was performed using the Blast2GO program v.2.6.0 [125, 126]. Annotated sequences were mapped to gene ontology (GO) terms using the following settings: annotation cut-off: 55; minimum GO weight: 5; and e-value: 1e-06. To augment the GO annotation, ANNEX analyses that gives manually curated univocal relationships between GO terms were also used (https://www.blast2go.com/; [127]). Enrichment of GO terms between different tissues was established using a Fisher’s exact test and applying a False Discovery Rate (FDR) adjusted p-value of 0.05. Metabolic KEGG pathway analysis was performed based on the Enzyme Code (EC) obtained for each GO term during the functional annotation step. Each EC was mapped to the corresponding metabolic pathway. The unique and specific tissue GO terms resulting from the 454 tissue transcriptome was analysed with REVIGO (http://revigo.irb.hr/, [128]), which uses a simple clustering algorithm that relies on semantic similarity measures to find representative subsets of GO terms.

SOLiD sequence reads

A backbone assembly was created using Newbler and CAP3 (with default parameters) combining all the 454 reads from the different tissues and stages together with the 21,018 ESTs available in public databases (July 2014) to produce 37,073 contigs. The contigs are available at http://ramadda.nerc-bas.ac.uk/repository in the folder: NERC-BAS datasets/Genomics/Transcriptomes/Hippoglossus_hippoglossus. This was then used as the reference for mapping the SOLiD sequences. The paired reads obtained for halibut samples were as follows: 9,532,993 for the pooled stage 5 sample; 11,277,613, 10,533,234 and 9,765,750 for the three stage 7 larvae; 9,658,624, 6,121,383 and 9,518,793 for the 3 stage 8 larvae; 107,948,772, 6,509,721 and 6,343,735 for the 3 stage 9A larvae; 11,101,201, 10,357,984 and 5,282,230 for the three stage 9B larvae; 12,950,027 and 9,195,626 for the two stage 9C larvae and 11,744,701 for the pooled juvenile sample. The reads were mapped to the Atlantic halibut contigs with Maq [129], using default parameters.

Expression analysis was carried out using pairwise comparisons rather than a factorial design. SOLiD sequences that failed to map to contigs were excluded from the analysis. Normalisation was carried out by dividing counts by library size. Differential expression was established using two approaches to increase stringency: a two-fold expression level difference, and the use of a linear model in Bayseq [130], with a Benjamini-Hochberg adjustment for multiple testing [131] with a cut-off set at 0.05. For the linear model, a proxy replication for mapping variance consisted of the separate mappings of the paired end reads to the contigs. Only mappings in which both paired end reads mapped to the same contig were used to generate expression levels and calculate significance of expression. The probable identity of the genes to which SOLiD reads mapped was established by sequence similarity searches using Blast [132], http://blast.ncbi.nlm.nih.gov/Blast.cgi) against the GenBank non redundant database [133] with a threshold score of < 1e-10.

Identification of tissue-specific and TH-responsive genes

Candidate genes underlying the tissue specific changes (454 transcriptomes and SOLiD analysis) that accompany metamorphosis, were identified by comparing the annotated contigs against 4 in-house databases of genes (setting: e-value ≤ 1e-20), identified using the QuickGO (http://www.ebi.ac.uk/QuickGO/) database and through literature searches. Databases included: 1) database enriched with skin-specific genes, eg. pigmentation, skin development and morphogenesis (141 sequences; Additional file 16); 2) a GI-tract-specific database enriched with genes involved in development, morphogenesis and acid secretion (179 sequences; Additional file 17); 3) genes involved in TH signalling and metabolism and genes involved in thyroid gland development (62 sequences; Additional file 18); 4) TH-responsive genes (189 sequences; Additional file 19).

SOLiD whole larval transcriptome data was used to identify transcripts with differential expression in three individuals per stage by carrying out pairwise comparisons. Two strategies were used for pairwise comparisons: i) pairwise comparisons between metamorphic stages (7, 8, 9A, 9B, 9C) or ii) pairwise comparisons between metamorphic stages (7, 8, 9A, 9B, 9C) and premetamorphic stage 5 larvae or juveniles. Since pairwise comparisons between metamorphic stages (7, 8, 9A, 9B and 9C) yielded very few differential transcripts only the results of the comparison against stage 5 or juveniles was analyzed in detail.

Quantitative real-time RT-PCR (qPCR)

A set of 16 transcripts was assessed by qPCR, using a subsample of the RNA extracted from whole larvae used for SOLiD transcriptome sequencing. Correlation analysis was carried out between transcripts with differential expression in SOLiD analysis (n = 5) and their expression determined by qPCR. The selected transcripts included genes with a constant expression in all developmental stages, transcripts with modified expression during development and genes involved in the TH cascade (Additional file 20). Specific primers for the target genes were designed based on the cDNA contig sequences (Additional file 20) using Beacon Design and Primer Premier 5.0 software (Premier Biosoft Int., Palo Alto, CA). The reference gene transcripts used to normalize the cDNA used for the PCR reactions were elongation factor I alpha (EF1A1) and 40S ribosomal protein S4 (RPS4) (Additional file 20).

For cDNA synthesis total RNA (10 μg) was first treated with Turbo DNA-Free kit (Ambion, Life Technologies, Carlsbad, USA) to remove contaminating genomic DNA. cDNA synthesis was performed with 500 ng of DNAse treated total RNA, 200 ng of random hexamers (GE Healthcare, Amersham, UK), 100 U of RevertAid M-MuLV Reverse Transcriptase (Fermentas, St Leon-Rot, Germany), 8 U of Ribolock RNase inhibitor (Fermentas, St Leon-Rot, Germany), and 0.5 mM dNTP’s. qPCR reactions were performed in duplicate using SsoFast™ EvaGreen® Supermix (Bio-Rad, Marnes La Coquette, France) chemistry in a StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, USA). The qPCR cycling conditions were: 30 s at 95 °C; 45 cycles of 5 s at 95 °C and 10 s at the optimal temperature for primer pairs (Additional file 18). A final melting curve over a range of 60–95 °C was performed for all reactions. Standard curves relating initial template quantity to amplification cycle were generated from the target gene cloned in pGEM®-T Easy (Promega, Madison, USA) using a 10-fold stepwise serial dilution series (initial concentration, 10⁸ copies amplicon/μl).

The qPCR efficiency for primer pairs ranged from 85 % and 100 % with an R² ≥ 0.99 (Additional file 20). The geometric mean of the reference genes RPS4 and EFIAI was used to normalize the qPCR data. Statistical analysis of the relative gene expression between stages was analyzed by one-way, ANOVA using SigmaPlot v10.0 (Systat Software, Inc., CA, USA) after checking for homogeneity. Tukey’s post-hoc test was used for pair wise multiple comparisons. The expression of transcripts in the halibut stages analyzed is presented as the mean ± standard error of the mean (SEM). Pearson correlation analysis was used to compare the qPCR relative gene expression levels and SOLiD differential count analysis. For correlation analysis six genes from SOLiD analysis with either a constant (RPL7 and FAU) or variable expression (Gloα1, Cpa2, ApoAI, and Krt1i2) during metamorphosis were selected. Statistical significance was established at p < 0.05.

Ethics and consent to participate

All experimental procedures involving animals complied with the Directive 98/58/EC (protection of animals kept for farming) and were authorized in accordance with Portuguese legislation for the use of laboratory animals under a Group-1 license from the Direcção-Geral de Veterinária, Ministério da Agricultura, do Desenvolvimento Rural e das Pescas.

Consent to publish

Not applicable

Availability of data and material

The 454 sequences for Atlantic halibut obtained in this study are available at the NCBI SRA under the accession number: SRP044664 and the consensus sequences of the contigs are available at http://ramadda.nerc-bas.ac.uk/repository in the folder: NERC-BAS datasets/Genomics/Transcriptomes/Hippoglossus_hippoglossus. All SOLiD sequence data were submitted to the NCBI SRA with the accession number: SRP073364.

The gene acronyms and respective full names are given in Additional files 15-18.