The transcriptome of metamorphosing flatfish
BMC Genomics volume 17, Article number: 413 (2016)
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.
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.
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.
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 . 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 . 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 .
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. ). 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)  and otolith mineralization in the Southern flounder (Paralichthys lethostigma) . 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. ). 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 ) 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).
454 transcriptome sequencing
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.
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).
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).
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).
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).
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).
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.
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).
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).
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−7) 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−4) was obtained when the juvenile post-metamorphic stage was included in the correlation analysis between SOLiD and qPCR data (Fig. 6b).
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 . 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 . 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 .
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 , in the common sole a pooled larval and adult liver and GI-tract transcriptome was established  and in Senegalese sole and common sole reference transcriptomes were derived by sequencing several tissue from juveniles and adults . 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 .
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 , 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 . 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  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  and fish under diverse challenges (e.g. viral challenge) , 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  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 . 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 . 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 . 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 , 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 . 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 . 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 , and suggests that modification of the thyroid tissue is essential for successful metamorphosis  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 . 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 . 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 . 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 . Furthermore, although in general metamorphosis is comparable between fish and amphibians their divergent evolution, biology and physiology  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  and is reminiscent of what occurs in Senegalese sole . 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  and Senegalese sole . 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 . 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 . 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 (T3RE) 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.  and Morvan-Dubois et al. ). 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 . 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.
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.
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  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) . 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.
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).
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 . All singletons were discarded. Files containing the reads have been submitted to the National Center for Biotechnology Information Short Read Archive (Accession number: SRP044664; ). 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.  (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/; ). 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/, ), 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 , 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 , with a Benjamini-Hochberg adjustment for multiple testing  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 , http://blast.ncbi.nlm.nih.gov/Blast.cgi) against the GenBank non redundant database  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, 108 copies amplicon/μl).
The qPCR efficiency for primer pairs ranged from 85 % and 100 % with an R2 ≥ 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
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.
Association of Animal Behaviour
Basic local alignment search tool
Bone morphogenic protein
Differentially expressed qPCR, Quantitative real-time polymerase chain reaction
Expressed sequence tags
False discovery rate
Kyoto encyclopedia of genes and genomes
National Center for Biotechnology Information
Next-generation sequencing technology
Reduce + Visual Gene Ontology
retinoid X receptor
standard error of the mean
Short Read Archive
transforming growth factor beta
Gilbert LI, Tata JR, Atkinson BG. Metamorphosis: postembryonic reprogramming of gene expression in amphibian and insect cells. San Diego: Academic; 1996.
Schreiber AM, Specker JL. Metamorphosis in the Summer Flounder, Paralichthys dentatus: Thyroidal Status Influences Gill Mitochondria-Rich Cells. Gen Comp Endocrinol. 2000;117(2):238–50.
Nakamura M, Ohki S, Suzuki A, Sakai K. Coral Larvae under Ocean Acidification: Survival, Metabolism, and Metamorphosis. PLoS One. 2011;6(1):e14521.
Rubio M, de Horna A, Belles X. MicroRNAs in metamorphic and non-metamorphic transitions in hemimetabolan insect metamorphosis. BMC Genomics. 2012;13:386.
Wong YH, Wang H, Ravasi T, Qian P-Y. Involvement of Wnt Signaling Pathways in the Metamorphosis of the Bryozoan Bugula neritina. PLoS One. 2012;7(3):e33323.
Huang JH, Lozano J, Belles X. Broad-complex functions in postembryonic development of the cockroach Blattella germanica shed new light on the evolution of insect metamorphosis. Biochim Biophys Acta. 2013;1830(1):2178–87.
Power DM, Einarsdóttir IE, Pittman K, Sweeney GE, Hildahl J, Campinho MA, et al. The Molecular and Endocrine Basis of Flatfish Metamorphosis. Rev Fish Sci. 2008;16 Suppl 1:95–111.
Brown DD, Cai L. Amphibian metamorphosis. Dev Biol. 2007;306(1):20–33.
Brown DD, Cai L, Das B, Marsh-Armstrong N, Schreiber AM, Juste R. Thyroid hormone controls multiple independent programs required for limb development in Xenopus laevis metamorphosis. Proc Natl Acad Sci U S A. 2005;102(35):12455–8.
Schreiber AM, Mukhi S, Brown DD. Cell–cell interactions during remodeling of the intestine at metamorphosis in Xenopus laevis. Dev Biol. 2009;331(1):89–98.
Ishizuya-Oka A, Ueda S. Apoptosis and cell proliferation in the Xenopus small intestine during metamorphosis. Cell Tissue Res. 1996;286(3):467–76.
Ishizuya-Oka A, Shi YB. Regulation of adult intestinal epithelial stem cell development by thyroid hormone during Xenopus laevis metamorphosis. Dev Dyn. 2007;236(12):3358–68.
K-i S, Machiyama F, Nishino S, Watanabe Y, Kashiwagi K, Kashiwagi A, et al. Molecular features of thyroid hormone-regulated skin remodeling in Xenopus laevis during metamorphosis. Dev Growth Differ. 2009;51(4):411–27.
Power DM, Llewellyn L, Faustino M, Nowell MA, Björnsson BT, Einarsdottir IE, et al. Thyroid hormones in growth and development of fish. Comp Biochem Physiol C Toxicol Pharmacol. 2001;130(4):447–59.
McMenamin SK, Parichy DM. Chapter Five - Metamorphosis in Teleosts. In: Yun-Bo S, editor. Curr Top Dev Biol, vol. 103. Cambridge: Academic; 2013. p. 127–65.
Inui Y, Yamano K, Miwa S. The role of thyroid hormone in tissue development in metamorphosing flounder. Aquaculture. 1995;135(1–3):87–98.
Klaren PHM, Wunderink YS, Yúfera M, Mancera JM, Flik G. The thyroid gland and thyroid hormones in Senegalese sole (Solea senegalensis) during early development and metamorphosis. Gen Comp Endocrinol. 2008;155(3):686–94.
Manchado M, Infante C, Rebordinos L, Cañavate JP. Molecular characterization, gene expression and transcriptional regulation of thyroid hormone receptors in Senegalese sole. Gen Comp Endocrinol. 2009;160(2):139–47.
Okada N, Takagi Y, Tanaka M, Tagawa M. Fine structure of soft and hard tissues involved in eye migration in metamorphosing Japanese flounder (Paralichthys olivaceus). Anat Rec A Discov Mol Cell Evol Biol. 2003;273A(1):663–8.
Tagawa M, Aritaki M. Production of symmetrical flatfish by controlling the timing of thyroid hormone treatment in spotted halibut Verasper variegatus. Gen Comp Endocrinol. 2005;141(2):184–9.
Miwa S, Yamano K, Inui Y. Thyroid hormone stimulates gastric development in flounder larvae during metamorphosis. J Exp Zool. 1992;261(4):424–30.
Schreiber AM, Wang X, Tan Y, Sievers Q, Sievers B, Lee M, et al. Thyroid hormone mediates otolith growth and development during flatfish metamorphosis. Gen Comp Endocrinol. 2010;169(2):130–7.
Campinho MA, Silva N, Nowell MA, Llewellyn L, Sweeney GE, Power DM. Troponin T isoform expression is modulated during Atlantic halibut metamorphosis. BMC Dev Biol. 2007;7:71.
Yamano K, Takano-Ohmuro H, Obinata T, Inui Y. Effect of Thyroid Hormone on Developmental Transition of Myosin Light Chains during Flounder Metamorphosis. Gen Comp Endocrinol. 1994;93(3):321–6.
Yamano K, Miwa S, Obinata T, Inui Y. Thyroid hormone regulates developmental changes in muscle during flounder metamorphosis. Gen Comp Endocrinol. 1991;81(3):464–72.
Campinho MA, Galay-Burgos M, Silva N, Costa RA, Alves RN, Sweeney GE, et al. Molecular and cellular changes in skin and muscle during metamorphosis of Atlantic halibut (Hippoglossus hippoglossus) are accompanied by changes in deiodinases expression. Cell Tissue Res. 2012;350(2):333–46.
Campinho MA, Silva N, Sweeney GE, Power DM. Molecular, cellular and histological changes in skin from a larval to an adult phenotype during bony fish metamorphosis. Cell Tissue Res. 2007;327(2):267–84.
Gomes A, Kamisaka Y, Harboe T, Power D, Ronnestad I. Functional modifications associated with gastrointestinal tract organogenesis during metamorphosis in Atlantic halibut (Hippoglossus hippoglossus). BMC Dev Biol. 2014;14(1):11.
Zambonino J-L, Gisbert E, Sarasquete C, Navarro I, Gutiérrez J, Cahu C. Ontogeny and physiology of the digestive system of marine fish larvae. In: Cyrino J, Bureau D, Kapoor B, editors. Feeding and Digestive Functions of Fishes. Enfield: Oxford & IBH Publishing Co. Pvt. Ltd; 2008. p. 281–348.
Crane HM, Pickford DB, Hutchinson TH, Brown JA. The Effects of Methimazole on Development of the Fathead Minnow, Pimephales promelas, from Embryo to Adult. Toxicol Sci. 2006;93(2):278–85.
Coady K, Marino T, Thomas J, Currie R, Hancock G, Crofoot J, et al. Evaluation of the amphibian metamorphosis assay: exposure to the goitrogen methimazole and the endogenous thyroid hormone L-thyroxine. Environ Toxicol Chem. 2010;29(4):869–80.
Liu YW, Chan WK. Thyroid hormones are important for embryonic to larval transitory phase in zebrafish. Differentiation. 2002;70(1):36–45.
Sæle Ø, Solbakken JS, Watanabe K, Hamre K, Power D, Pittman K. Staging of Atlantic halibut (Hippoglossus hippoglossus L.) from first feeding through metamorphosis, including cranial ossification independent of eye migration. Aquaculture. 2004;239(1-4):445–65.
Rosenkilde P, Ussing AP. What mechanisms control neoteny and regulate induced metamorphosis in urodeles? Int J Dev Biol. 1996;40(4):665–73.
Vlaeminck-Guillem V, Safi R, Guillem P, Leteurtre E, Duterque-Coquillaud M, Laudet V. Thyroid hormone receptor expression in the obligatory paedomorphic salamander Necturus maculosus. Int J Dev Biol. 2006;50(6):553–60.
Callery EM, Elinson RP. Thyroid hormone-dependent metamorphosis in a direct developing frog. Proc Natl Acad Sci U S A. 2000;97(6):2615–20.
Galay-Burgos M, Power DM, Llewellyn L, Sweeney GE. Thyroid hormone receptor expression during metamorphosis of Atlantic halibut (Hippoglossus hippoglossus). Mol Cell Endocrinol. 2008;281(1–2):56–63.
Pereiro P, Balseiro P, Romero A, Dios S, Forn-Cuni G, Fuste B, et al. High-Throughput Sequence Analysis of Turbot (Scophthalmus maximus) Transcriptome Using 454-Pyrosequencing for the Discovery of Antiviral Immune Genes. PLoS One. 2012;7(5):e35369.
Ferraresso S, Bonaldo A, Parma L, Cinotti S, Massi P, Bargelloni L, et al. Exploring the larval transcriptome of the common sole (Solea solea L.). BMC Genomics. 2013;14:315.
Benzekri H, Armesto P, Cousin X, Rovira M, Crespo D, Merlo M, et al. De novo assembly, characterization and functional annotation of Senegalese sole (Solea senegalensis) and common sole (Solea solea) transcriptomes: integration in a database and design of a microarray. BMC Genomics. 2014;15(1):952.
Douglas SE, Knickle LC, Kimball J, Reith ME. Comprehensive EST analysis of Atlantic halibut (Hippoglossus hippoglossus), a commercially relevant aquaculture species. BMC Genomics. 2007;8:144.
Mommens M, Fernandes J, Bizuayehu T, Bolla S, Johnston I, Babiak I. Maternal gene expression in Atlantic halibut (Hippoglossus hippoglossus L.) and its relation to egg quality. BMC Research Notes. 2010;3(1):1–11.
Yúfera M, Moyano FJ, Astola A, Pousão-Ferreira P, Martínez-Rodríguez G. Acidic Digestion in a Teleost: Postprandial and Circadian Pattern of Gastric pH, Pepsin Activity, and Pepsinogen and Proton Pump mRNAs Expression. PLoS One. 2012;7(3):e33687.
Kristiansson E, Asker N, Forlin L, Larsson DG. Characterization of the Zoarces viviparus liver transcriptome using massively parallel pyrosequencing. BMC Genomics. 2009;10:345.
Coppe A, Pujolar JM, Maes GE, Larsen PF, Hansen MM, Bernatchez L, et al. Sequencing, de novo annotation and analysis of the first Anguilla anguilla transcriptome: EeelBase opens new perspectives for the study of the critically endangered European eel. BMC Genomics. 2010;11:635.
Fu B, He S. Transcriptome Analysis of Silver Carp (Hypophthalmichthys molitrix) by Paired-End RNA Sequencing. DNA Res. 2012;19(2):131–42.
Galay-Burgos M, Llewellyn L, Bjornsson BT, Pittman K, Power DM, Smaradottir H, et al. Isolation of subtractive clones enriched in sequences that show differential expression during metamorphosis of Atlantic halibut (Hippoglossus hippoglossus). GenBank. 2006 [http://www.ncbi.nlm.nih.gov/nucest/?term=LIBEST_020648].
Cerda J, Mercade J, Lozano JJ, Manchado M, Tingaud-Sequeira A, Astola A, et al. Genomic resources for a commercial flatfish, the Senegalese sole (Solea senegalensis): EST sequencing, oligo microarray design, and development of the Soleamold bioinformatic platform. BMC Genomics. 2008;9:508.
Salem M, Rexroad 3rd CE, Wang J, Thorgaard GH, Yao J. Characterization of the rainbow trout transcriptome using Sanger and 454-pyrosequencing approaches. BMC Genomics. 2010;11:564.
Zhang Z, Wang Y, Wang S, Liu J, Warren W, Mitreva M, et al. Transcriptome Analysis of Female and Male Xiphophorus maculatus Jp 163 A. PLoS One. 2011;6(4):e18379.
Zouboulis CC. Human skin: an independent peripheral endocrine organ. Horm Res. 2000;54(5-6):230–42.
Roosterman D, Goerge T, Schneider SW, Bunnett NW, Steinhoff M. Neuronal control of skin function: the skin as a neuroimmunoendocrine organ. Physiol Rev. 2006;86(4):1309–79.
Shoulders MD, Raines RT. Collagen structure and stability. Annu Rev Biochem. 2009;78:929–58.
Le Guellec D, Morvan-Dubois G, Sire JY. Skin development in bony fish with particular emphasis on collagen deposition in the dermis of the zebrafish (Danio rerio). Int J Dev Biol. 2004;48(2-3):217–31.
Olaso E, Lin HC, Wang LH, Friedman SL. Impaired dermal wound healing in discoidin domain receptor 2-deficient mice associated with defective extracellular matrix remodeling. Fibrogenesis & tissue repair. 2011;4(1):5.
Vogel WF, Abdulhussein R, Ford CE. Sensing extracellular matrix: An update on discoidin domain receptor function. Cell Signal. 2006;18(8):1108–16.
Page RB, Voss SR, Samuels AK, Smith JJ, Putta S, Beachy CK. Effect of thyroid hormone concentration on the transcriptional response underlying induced metamorphosis in the Mexican axolotl (Ambystoma). BMC Genomics. 2008;9:78.
Page RB, Monaghan JR, Walker JA, Voss SR. A model of transcriptional and morphological changes during thyroid hormone-induced metamorphosis of the axolotl. Gen Comp Endocrinol. 2009;162(2):219–32.
Ishizuya-Oka A. Amphibian organ remodeling during metamorphosis: insight into thyroid hormone-induced apoptosis. Dev Growth Differ. 2011;53(2):202–12.
Brunelli E, Bernabò I, Coscarelli F, La Russa D, Tripepi S. Remodelling of the skin during metamorphosis in the Italian newt (Lissotriton italicus) (Amphibia, Urodela): localization pattern of keratins, stromelysin-3 (MMP-11), and pan-cadherin. Zoomorphology. 2015;134(1):135–47.
Watt FM, Fujiwara H. Cell-extracellular matrix interactions in normal and diseased skin. Cold Spring Harb Perspect Biol. 2011;3(4):a005124.
Darias MJ, Andree KB, Boglino A, Fernández I, Estévez A, Gisbert E. Coordinated Regulation of Chromatophore Differentiation and Melanogenesis during the Ontogeny of Skin Pigmentation of Solea senegalensis (Kaup, 1858). PLoS One. 2013;8(5):e63005.
Ottesen OH, Strand HK. Growth, development, and skin abnormalities of halibut (Hippoglossus hippoglossus L) juveniles kept on different bottom substrates. Aquaculture. 1996;146(1-2):17–25.
Luizi FS, Gara B, Shields RJ, Bromage NR. Further description of the development of the digestive organs in Atlantic halibut (Hippoglossus hippoglossus) larvae, with notes on differential absorption of copepod and Artemia prey. Aquaculture. 1999;176(1–2):101–16.
Gomes AS, Alves RN, Rønnestad I, Power DM. Orchestrating change: The thyroid hormones and GI-tract development in flatfish metamorphosis. Gen Comp Endocrinol. 2015;220:2–12.
Zambonino Infante JL, Cahu CL. High dietary lipid levels enhance digestive tract maturation and improve dicentrarchus labrax larval development. J Nutr. 1999;129(6):1195–200.
Ozkizilcik S, Chu F-LE, Place AR. Ontogenetic changes of lipolytic enzymes in striped bass (Morone saxatilis). Comp Biochem Physiol B Biochem Mol Biol. 1996;113(3):631–7.
Heimeier R, Das B, Buchholz D, Fiorentino M, Shi Y-B. Studies on Xenopus laevis intestine reveal biological pathways underlying vertebrate gut adaptation from embryo to adult. Genome Biol. 2010;11(5):R55.
Douglas SE, Knickle LC, Williams J, Flight RM, Reith ME. A first generation Atlantic halibut Hippoglossus hippoglossus (L.) microarray: application to developmental studies. J Fish Biol. 2008;72(9):2391–406.
Wang X, Tan Y, Sievers Q, Sievers B, Lee M, Burrall K, et al. Thyroid hormone-responsive genes mediate otolith growth and development during flatfish metamorphosis. Comp Biochem Physiol A Mol Integr Physiol. 2011;158(1):163–8.
Fernandez LP, Lopez-Marquez A, Santisteban P. Thyroid transcription factors in development, differentiation and disease. Nat Rev Endocrinol. 2015;11(1):29–42.
Logan CY, Nusse R. The Wnt signaling pathway in development and disease. Annu Rev Cell Dev Biol. 2004;20:781–810.
Botchkarev VA, Sharov AA. BMP signaling in the control of skin development and hair follicle growth. Differentiation. 2004;72(9-10):512–26.
Ingham PW, McMahon AP. Hedgehog signaling in animal development: paradigms and principles. Genes Dev. 2001;15(23):3059–87.
Janicke M, Carney TJ, Hammerschmidt M. Foxi3 transcription factors and Notch signaling control the formation of skin ionocytes from epidermal precursors of the zebrafish embryo. Dev Biol. 2007;307(2):258–71.
Paridaen JT, Huttner WB. Neurogenesis during development of the vertebrate central nervous system. EMBO Rep. 2014;15(4):351–64.
Pascual A, Aranda A. Thyroid hormone receptors, cell growth and differentiation. Biochim Biophys Acta. 2013;1830(7):3908–16.
Huggins P, Johnson CK, Schoergendorfer A, Putta S, Bathke AC, Stromberg AJ, et al. Identification of differentially expressed thyroid hormone responsive genes from the brain of the Mexican Axolotl (Ambystoma mexicanum). Comp Biochem Physiol C Toxicol Pharmacol. 2012;155(1):128–35.
Ishizuya-Oka A, Hasebe T, Buchholz DR, Kajita M, Fu L, Shi YB. Origin of the adult intestinal stem cells induced by thyroid hormone in Xenopus laevis. FASEB J. 2009;23(8):2568–75.
Sun G, Heimeier RA, Fu L, Hasebe T, Das B, Ishizuya-Oka A, et al. Expression Profiling of Intestinal Tissues Implicates Tissue-Specific Genes and Pathways Essential for Thyroid Hormone-Induced Adult Stem Cell Development. Endocrinology. 2013;154(11):4396–407.
Veldhoen N, Crump D, Werry K, Helbing CC. Distinctive gene profiles occur at key points during natural metamorphosis in the Xenopus laevis tadpole tail. Dev Dynam. 2002;225(4):457–68.
Buchholz DR, Paul BD, Fu L, Shi Y-B. Molecular and developmental analyses of thyroid hormone receptor function in Xenopus laevis, the African clawed frog. Gen Comp Endocrinol. 2006;145(1):1–19.
Furlow JD, Kanamori A. The transcription factor basic transcription element-binding protein 1 is a direct thyroid hormone response gene in the frog Xenopus laevis. Endocrinology. 2002;143(9):3295–305.
Das B, Heimeier RA, Buchholz DR, Shi YB. Identification of direct thyroid hormone response genes reveals the earliest gene regulation programs during frog metamorphosis. J Biol Chem. 2009;284(49):34167–78.
Yúfera M, Halm S, Beltran S, Fuste B, Planas JV, Martinez-Rodriguez G. Transcriptomic characterization of the larval stage in gilthead seabream (Sparus aurata) by 454 pyrosequencing. Mar Biotechnol. 2012;14(4):423–35.
Helbing CC, Werry K, Crump D, Domanski D, Veldhoen N, Bailey CM. Expression Profiles of Novel Thyroid Hormone-Responsive Genes and Proteins in the Tail of Xenopus laevis Tadpoles Undergoing Precocious Metamorphosis. Mol Endocrinol. 2003;17(7):1395–409.
Das B, Cai L, Carter MG, Piao YL, Sharov AA, Ko MS, et al. Gene expression changes at metamorphosis induced by thyroid hormone in Xenopus laevis tadpoles. Dev Biol. 2006;291(2):342–55.
Buchholz DR, Heimeier RA, Das B, Washington T, Shi Y-B. Pairing morphology with gene expression in thyroid hormone-induced intestinal remodeling and identification of a core set of TH-induced genes across tadpole tissues. Dev Biol. 2007;303(2):576–90.
Brekken RA, Sage EH. SPARC, a matricellular protein: at the crossroads of cell-matrix communication. Matrix Biol. 2001;19(8):816–27.
Tremble PM, Lane TF, Sage EH, Werb Z. SPARC, a secreted protein associated with morphogenesis and tissue remodeling, induces expression of metalloproteinases in fibroblasts through a novel extracellular matrix-dependent pathway. J Cell Biol. 1993;121(6):1433–44.
Yan Q, Blake D, Clark JI, Sage EH. Expression of the matricellular protein SPARC in murine lens: SPARC is necessary for the structural integrity of the capsular basement membrane. J Histochem Cytochem. 2003;51(4):503–11.
Torres-Nunez E, Suarez-Bregua P, Cal L, Cal R, Cerda-Reverter JM, Rotllant J. Molecular cloning and characterization of the matricellular protein Sparc/osteonectin in flatfish, Scophthalmus maximus, and its developmental stage-dependent transcriptional regulation during metamorphosis. Gene. 2015;568(2):129–39.
Matsumoto R, Sugimoto M. Dermal matrix proteins initiate re-epithelialization but are not sufficient for coordinated epidermal outgrowth in a new fish skin culture model. Cell Tissue Res. 2007;327(2):249–65.
Plow EF, Haas TA, Zhang L, Loftus J, Smith JW. Ligand binding to integrins. J Biol Chem. 2000;275(29):21785–8.
Marchand O, Duffraisse M, Triqueneaux G, Safi R, Laudet V. Molecular cloning and developmental expression patterns of thyroid hormone receptors and T3 target genes in the turbot (Scophtalmus maximus) during post-embryonic development. Gen Comp Endocrinol. 2004;135(3):345–57.
Miwa S, Tagawa M, Inui Y, Hirano T. Thyroxine surge in metamorphosing flounder larvae. Gen Comp Endocrinol. 1988;70(1):158–63.
Yamano K. The Role of Thyroid Hormone in Fish Development with Reference to Aquaculture. JARQ. 2005;39:161–8.
Schreiber AM, Specker JL. Metamorphosis in the summer flounder (Paralichthys dentatus): stage-specific developmental response to altered thyroid status. Gen Comp Endocrinol. 1998;111(2):156–66.
Einarsdóttir I, Silva N, Power D, Smáradóttir H, Björnsson B. Thyroid and pituitary gland development from hatching through metamorphosis of a teleost flatfish, the Atlantic halibut. Anat Embryol. 2006;211(1):47–60.
Manchado M, Infante C, Asensio E, Planas JV, Cañavate JP. Thyroid hormones down-regulate thyrotropin β subunit and thyroglobulin during metamorphosis in the flatfish Senegalese sole (Solea senegalensis Kaup). Gen Comp Endocrinol. 2008;155(2):447–55.
Cai L, Brown DD. Expression of type II iodothyronine deiodinase marks the time that a tissue responds to thyroid hormone-induced metamorphosis in Xenopus laevis. Dev Biol. 2004;266(1):87–95.
Campinho MA, Galay-Burgos M, Sweeney GE, Power DM. Coordination of deiodinase and thyroid hormone receptor expression during the larval to juvenile transition in sea bream (Sparus aurata, Linnaeus). Gen Comp Endocrinol. 2010;165(2):181–94.
Darras VM, Houbrechts AM, Van Herck SLJ. Intracellular thyroid hormone metabolism as a local regulator of nuclear thyroid hormone receptor-mediated impact on vertebrate development. Biochim Biophys Acta. 2015;1849(2):130–41.
Nakajima K, Fujimoto K, Yaoita Y. Programmed cell death during amphibian metamorphosis. Semin Cell Dev Biol. 2005;16(2):271–80.
Van der Geyten S, Van den Eynde I, Segers IB, Kühn ER, Darras VM. Differential expression of iodothyronine deiodinases in chicken tissues during the last week of embryonic development. Gen Comp Endocrinol. 2002;128(1):65–73.
Isorna E, Obregon MJ, Calvo RM, Vazquez R, Pendon C, Falcon J, et al. Iodothyronine deiodinases and thyroid hormone receptors regulation during flatfish (Solea senegalensis) metamorphosis. J Exp Zool B Mol Dev Evol. 2009;312b(3):231–46.
Itoh K, Watanabe K, Wu X, Suzuki T. Three Members of the Iodothyronine Deiodinase Family, dio1, dio2 and dio3, are Expressed in Spatially and Temporally Specific Patterns During Metamorphosis of the Flounder, Paralichthys olivaceus. Zool Sci. 2010;27(7):574–80.
Friesema EC, Jansen J, Milici C, Visser TJ. Thyroid hormone transporters. Vitam Horm. 2005;70:137–67.
Visser WE, Friesema EC, Visser TJ. Minireview: thyroid hormone transporters: the knowns and the unknowns. Mol Endocrinol. 2011;25(1):1–14.
Schweizer U, Köhrle J. Function of thyroid hormone transporters in the central nervous system. Biochim Biophys Acta. 2013;1830(7):3965–73.
Arjona FJ, de Vrieze E, Visser TJ, Flik G, Klaren PHM. Identification and Functional Characterization of Zebrafish Solute Carrier Slc16a2 (Mct8) as a Thyroid Hormone Membrane Transporter. Endocrinology. 2011;152(12):5065–73.
Campinho MA, Saraiva J, Florindo C, Power DM. Maternal Thyroid Hormones Are Essential for Neural Development in Zebrafish. Mol Endocrinol. 2014;28(7):1136–49.
Vatine GD, Zada D, Lerer-Goldshtein T, Tovin A, Malkinson G, Yaniv K, et al. Zebrafish as a model for monocarboxyl transporter 8-deficiency. J Biol Chem. 2013;288(1):169–80.
Connors KA, Korte JJ, Anderson GW, Degitz SJ. Characterization of thyroid hormone transporter expression during tissue-specific metamorphic events in Xenopus tropicalis. Gen Comp Endocrinol. 2010;168(1):149–59.
Geysens S, Ferran JL, Van Herck SLJ, Tylzanowski P, Puelles L, Darras VM. Dynamic mRNA distribution pattern of thyroid hormone transporters and deiodinases during early embryonic chicken brain development. Neuroscience. 2012;221:69–85.
Taylor PM, Ritchie JWA. Tissue uptake of thyroid hormone by amino acid transporters. Best Pract Res Clin Endocrinol Metab. 2007;21(2):237–51.
Van Herck SLJ, Geysens S, Delbaere J, Darras VM. Regulators of thyroid hormone availability and action in embryonic chicken brain development. Gen Comp Endocrinol. 2013;190:96–104.
Van Herck SLJ, Delbaere J, Bourgeois NMA, McAllan BM, Richardson SJ, Darras VM. Expression of thyroid hormone transporters and deiodinases at the brain barriers in the embryonic chicken: Insights into the regulation of thyroid hormone availability during neurodevelopment. Gen Comp Endocrinol. 2015;214:30–9.
Yamano K, Miwa S. Differential gene expression of thyroid hormone receptor alpha and beta in fish development. Gen Comp Endocrinol. 1998;109(1):75–85.
Buchholz DR. More similar than you think: Frog metamorphosis as a model of human perinatal endocrinology. Dev Biol. 2015;408(2):188–95.
Grimaldi AG, Buisine N, Bilesimo P, Sachs LM. High-throughput sequencing will metamorphose the analysis of thyroid hormone receptor function during amphibian development. Curr Top Dev Biol. 2013;103:277–303.
Morvan-Dubois G, Demeneix BA, Sachs LM. Xenopus laevis as a model for studying thyroid hormone signalling: From development to metamorphosis. Mol Cell Endocrinol. 2008;293(1–2):71–9.
Chevreux B, Wetter T, Suhai S. Genome sequence assembly using trace signals and additional sequence information. Journal of Computer Science & Systems Biology. 1999;99:45–56.
Gomes AS, Alves RN, Stueber K, Thorne MAS, Smáradóttir H, Reinhard R, et al. Transcriptome of the Atlantic halibut (Hippoglossus hippoglossus). Marine Genomics. 2014;18, Part B(0):101–3.
Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.
Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36(10):3420–35.
Myhre S, Tveit H, Mollestad T, Laegreid A. Additional gene ontology structure for improved biological reasoning. Bioinformatics. 2006;22(16):2020–7.
Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLoS One. 2011;6(7):e21800.
Li H, Ruan J, Durbin R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 2008;18(11):1851–8.
Hardcastle T, Kelly K. baySeq: Empirical Bayesian methods for identifying differential expression in sequence count data. BMC Bioinformatics. 2010;11(1):422.
Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Ser B Stat Methodol. 1995;57(1):289–300.
Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389–402.
Benson DA, Cavanaugh M, Clark K, Karsch-Mizrachi I, Lipman DJ, Ostell J, et al. GenBank. Nucleic Acids Res. 2013;41(Database issue):D36–42.
Searcy BT, Beckstrom-Sternberg SM, Beckstrom-Sternberg JS, Stafford P, Schwendiman AL, Soto-Pena J, et al. Thyroid hormone-dependent development in Xenopus laevis: a sensitive screen of thyroid hormone signaling disruption by municipal wastewater treatment plant effluent. Gen Comp Endocrinol. 2012;176(3):481–92.
The authors thank F. Zimmermann for help with R scripts and comments.
This research study was funded by the European Community FP7 project LIFECYCLE (FP7 222719, http://www.lifecycle.gu.se/, 17.08.2012). Ricardo N. Alves was funded by FCT (SFRH/BD/69209/2010). MSC and MAST were funded under the British Antarctic Survey Polar Sciences for Planet Earth programme.
DMP conceived and planned the project. HS developed and ran the culture system, RNA, ASG, MT performed the practical work including molecular biology. RNA, ASG, MAST, KS and RR performed the bioinformatics analyses. DMP, RNA, ASG, IR and MSC analyzed and interpreted the data. DMP, RNA, ASG drafted the manuscript. All authors revised it critically for important intellectual content. All authors have given their final approval of the version to be published.
The authors declare that they have no competing interests.
Scheme of the data processing pipeline for de novo transcriptome assembly, annotation and Gene Ontology analysis of Atlantic halibut skin, GI-tract and head transcriptomes. (PDF 263 kb)
Annotation of the skin transcriptome. Contig ID, transcript name, number of GOs and their description and enzymatic codes are shown for each contig. (XLSX 3090 kb)
Annotation of the GI-tract transcriptome. Contig ID, transcript name, number of GOs and their description and enzymatic codes are shown for each contig. (XLSX 1784 kb)
Annotation of the head transcriptome. Contig ID, transcript name, number of GOs and their description and enzymatic codes are shown for each contig. (XLSX 4280 mb)
Schematic representation of the functional annotation obtained after analysis of the transcriptomes for skin, GI-tract and head. The GO terms (level 2) used for classification were biological process (BP), molecular function (MF) and cellular component (CC). (PDF 107 kb)
Significantly overrepresented Biological Process GO terms identified for the skin transcriptome (FDR < 0.05). (DOC 89 kb)
Significantly overrepresented Biological Process GO terms identified for the GI-tract transcriptome (FDR < 0.05). (DOC 71 kb)
Significantly overrepresented Biological Process GO terms identified for the head transcriptome (FDR < 0.05). (DOC 116 kb)
List of the most representative metabolic pathways in the skin, GI-tract and head transcriptomes using KEGG analysis. (DOCX 16 kb)
The list of the top most significantly up-regulated genes between premetamorphic stage 5 and the Juvenile stage determined using a linear model in Bayseq with a Benjamin-Hochberg adjustment for multiple testing analysis with a cut-off set at 0.05 (FDR < 0.05). Contig name, Gene name, Accession number (no.), Organism and E-value are shown for each gene. (DOCX 75 kb)
Expression profile of putative TH-responsive transcripts that had a differential expression between Atlantic halibut metamorphic stages (7–9C) and the juvenile. Expression levels are represented as log2 of fold-change (juvenile/metamorphic stages expression levels). (DOCX 41 kb)
Heat map with the expression profile (log2 of fold-change) of putative thyroid hormones (TH) responsive transcripts with differential expression between juvenile and metamorphic stages of Atlantic halibut. (PDF 312 kb)
Reactome pathway analysis for the 145 differential expressed TH-responsive transcripts identified when Atlantic halibut metamorphic stages were compared with the juvenile. Reactome analysis was performed using INTREPROSCAN accession numbers obtained from the functional annotation of the Atlantic halibut transcriptome with Blast2GO. (DOCX 18 kb)
Relative gene expression analysis of transcripts involved in the TH cascade during Atlantic halibut development (stage 5 to juvenile; n = 5 per stage) using quantitative RT-PCR (qPCR). Thyroglobulin (Tg), TH receptor alpha A (TRαA), TH receptor alpha B (TRαB), TH receptor beta (TRβ), deiodinase 3 (DIO3), deiodinase 2 (DIO2), deiodinase 1 (DIO1), monocarboxylate transporter 8 (MCT8), and monocarboxylate transporter 10 (MCT10) gene expression. Results are presented as mean ± SEM of the candidate gene expression, normalized using the geometric mean of the reference genes RPS4 and EF1A1. Significant difference (p < 0.05; one-way ANOVA) of normalized transcript expression between stages are indicated by different letters. (PDF 200 kb)
Quantitative RT-PCR (qPCR) of the relative expression of ribosomal protein L7 (RPL7), 40S ribosomal protein S30 (FAU); alpha-globin 1 (Gloα1), carboxypeptidase A2 (Cpa2), apolipoprotein A-I (ApoAI) and type I keratin isoform 2 (Krt1i2). Analysis of the indicated transcripts was performed in whole Atlantic halibut larvae during development (stage 5 to juvenile; n = 5). The results are presented as mean ± SEM of the normalized expression, using the geometric mean of the reference genes RPS4 and EFIAI. Different letters represent significantly different mean values (p < 0.05; one-way ANOVA). (PDF 188 kb)
In-house skin-specific database enriched with candidate genes involved in vertebrate skin development, morphogenesis and pigmentation. Protein name, Symbol, Accession number (no.), Organism and Biological role are shown. (XLSX 19 kb)
In house GI-tract-specific database enriched with genes involved in GI-tract development, morphogenesis and acid secretion. Protein name, Symbol, Accession number (no.), Organism and Biological role (when associated) are shown. (XLSX 24 kb)
In-house database of candidate genes involved in thyroid gland development/thyroid hormone (TH) metabolism and signaling, including transcripts with a relevant role in TH synthesis, transport and activity. Protein name, Symbol, Accession number (no.), Organism and Biological role are shown. (XLSX 13 kb)
Specific primers used for qPCR gene expression analysis. Gene symbol, name and function are shown. The annealing temperature (Ta °C), amplicon length (bp), R2 and qPCR efficiency (%) are indicated for each primer pair. (DOCX 17 kb)
About this article
Cite this article
Alves, R.N., Gomes, A.S., Stueber, K. et al. The transcriptome of metamorphosing flatfish. BMC Genomics 17, 413 (2016). https://doi.org/10.1186/s12864-016-2699-x