Brains in Metamorphosis: Temporal Transcriptome Dynamics in Hatchery-Reared Flatfishes

biology Article Brains in Metamorphosis: Temporal Transcriptome Dynamics in Hatchery-Reared Flatfishes Laura Guerrero-Peña 1, Paula Suarez-Bregua 1,*, Luis Méndez-Martínez 1, Pablo García-Fernández 2, Ricardo Tur 2, Juan A. Rubiolo 3,4, Juan J. Tena 5and Josep Rotllant 1,*   Citation: Guerrero-Peña, L.; Suarez-Bregua, P.; Méndez-Martínez, L.; García-Fernández, P.; Tur, R.; Rubiolo, J.A.; Tena, J.J.; Rotllant, J. Brains in Metamorphosis: Temporal Transcriptome Dynamics in Hatchery-Reared Flatfishes. Biology 2021,10, 1256. https://doi.org/ 10.3390/biology10121256 Academic Editor: Arne Ludwig Received: 28 October 2021 Accepted: 29 November 2021 Published: 2 December 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Aquatic Biotechnology Lab., Institute of Marine Research, Spanish National Research Council (IIM-CSIC), 36208 Vigo, Spain; lguerrer[email protected] (L.G.-P.); [email protected] (L.M.-M.) 2Nueva Pescanova Biomarine Center, S.L., 36980 O Grove, Spain; [email protected] (P.G.-F.); [email protected] (R.T.) 3Facultad de Ciencias Bioquímicas y Farmacéuticas-Centro Científico y Tecnológico Acuario del Río Paraná, Universidad Nacional de Rosario, Rosario S2002LRK, Argentina; ja.r[email protected] 4Departamento de Genética, Facultad de Veterinaria, Universidad de Santiago de Compostela, 27002 Lugo, Spain 5 Centro Andaluz de Biología del Desarrollo (CABD), CSIC-Universidad Pablo de Olavide, 41013 Sevilla, Spain; [email protected] *Correspondence: [email protected] (P.S.-B.); [email protected] (J.R.) Simple Summary: Metamorphosis, the process by which the young forms of some animals develop into adult forms, has long kept scientists on tenterhooks. Flatfish undergo one of the most dramatic metamorphoses described in the entire animal kingdom, in which a symmetrical larva that swims upright and looks just like a typical baby fish becomes a completely flat, asymmetrical juvenile that will live its entire adult life attached to the bottom. To answer the question of how the same organism can generate two completely different body plans associated with different phases of the life cycle, we investigated the dynamics of the brain transcriptome, which is the regulatory center of specific endocrine-activated developmental processes during metamorphosis. Our results show, for the first time, a temporary immune system reorganization during flatfish metamorphic remodelling process. Therefore, characterizing and understanding all the developmental changes that take place during metamorphosis will assist in the understanding the importance of each of these processes in the normal development of an individual and therefore, facilitate the transfer of knowledge to prevent abnormal development or developmental pathologies. Abstract: Metamorphosis is a captivating process of change during which the morphology of the larva is completely reshaped to face the new challenges of adult life. In the case of fish, this process initiated in the brain has traditionally been considered to be a critical rearing point and despite the pioneering molecular work carried out in other flatfishes, the underlying molecular basis is still relatively poorly characterized. Turbot brain transcriptome of three developmental stages (pre-metamorphic, climax of metamorphosis and post-metamorphic) were analyzed to study the gene expression dynamics throughout the metamorphic process. A total of 1570 genes were differentially expressed in the three developmental stages and we found a specific pattern of gene expression at each stage. Unexpectedly, at the climax stage of metamorphosis, we found highly expressed genes related to the immune response, while the biological pathway enrichment analysis in pre-metamorphic and post-metamorphic were related to cell differentiation and oxygen carrier activity, respectively. In addition, our results confirm the importance of thyroid stimulating hormone, increasing its expression during metamorphosis. Based on our findings, we assume that immune system activation during the climax of metamorphosis stage could be related to processes of larval tissue inflammation, resorption and replacement, as occurs in other vertebrates. Keywords: turbot; metamorphosis; brain; RNA; sequencing; transcriptome Biology 2021,10, 1256. https://doi.org/10.3390/biology10121256 https://www.mdpi.com/journal/biology Biology 2021,10, 1256 2 of 16 1. Introduction Metamorphosis is a post-embryonic process that involves radical changes in morphology, physiology and habitat, leading to stage-specific organs and structures [ 1 ]. Metamorphosis is widely distributed throughout the animal kingdom [ 2 ], which results in varying types of metamorphosis that probably do not share the same evolutionary mechanisms [ 3 ]. The well-known metamorphoses are those carried out by insects [ 4 , 5 ] and amphibians [ 6 ], but many other invertebrates and vertebrates have to face this process [7]. Most teleost fish undergo metamorphosis in the transition from larva to juvenile stage [ 8 , 9 ]; however, the most dramatic metamorphosis known in fish is that of flatfish (order Pleuronectiformes) [ 10 ], which includes species of high economic value such as turbot (Scophthalmus maximus), Senegalese sole (Solea senegalensis) or Atlantic halibut (Hippoglossus hippoglossus). During flatfish metamorphosis, a pelagic and bilateral symmetric larva becomes a benthonic asymmetric juvenile. Profound internal and external changes, such as eye migration to the opposite side of the body, remodelling of the craniofacial complex [ 11 – 13 ] and redistribution of the skin pigmentation [ 14 ], lead to an asymmetric juvenile fish adapted to benthic life. Metamorphosis is an energy demanding process that has an important impact on the feeding, growth and, in some cases, a higher mortality rate in flatfishes [ 15 ]. The unsuccessful larva-to-juvenile transition can lead to malformed tissues due to incomplete migration of the eye, bone deformity and/or abnormal pigmentation deposition, ultimately affecting the commercial value of flatfishes [16]. As occurs in amphibians [ 17 ], thyroid hormones (THs) have an important role as endocrine regulators of flatfish metamorphosis. In fact, TH action has been shown as mandatory to successfully carry out this process [ 18 ]. Treatments with inhibitors of THs arrest metamorphosis, interrupting eye migration alongside development and growth rate [ 19 ]. TH production is determined by the hypothalamus-pituitary-thyroid (HPT) axis. The pituitary, located in the brain and ventral to the hypothalamus [ 20 ], prompts the thyroid gland by releasing thyroid stimulating hormone (TSH). Thyroid glands are the only component of the HPT axis that is located outside the brain and secretes THs (triiodothyronine [T3] and tetraiodothyronine [T4]) to act on target tissues [ 21 , 22 ] via thyroid hormone receptors (i.e., TH receptor alpha [TR α ] and TH receptor beta [TR β ]) [ 9 ]. Then, the functional adaptive changes associated with post-embryonic development are supposed to be tightly regulated by the expression of specific genes. Although one of the two major components of the HPT axis is found in the brain, this organ does not seem to suffer an obvious morphological remodelling during the flatfish metamorphosis. Despite the fact that the horizontal semicircular canals, with respect to the eyes, change from parallel to perpendicular in the transition from larva to juvenile, the symmetry of the brain is not affected [ 23 ]. The only regions of the brain known to become asymmetric after metamorphosis are the olfactory lobes [ 24 ] and telencephalic hemispheres [ 25 ]. Neither has evidence of asymmetric remodelling been found in the neuronal nerve pathways [10]. The present study aims to characterize the gene expression profile and molecular mechanisms involved in the flatfish brain during metamorphosis. Turbot brain transcriptome of three developmental stages (pre-metamorphic, climax of metamorphosis and post-metamorphic) were analyzed to study the gene expression dynamics throughout the metamorphic process. The brain was selected as the target tissue since it is the regulatory center of specific endocrine-activated developmental processes during metamorphosis. We focused on the differential gene expression profile throughout metamorphosis and the enriched pathways in each developmental stage. 2. Materials and Methods 2.1. Fish Collection and Sampling Newborn turbots (Scophthalmus maximus) reared under a standard commercial production cycle were supplied by the company Insuiña SL, Grupo Nueva Pescanova (Pontevedra, Biology 2021,10, 1256 3 of 16 Spain). Fish from a single-pair mating were collected at different stages by experienced company staff. The number of fish sampled for all experimental procedures was estimated according to the minimum number of animals necessary to provide reliable and robust statistical results. The metamorphic stages were defined following the morphological criteria described by Al-Maghazachi and Gibson [ 26 ] and Suarez-Bregua [ 27 ], based on eye migration, whole body symmetry and rearing temperature (18 ◦ C): pre-metamorphic stage (stage 3b: 15 days post fertilization [dpf]), before eye migration on a symmetrical larva; metamorphic climax (stage 4d: 30 dpf), larva exhibiting asymmetrical features with upper edge of right eye visible from left side; and post-metamorphic stage (stage 5c: 57 dpf), asymmetric juveniles that achieved complete eye migration. Individual fish samples at each metamorphic stage (N= 3 independent biological replicates per stage) were euthanized using a lethal dose of MS-222 (250 mg/L for 30–40 min) [ 28 ] (Sigma-Aldrich, Saint Louis, MO, USA), photographed and dissected with a Leica M165FC stereomicroscope equipped with a DFC310FX camera (Leica, Wetzlar, Germany). Ethical approval (AGL2017-89648P) for all studies was obtained from the Institutional Animal Care and Use Committee of the IIM-CSIC Institute in accordance with the National Advisory Committee for Laboratory Animal Research Guidelines licensed by the Spanish Authority (RD53/2013). This work was in conformance with the European animal directive (2010/63/UE) for the protection of experimental animals. 2.2. RNA Isolation and Sequencing Fish brains (N = 9) were dissected. Briefly, the turbot head was cut off and the top of skull was opened to remove all brain tissue, including the pituitary gland. Samples were then fixed in RNAlater (Thermo Fisher Scientific, Waltham, MA, USA) for 24 h at 4 ◦ C and stored at − 80 ◦ C until use. Brain samples were removed from RNAlater solution (Invitrogen, Waltham, MA, USA) and homogenized in RLT buffer (RNeasy Mini Kit, Qiagen, Venlo, Germany). Total RNA was extracted and purified using the RNeasy Mini Kit (Qiagen) with on-column DNase digestion (Qiagen) according to the manufacturer’s instructions. RNA concentration was quantified on a Qubit 4 fluorometer (Thermo Fisher Scientific) and RNA integrity (RIN (RNA integrity number) > 8) was verified on an Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Approximately 1 µ g of total RNA was initially used for BGISEQ-500 standard library construction at BGI (Beijing Genomics Institute, Shenzhen, China). Prepared cDNA libraries were sequenced on a BGISEQ-500 platform and single-end reads of 50 base pairs (bp) length were generated per sample. 2.3. Transcriptome Analysis and Annotation Reads were quality checked (phred score > 30) using FastQC v0.11.8 (http://www.bi oinformatics.babraham.ac.uk/projects/fastqc/; accesed 18 November 2019) and mapped to the turbot genome assembly (ASM318616v1) [ 29 ] using STAR v2.7.0e alignment software [ 30 ]. The turbot reference genome and respective annotation file were downloaded from Ensembl Genome Browser (ftp://ftp.ensembl.org/pub/release-101/fasta/scophtha lmus_maximus/; accesed 25 November 2019). HTseq v0.10.0 [ 31 ] was used to transform uniquely mapped reads into counts and assign them to genes. We performed a functional annotation of the whole turbot genome using Sma3s v2 software [ 32 ]. We first obtained the predicted amino acid sequences from Ensembl REST server through the API with a custom script, and then, query sequences were compared with non-redundant protein sequences of the Swissprot and TrEMBL vertebrate databases using an E-value threshold of 1 × 10 −6 . Subsequently, the annotated genes were assigned to the three main categories of Gene Ontology (GO): biological process (BP), molecular function (MF), and cell component (CC). The large computational operations were carried out using the resources of the Supercomputing Center of Galicia (CESGA). Biology 2021,10, 1256 4 of 16 2.4. Differential Expression Analysis, Clustering and GO Enrichment Gene count data were normalized and pairwise comparisons were performed to identify differentially expressed genes (DEGs) with the DESeq2 R package v1.26.0 [ 33 ]. p-values were adjusted (padj) by false discovery rate (FDR) [ 34 ]. Only genes with a padj < 0.01 and Log2 fold change (Log2FC) ≤ − 2 or ≥ 2 were considered as DEGs. DEGs were analyzed according to their expression pattern throughout three different approaches: (1) Hierarchical clustering via heatmap; (2) Soft clustering using Mfuzz software [ 35 ]; (3) overlapping clustered up and downregulated genes using custom Venn diagram. GO enrichment analysis was carried out using clusterProfiler R package v3.14.3 [ 36 ], based on hypergeometric distribution and FDR control [ 34 ]. We used the previously functionally annotated turbot genome as a background. 2.5. Quantitative Real-Time PCR (qRT-PCR) To validate the RNA sequencing and transcriptome analysis, a set of five selected genes were evaluated by qRT-PCR. Genes related to immune system response (chitinase 3, chit3; and Interferon-induced helicase C domain-containing protein, ifih1) and braindriven metamorphic remodeling (thyroid hormone receptor alpha, thra; thyroid-stimulating hormone subunit beta a, tshba; and ependymin, epd) were selected. For cDNA synthesis, 200 ng of the total RNA isolated from each sample was reversetranscribed according to the Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) protocol. Samples were amplified in duplicate containing 10 µ L of PowerUp SYBR Green Master Mix (2 × ) (Thermo Fisher Scientific), 1 µ L of 0.5 µ M of each primer, 7 µ L nuclease-free water, and 1 µ L of cDNA template. qRT-PCR reactions were analyzed with a QuantStudio3 Real-Time PCR System (Thermo Fisher Scientific) under the following cycling conditions: initial uracil-DNA-glycosylase step at 50 ◦ C for 2 min, Dual-Lock ™ DNA polymerase activation at 95 ◦ C for 2 min, followed by 40 cycles of denaturation at 95 ◦ C for 15 s and annealing/extension at 60 ◦ C for 1 min. Gene expression of chit3, epd,thra,tshba and ifih1 genes was assessed in two independent experiments by using the efficiency-calibrated method, as previously described [ 37 ]. Relative mRNA expression levels were normalized to the housekeeping 18S ribosomal gene. Primer sets used for each gene are listed in Supplementary Table S1. 3. Results 3.1. Transcriptome Assembly and Annotation We sequenced the brain transcriptomes in three key developmental stages across turbot metamorphosis (pre-metamorphic, climax and post-metamorphic stages) (Figure 1). Nine cDNA libraries (three replicates per each metamorphic stage) were sequenced and more than 208 million 50 bp single-end reads were generated. Reads were processed for subsequent transcriptome analysis (Table 1). The reads were mapped to the turbot reference genome, obtaining average mapping rates of 88.98%, 88.78% and 89.23% for premetamorphic, climax and post-metamorphic stages, respectively. The average number of genes assigned was 19,119, 19,045 and 18,996 for pre-metamorphic, climax and postmetamorphic stages, respectively (Table 1, Figure 2). Reads assigned to each gene of turbot genome are listed in Supplementary Table S2. Table 1. Summary of RNA sequencing, assembly and annotation data of the turbot brain samples at three key metamorphic stages. Pre-Metamorphic Stage Climax Stage Post-Metamorphic Stage Total R1 R2 R3 R1 R2 R3 R1 R2 R3 Reads 24,000,160 22,065,907 23,172,275 23,261,839 23,167,451 23,191,562 23,275,068 23,173,993 23,358,252 208,666,507 Mapped reads 21,366,876 19,668,853 20,569,271 20,699,424 20,540,205 20,569,248 20,731,453 20,676,181 20,882,737 185,704,248 Mapping rate (%) 89.03 89.14 88.77 88.98 88.66 88.69 89.07 89.22 89.40 - Number of genes 19,137 19,058 19,162 19,042 19,045 19,048 18,977 18,993 19,017 - Biology 2021,10, 1256 5 of 16 Biology 2021, 10, x FOR PEER REVIEW 5 of 18 Figure 1. Turbot developmental stages selected during the metamorphosis process for transcriptome analysis. From left to right: (a) pre-metamorphic symmetrical larva prior to eye migration (stage 3b: 15dpf) (b) asymmetrical larva at the metamorphic climax with upper edge of migrating right eye visible from left side (stage 4d: 30dpf) (c) post-metamorphic asymmetrical juvenile with upper eye entirely placed on the left side (stage 5c: 57dpf). Fish were reared at 18 °C. Scale bars (a) 2mm, (b,c) 5 mm. Table 1. Summary of RNA sequencing, assembly and annotation data of the turbot brain samples at three key metamorphic stages. Pre-Metamorphic Stage Climax Stage Post-Metamorphic Stage Total R1 R2 R3 R1 R2 R3 R1 R2 R3 Reads 24,000,160 22,065,907 23,172,275 23,261,839 23,167,451 23,191,562 23,275,068 23,173,993 23,358,252 208,666,507 Mapped reads 21,366,876 19,668,853 20,569,271 20,699,424 20,540,205 20,569,248 20,731,453 20,676,181 20,882,737 185,704,248 Mapping rate (%) 89.03 89.14 88.77 88.98 88.66 88.69 89.07 89.22 89.40 - Number of genes 19,137 19,058 19,162 19,042 19,045 19,048 18,977 18,993 19,017 - On the lookout for a subsequent successful GO enrichment analysis, a functional annotation of the turbot genome was carried out. 95.30% of the genes in the turbot genome were functionally annotated, summarized in three categories: BP, MF and CC. A total of 286,299 different GO terms were assigned, with 53.15%, 22.46% and 24.38%, corresponding to the BP, MF and CC categories, respectively. This functionally annotated turbot genome was used as a background to perform GO enrichment for selected genes of interest. Figure 1. Turbot developmental stages selected during the metamorphosis process for transcriptome analysis. From left to right: ( a ) pre-metamorphic symmetrical larva prior to eye migration (stage 3b: 15dpf) ( b ) asymmetrical larva at the metamorphic climax with upper edge of migrating right eye visible from left side (stage 4d: 30dpf) ( c ) post-metamorphic asymmetrical juvenile with upper eye entirely placed on the left side (stage 5c: 57dpf). Fish were reared at 18 ◦ C. Scale bars (a) 2mm, (b,c) 5 mm. Biology 2021, 10, x FOR PEER REVIEW 6 of 18 Figure 2. Overview of pipeline strategy followed in this work, from brain dissection to GO analysis of DEGs. 3.2. Gene Expression Dynamics in Metamorphosing Turbot Brains To identify DEGs in the turbot brain during metamorphosis, we performed pairwise comparisons among three postembryonic developmental stages (pre-metamorphic vs. climax, climax vs. post-metamorphic and pre-metamorphic vs. post-metamorphic). A total of 1570 genes were differentially expressed in the three metamorphic developmental stages. A high proportion of DEGs were found when pre-metamorphic vs. climax stages (338 up- and 221 down-regulated) and pre-metamorphic vs. post-metamorphic stages (415 up- and 487 down-regulated) were compared (Figures 3a and 3c, respectively; Supplementary Tables S3 and S4). However, a significantly lower number of DEGs (33 up- and 76 down-regulated) was found after comparison between the climax and post-metamor- phic stages (Figure 3b; Supplementary Table S5). All DEGs from the three developmental stages were combined into a single set and hierarchically clustered within a heatmap in order to produce an overview of the gene expression profiles across metamorphosis (Figure 3d). In the heatmap, DEGs were clustered according to gene expression level. Overall, two major clusters could be observed. Most genes clustered on the top half of the heatmap displayed decreased expression throughout the metamorphosis process (Figure 3d), while groups of genes on the bottom half showed increased expression over time. Specifically, a set of tightly clustered DEGs showed a marked expression peak at the metamorphic climax, while exhibiting down-regulated gene expression at the pre- and post-metamor- phic stages (Figure 3d). Figure 2. Overview of pipeline strategy followed in this work, from brain dissection to GO analysis of DEGs. On the lookout for a subsequent successful GO enrichment analysis, a functional annotation of the turbot genome was carried out. 95.30% of the genes in the turbot genome were functionally annotated, summarized in three categories: BP, MF and CC. A total of Biology 2021,10, 1256 6 of 16 286,299 different GO terms were assigned, with 53.15%, 22.46% and 24.38%, corresponding to the BP, MF and CC categories, respectively. This functionally annotated turbot genome was used as a background to perform GO enrichment for selected genes of interest. 3.2. Gene Expression Dynamics in Metamorphosing Turbot Brains To identify DEGs in the turbot brain during metamorphosis, we performed pairwise comparisons among three postembryonic developmental stages (pre-metamorphic vs. climax, climax vs. post-metamorphic and pre-metamorphic vs. post-metamorphic). A total of 1570 genes were differentially expressed in the three metamorphic developmental stages. A high proportion of DEGs were found when pre-metamorphic vs. climax stages (338 up- and 221 down-regulated) and pre-metamorphic vs. post-metamorphic stages (415 up- and 487 down-regulated) were compared (Figure 3a,c, respectively; Supplementary Tables S3 and S4). However, a significantly lower number of DEGs (33 up- and 76 down-regulated) was found after comparison between the climax and post-metamorphic stages (Figure 3b; Supplementary Table S5). All DEGs from the three developmental stages were combined into a single set and hierarchically clustered within a heatmap in order to produce an overview of the gene expression profiles across metamorphosis (Figure 3d). In the heatmap, DEGs were clustered according to gene expression level. Overall, two major clusters could be observed. Most genes clustered on the top half of the heatmap displayed decreased expression throughout the metamorphosis process (Figure 3d), while groups of genes on the bottom half showed increased expression over time. Specifically, a set of tightly clustered DEGs showed a marked expression peak at the metamorphic climax, while exhibiting down-regulated gene expression at the pre- and post-metamorphic stages (Figure 3d). We next focused on the selected gene clusters that exhibited specific expression at pre-metamorphic stage (Supplementary Table S6), climax stage (Supplementary Table S7) and post-metamorphic stage (Supplementary Table S8) and a GO term enrichment analysis was performed to investigate the biological functions associated with each gene dataset (Figure 4), using the previous functional annotation of the turbot genome as a background. In this analysis, the most significant GO terms (padj < 0.05) were identified. Genes from the cluster showing overexpression at the pre-metamorphic stage revealed a significant enrichment of BP GO terms related to the developmental process and the development of anatomical structures, such as anterior/posterior pattern specification (GO:0009952), embryonic skeletal system (GO:0048704) and multicellular organism development (GO:0007275). Regarding the MF ontology, the most significant GO terms included heme binding (GO:0020037) and monooxygenase activity (GO:0004497) and in the CC category, we found ontologies such as RNA polymerase II transcription regulator complex (GO:009075) and apical plasma membrane (GO:0016324) (Figure 4a; Supplementary Table S9). Genes that were predominantly expressed during the climax of metamorphosis enrich BP ontologies related to immune system processes, such as antigen processing and presentation of endogenous peptide antigen via MHC class I (GO:00019885), innate immune response (GO:0045087) and immune response (GO:0006958), in MF category the most significant GO terms were NAD+ ADP-ribosyltransferase activity (GO:0003950), extracellular matrix structural constituent (GO:0005201) and protein transmembrane transporter activity (GO:0008320) and in the CC category we found ontologies such as collagen trimer (GO:0005581), collagen type I primer (GO:0005584) and extracellular matrix (GO:0031012) (Figure 4b; Supplementary Table S10). We observed that clustered genes with a higher expression in the post-metamorphosis stage also enrich the BP ontologies related to immune system process, such as immune response (GO:0006955), although we can also highlight ontologies related to cellular components organization, such as extracellular matrix organization (GO:0030198) or muscle system processes, like the cardiac muscle contraction (GO:0060048) and skeletal muscle contraction (GO:0003009) (Figure 4c; Supplementary Table S11). Biology 2021,10, 1256 7 of 16 Biology 2021, 10, x FOR PEER REVIEW 7 of 18 Figure 3. Visualization of gene expression datasets in turbot brains throughout the developmental metamorphosis process. MA plot of all transcriptome genes from pairwise comparisons: (a) pre-metamorphic vs. climax, (b) climax vs. post-meta- morphic and (c) pre-metamorphic vs. post-metamorphic. The red dots plotted represent genes with an adjusted p-value < 0.1, while gray dots are those genes that do not show the established significance between the different stages. (d) Heatmap displaying the DEGs hierarchically clustered according to the expression profiles throughout metamorphosis. Each column represents an individual triplicate from each metamorphic stage (pre-metamorphic, climax and post-metamorphic) and each row represents different DEGs. The colors from light blue to dark blue indicate gene expression from low to high, respectively. We next focused on the selected gene clusters that exhibited specific expression at pre-metamorphic stage (Supplementary Table S6), climax stage (Supplementary Table S7) and post-metamorphic stage (Supplementary Table S8) and a GO term enrichment analysis was performed to investigate the biological functions associated with each gene dataset (Figure 4), using the previous functional annotation of the turbot genome as a background. In this analysis, the most significant GO terms (padj < 0.05) were identified. Genes from the cluster showing overexpression at the pre-metamorphic stage revealed a significant enrichment of BP GO terms related to the developmental process and the development of anatomical structures, such as anterior/posterior pattern specification (GO:0009952), embryonic skeletal system (GO:0048704) and multicellular organism development (GO:0007275). Regarding the MF ontology, the most significant GO terms included heme binding (GO:0020037) and monooxygenase activity (GO:0004497) and in the CC category, we found ontologies such as RNA polymerase II transcription regulator complex (GO:009075) and apical plasma membrane (GO:0016324) (Figure 4a; Supplementary Table S9). Genes that were predominantly expressed during the climax of metamorphosis enrich BP ontologies related to immune system processes, such as antigen processing and presentation of endogenous peptide antigen via MHC class I (GO:00019885), Figure 3. Visualization of gene expression datasets in turbot brains throughout the developmental metamorphosis process. MA plot of all transcriptome genes from pairwise comparisons: ( a ) pre-metamorphic vs. climax, ( b ) climax vs. postmetamorphic and ( c ) pre-metamorphic vs. post-metamorphic. The red dots plotted represent genes with an adjusted p-value < 0.1, while gray dots are those genes that do not show the established significance between the different stages. ( d ) Heatmap displaying the DEGs hierarchically clustered according to the expression profiles throughout metamorphosis. Each column represents an individual triplicate from each metamorphic stage (pre-metamorphic, climax and post-metamorphic) and each row represents different DEGs. The colors from light blue to dark blue indicate gene expression from low to high, respectively. 3.3. qPCR Validation of Gene Expression Patterns Real-time qPCR was used to validate the expression levels of chit3,epd,thra,tshba and ifih1 genes across turbot metamorphosis (Figure 5d). The qPCR results of the analyzed genes were consistent with the RNA sequencing data. We found that chit3,ifih1 and thra expression levels peaked at the metamorphic climax and decreased after overcoming metamorphosis (Figure 5d). In addition, the expression of tshba increased at the climax stage, but remained high in the post-metamorphic stage. The only gene that exhibited a gradual increase of expression from the pre-metamorphic to post-metamorphic stage was epd. (Figure 5d). Biology 2021,10, 1256 8 of 16 Biology 2021, 10, x FOR PEER REVIEW 8 of 18 innate immune response (GO:0045087) and immune response (GO:0006958), in MF category the most significant GO terms were NAD+ ADP-ribosyltransferase activity (GO:0003950), extracellular matrix structural constituent (GO:0005201) and protein transmembrane transporter activity (GO:0008320) and in the CC category we found ontologies such as collagen trimer (GO:0005581), collagen type I primer (GO:0005584) and extracellular matrix (GO:0031012) (Figure 4b; Supplementary Table S10). We observed that clustered genes with a higher expression in the post-metamorphosis stage also enrich the BP ontologies related to immune system process, such as immune response (GO:0006955), although we can also highlight ontologies related to cellular components organization, such as extracellular matrix organization (GO:0030198) or muscle system processes, like the cardiac muscle contraction (GO:0060048) and skeletal muscle contraction (GO:0003009) (Figure 4c; Supplementary Table S11). Figure 4. Main biological processes, molecular functions and cellular components boosted during metamorphosis in turbot brain. Gene Ontology (GO) enrichment of genes from hierarchical clustering showing high expression at (a) pre-metamorphic, (b) climax or (c) post-metamorphic stages. The most enriched GO terms (top 10, adjusted p-value < 0.05) belonging to biological process (BP), molecular function (MF) and cellular component (CC) categories are represented. Figure 4. Main biological processes, molecular functions and cellular components boosted during metamorphosis in turbot brain. Gene Ontology (GO) enrichment of genes from hierarchical clustering showing high expression at ( a ) premetamorphic, ( b ) climax or ( c ) post-metamorphic stages. The most enriched GO terms (top 10, adjusted p-value < 0.05) belonging to biological process (BP), molecular function (MF) and cellular component (CC) categories are represented. Biology 2021,10, 1256 9 of 16 Biology 2021, 10, x FOR PEER REVIEW 9 of 18 3.3. qPCR Validation of Gene Expression Patterns Real-time qPCR was used to validate the expression levels of chit3, epd, thra, tshba and ifih1 genes across turbot metamorphosis (Figure 5d). The qPCR results of the analyzed genes were consistent with the RNA sequencing data. We found that chit3, ifih1 and thra expression levels peaked at the metamorphic climax and decreased after overcoming metamorphosis (Figure 5d). In addition, the expression of tshba increased at the climax stage, but remained high in the post-metamorphic stage. The only gene that exhibited a gradual increase of expression from the pre-metamorphic to post-metamorphic stage was epd. (Figure 5d). Figure 5. Selected clusters from Mfuzz soft clustering analysis of DEGs in turbot brains throughout the developmental metamorphosis process. DEGs show up-regulated genes corresponding to the: (a) pre-metamorphic stage, (b) climax stage, (c) post-metamorphic stage. Color code, from magenta to yellow, denote high or low Mfuzz membership values, respectively. Time 1, 2 and 3 in X-axis corresponds to pre-metamorphic, climax and post-metamorphic stage, respectively. (d) Validation of expression patterns of chit3, epd, thra, tshba, and ifih1 genes in the turbot brain using qPCR. In the bar plots, the trend lines represent the fold change obtained by analyzing the RNAseq values during the metamorphic stages. The bars represent the fold change values obtained by qPCR. Results were normalized to 18S gene and expressed as the mean ± SEM of two independent experiments. Data from climax stage was set at 1. 3.4. Stage-Specific Gene Expression and GO Enrichment Analysis To further identify genes in the turbot brain with potentially important roles throughout the developmental metamorphosis process, we first explored the temporal expression patterns of DEGs by soft clustering analysis. The 1570 DEGs were divided into 28 clusters according to the similarity of their expression patterns throughout the three selected key stages. We then selected clusters containing genes that showed consistent and specific up-regulated expression at each metamorphic stage (Figure 5a–c). Clusters with genes highly expressed in the pre-metamorphic stage and then down-regulated displayed a total of 236 genes (Figure 5a; Supplementary Table S12), while 23 genes from a single cluster were specifically expressed in the post-metamorphic stage (Figure 5c; Supplementary Table S13). Interestingly, we found two clusters contained a total of 63 genes Figure 5. Selected clusters from Mfuzz soft clustering analysis of DEGs in turbot brains throughout the developmental metamorphosis process. DEGs show up-regulated genes corresponding to the: ( a ) pre-metamorphic stage, ( b ) climax stage, ( c ) post-metamorphic stage. Color code, from magenta to yellow, denote high or low Mfuzz membership values, respectively. Time 1, 2 and 3 in X-axis corresponds to pre-metamorphic, climax and post-metamorphic stage, respectively. ( d ) Validation of expression patterns of chit3, epd, thra, tshba, and ifih1 genes in the turbot brain using qPCR. In the bar plots, the trend lines represent the fold change obtained by analyzing the RNAseq values during the metamorphic stages. The bars represent the fold change values obtained by qPCR. Results were normalized to 18S gene and expressed as the mean ± SEM of two independent experiments. Data from climax stage was set at 1. 3.4. Stage-Specific Gene Expression and GO Enrichment Analysis To further identify genes in the turbot brain with potentially important roles throughout the developmental metamorphosis process, we first explored the temporal expression patterns of DEGs by soft clustering analysis. The 1570 DEGs were divided into 28 clusters according to the similarity of their expression patterns throughout the three selected key stages. We then selected clusters containing genes that showed consistent and specific up-regulated expression at each metamorphic stage (Figure 5a–c). Clusters with genes highly expressed in the pre-metamorphic stage and then down-regulated displayed a total of 236 genes (Figure 5a; Supplementary Table S12), while 23 genes from a single cluster were specifically expressed in the post-metamorphic stage (Figure 5c; Supplementary Table S13). Interestingly, we found two clusters contained a total of 63 genes whose expression peaks at the metamorphic climax (Figure 5b; Supplementary Table S14), which supports the data represented on the heatmap (Figure 3d). Additionally, we performed a GO enrichment analyses of the selected clusters (Figure 6). In the pre-metamorphic stage, the gene set from clusters selected showed significant enrichment of early development-related BP GO terms, such as embryonic skeletal system morphogenesis (GO:0048704), cholesterol homeostasis (GO:0042632) and anterior/posterior pattern specification (GO:0009952). For MF ontology, the most significant GO terms were UDP-glycosyltransferase activity (GO:0008194), hormone activity (GO:0005179) and heme Biology 2021,10, 1256 16 of 16 30. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013,29, 15–21. [CrossRef] [PubMed] 31. Anders, S.; Pyl, P.T.; Huber, W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics 2015 , 31, 166–169. [CrossRef] 32. Muñoz-Mérida, A.; Viguera, E.; Claros, M.G.; Trelles, O.; Pérez-Pulido, A.J. Sma3s: A three-step modular annotator for large sequence datasets. DNA Res. 2014,21, 341–353. [CrossRef] 33. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014,15, 550. [CrossRef] [PubMed] 34. Benjamini, Y.; Hochberg, Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J. R. Stat. Soc. Ser. B 1995,57, 289–300. [CrossRef] 35. Kumar, L.; Futschik, M.E. Mfuzz: A software package for soft clustering of microarray data. Bioinformation 2007 ,2, 5–7. [CrossRef] 36. Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. ClusterProfiler: An R package for comparing biological themes among gene clusters. Omi. A J. Integr. Biol. 2012,16, 284–287. [CrossRef] 37. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001 ,29, 2002–2007. [CrossRef] 38. Zoologie, R.C.; De Mazurais, D.; Darias, M.; Cahu, C.L. Transcriptomics for understanding marine fish larval development. Can. J. Zool. 2011,89, 599–611. 39. Kavouras, M.; Malandrakis, E.E.; Golomazou, E.; Konstantinidis, I.; Blom, E.; Palstra, A.P.; Anastassiadis, K.; Panagiotaki, P.; Exadactylos, A. Hox gene expression profiles during embryonic development of common sole. Anim. Biol. 2019 ,69, 183–198. [CrossRef] 40. Mallo, M.; Wellik, D.M.; Deschamps, J. Hox genes and regional patterning of the vertebrate body plan. Dev. Biol. 2010 ,344, 7–15. [CrossRef] [PubMed] 41. Cunha, I.; Galante-Oliveira, S.; Rocha, E.; Urbatzka, R.; Castro, L.F.C. Expression of intercellular lipid transport and cholesterol metabolism genes in eggs and early larvae stages of turbot, Scophthalmus maximus, a marine aquaculture species. Mar. Biol. 2015,162, 1673–1683. [CrossRef] 42. Cunha, I.; Galante-Oliveira, S.; Rocha, E.; Planas, M.; Urbatzka, R.; Castro, L.F.C. Dynamics of PPARs, fatty acid metabolism genes and lipid classes in eggs and early larvae of a teleost. Comp. Biochem. Physiol. B Biochem. Mol. Biol. 2013 ,164, 247–258. [CrossRef] 43. Larsen, D.A.; Swanson, P.; Dickey, J.T.; Rivier, J.; Dickhoff, W.W. In vitro thyrotropin-releasing activity of corticotropin-releasing hormone-family peptides in coho salmon, Oncorhynchus kisutch. Gen. Comp. Endocrinol. 1998,109, 276–285. [CrossRef] 44. Gao, C.; Cai, X.; Zhang, Y.; Su, B.; Song, H.; Wenqi, W.; Li, C. Characterization and expression analysis of chitinase genes (CHIT1, CHIT2 and CHIT3) in turbot (Scophthalmus maximus L.) following bacterial challenge. Fish Shellfish Immunol. 2017 ,64, 357–366. [CrossRef] 45. Teng, Z.; Sun, C.; Liu, S.; Wang, H.; Zhang, S. Functional characterization of chitinase-3 reveals involvement of chitinases in early embryo immunity in zebrafish. Dev. Comp. Immunol. 2014,46, 489–498. [CrossRef] 46. Brisse, M.; Ly, H. Comparative structure and function analysis of the RIG-I-like receptors: RIG-I and MDA5. Front. Immunol. 2019 , 10, 1–27. [CrossRef] 47. Sadler, A.J. The role of MDA5 in the development of autoimmune disease. J. Leukoc. Biol. 2017,103, 185–192. [CrossRef] 48. Morimoto, K.; Nakajima, K. Role of the Immune System in the Development of the Central Nervous System. Front. Neurosci. 2019,13. [CrossRef] [PubMed] 49. Orozco, S.; Oberst, A. RIPK3 in cell death and inflammation: The Good, the Bad, and the Ugly. Immunol. Rev. 2017 ,277, 102–112. [CrossRef] [PubMed] 50. Numajiri Haruki, A.; Naito, T.; Nishie, T.; Saito, S.; Nagata, K. Interferon-inducible antiviral protein MxA enhances cell death triggered by endoplasmic reticulum stress. J. Interf. Cytokine Res. 2011,31, 847–856. [CrossRef] 51. Andersen, M.H.; Schrama, D.; Thor Straten, P.; Becker, J.C. Cytotoxic T cells. J. Invest. Dermatol. 2006 ,126, 32–41. [CrossRef] [PubMed] 52. Yu, J.; Zhou, X.; Chang, M.; Nakaya, M.; Chang, J.-H.; Xiao, Y.J.; Lindsey, W.; Dorta-Estremera, S.; Cao, W.; Zal, N.; et al. Regulation of T-cell activation and migration by the kinase TBK1 during neuroinflammation. Nat. Commun. 2015,6, 6074. [CrossRef] 53. Rollins-Smith, L.A. Metamorphosis and the amphibian immune system. Immunol. Rev. 1998 ,166, 221–230. [CrossRef] [PubMed] 54. Pollet, N. Expression of immune genes during metamorphosis of Xenopus: A survey. Front. Biosci. 2010 ,15, 348–358. [CrossRef] 55. Roelfsema, F.; Boelen, A.; Kalsbeek, A.; Fliers, E. Regulatory aspects of the human hypothalamus-pituitary-thyroid axis. Best Pract. Res. Clin. Endocrinol. Metab. 2017,31, 487–503. [CrossRef] [PubMed] 56. Khansari, N.; Shakiba, Y.; Mahmoudi, M. Chronic Inflammation and Oxidative Stress as a Major Cause of Age- Related Diseases and Cancer. Recent Pat. Inflamm. Allergy Drug Discov. 2009,3, 73–80. [CrossRef] 57. Carbone, D.; Faggio, C. Importance of prebiotics in aquaculture as immunostimulants. Effects on immune system of Sparus aurata and Dicentrarchus labrax. Fish Shellfish Immunol. 2016,54, 172–178. [CrossRef] [PubMed] 58. Remen, M.; Nederlof, M.A.J.; Folkedal, O.; Thorsheim, G.; Sitjà-Bobadilla, A.; Pérez-Sánchez, J.; Oppedal, F.; Olsen, R.E. Effect of temperature on the metabolism, behaviour and oxygen requirements of Sparus aurata. Aquac. Environ. Interact. 2015 ,7, 115–123. [CrossRef]