A comprehensive coding and microRNA transcriptome of vertebral bone in postlarvae and juveniles of Senegalese sole (Solea senegalensis)

Genomics 116 (2024) 110802 Available online 28 January 2024 0888-7543/© 2024 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). A comprehensive coding and microRNA transcriptome of vertebral bone in postlarvae and juveniles of Senegalese sole (Solea senegalensis) Carmen Bouza a , * , Ana P. Losada b , Carlos Fern´ andez a , Jos´ e A. ´ Alvarez-Dios c , Ana Manuela de Azevedo b , Andr´ es Barreiro b , Dami´ an Costas d , María Isabel Quiroga b , Paulino Martínez a , Sonia V´ azquez b a Department of Zoology, Genetics and Physical Anthropology, Faculty of Veterinary, Campus Terra, Universidade de Santiago de Compostela, 27002 Lugo, Spain b Department of Anatomy, Animal Production and Veterinary Clinical Sciences, Faculty of Veterinary, Campus Terra, Universidade de Santiago de Compostela, 27002 Lugo, Spain c Department of Applied Mathematics, Faculty of Mathematics, Universidade de Santiago de Compostela, 15705 Santiago de Compostela, Spain d Centro de Investigaci´ on Mari˜ na, Universidade de Vigo, ECIMAT, Vigo 36331, Spain ARTICLE INFO Keywords: Solea senegalensis Vertebral bone Spinal anomalies Transcriptomics mRNA miRNA ABSTRACT Understanding vertebral bone development is essential to prevent skeletal malformations in farmed fish related to genetic and environmental factors. This is an important issue in Solea senegalensis, with special impact of spinal anomalies in postlarval and juvenile stages. Vertebral bone transcriptomics in farmed fish mainly comes from coding genes, and barely on miRNA expression. Here, we used RNA-seq of spinal samples to obtain the first comprehensive coding and miRNA transcriptomic repertoire for postlarval and juvenile vertebral bone, covering different vertebral phenotypes and egg-incubation temperatures related to skeleton health in S. senegalensis. Coding genes, miRNA and pathways regulating bone development and growth were identified. Differential transcriptomic profiles and suggestive mRNA-miRNA interactions were found between postlarvae and juveniles. Bone-related genes and functions were associated with the extracellular matrix, development and regulatory processes, calcium binding, retinol and lipid metabolism or response to stimulus, including those revealed by the miRNA targets related to signaling, cellular and metabolic processes, growth, cell proliferation and biological adhesion. Pathway enrichment associated with fish skeleton were identified when comparing postlarvae and juveniles: growth and bone development functions in postlarvae, while actin cytoskeleton, focal adhesion and proteasome related to bone remodeling in juveniles. The transcriptome data disclosed candidate coding and miRNA gene markers related to bone cell processes, references for functional studies of the anosteocytic bone of S. senegalensis. This study establishes a broad transcriptomic foundation to study healthy and anomalous spines under early thermal conditions across life-stages in S. senegalensis, and for comparative analysis of skeleton homeostasis and pathology in fish and vertebrates. 1. Introduction The Senegalese sole (Solea senegalensis) is a valuable flatfish for European aquaculture diversification. Farming production raised up to 2700 t in the last years [1,2]. The main priorities for Senegalese sole farming industry include the optimization of reproduction, nutrition, growth, disease control, and, in general, fish quality. The incidence of skeletal malformations represents a major constrain for its production, ranging between 40% and 100% from experimental and farm data [1,3–5], which have been associated with decreasing fish quality and welfare [1,6–8]. Spinal disorders are among the most prevalent skeletal anomalies in Senegalese sole affecting quality at postlarval and juvenile stages just before transfer to on-growing farms [9–11]. Pathology methods have allowed assessment of skeleton health and diagnosis of different types of vertebral anomalies (i.e., fusions, compressions, axis deviations), which can remain externally undetected at earlier stages that may progress later to more severe pathologies or show anomalous vertebral tissue organization [3,4,9,10,12–14]. Studies on bone metabolism and development underlying skeletal homeostasis revealed a multifactorial aetiology of anomalies in fish * Corresponding author. E-mail address: [email protected] (C. Bouza). Contents lists available at ScienceDirect Genomics journal homepage: www.elsevier.com/locate/ygeno https://doi.org/10.1016/j.ygeno.2024.110802 Received 15 September 2023; Received in revised form 22 January 2024; Accepted 26 January 2024 Genomics 116 (2024) 110802 2 aquaculture, involving complex interactions between environmental, nutrition and genetic factors [1,5,7,15–18]. Among abiotic parameters, suboptimal temperatures during early development and/or rearing stages have been reported as major risk factors [19–24], including the Senegalese sole [25,26]. Genetic and genomic studies have provided information on the heritable component of skeleton malformations along with quantitative trait loci and candidate genes, pinpointing to differential associated markers and genomic regions across different types of skeletal anomalies [7,16,17,27,28]. Functional analyses have been focused on candidate coding genes related to fish skeletogenesis under different conditions, including thermal and nutritional factors (e. g., [19,29–31]). RNA sequencing (RNA-seq) for quantifying global expression profiles of coding and non-coding genes in aquaculture [32] has been scarcely applied to date for unveiling molecular mechanisms underlying vertebral bone development. Transcriptomes of spinal tissues have provided useful resources to unravel vertebral bone-related coding genes in farmed fish, revealing conserved regulatory mechanisms within vertebrates. Coding genes, signaling pathways and transcription factors related to osteoblastogenesis, chondrogenesis, matrix mineralization or bone resorption have been identified in a few farmed teleosts by stablishing reference transcriptomes of vertebral tissues from normal and/or anomalous fish under different rearing conditions [33–37]. Among these, fish rearing temperature has been related to transcriptomic changes in the vertebral column and bone-associated cells from subsets of candidate genes [19,20,22,23] or large-scale RNA-seq data [36]. Non-coding microRNAs (miRNA) represent major transcriptional regulators of biological processes related to bone development and pathology in vertebrates, including model teleosts [38–42]. In farmed fish, large-scale miRNA transcriptomic inventories have provided useful references across life-cycle stages laying the foundation for functional and comparative studies of productive traits, such as growth, myogenesis and immunology [43–47], but to a much lesser extent for vertebral bone development and pathology. The increasing genomic and transcriptomic resources in Senegalese sole have provided insights into the genetic basis of production traits (e. g., larval development, reproduction, growth, response to pathogens) for industrial application [48–51]. Some skeleton-related coding genes have been assayed from larval stages to adulthood under different conditions in this flatfish with acellular (anosteocytic) bone. Liposoluble vitamins and signaling pathways during metamorphosis has been explored on the origin of bone anomalies [52–55], as well as head transcriptomics on retinoid signaling pathway during metamorphosis [49,56]. However, to our knowledge, no large-scale RNA-seq data have been generated from spinal samples of Senegalese sole for studying the vertebral bone development and homeostasis in postlarval and juvenile stages, susceptible to vertebral anomalies [5,10], which have been associated with exposition to different embryonic thermal conditions in this species [25,26]. Temperature during early development has been related to long-term effects on gene expression in fish [44] and also to changes in the expression of growth-related coding genes and miRNA in the Senegalese sole [57,58]. This study aimed to generate large-scale transcriptome resources of vertebral bone in postlarval and juvenile stages of the Senegalese sole to establish the repertory of expressed bone-related coding and miRNA genes in this species. Different spinal phenotypes and egg-incubation temperatures related to skeleton health were analysed to broaden transcriptome RNA-seq data. Candidate bone-related genes, functions and pathways were explored and compared between postlarvae and juveniles to characterize the vertebral transcriptome in early life stages of the species, pursued to facilitate future studies on vertebral bone development and pathology in response to early thermal conditions in Senegalese sole and fish. 2. Material and methods 2.1. Biological sampling Hatchery-reared Senegalese sole was sampled at postlarval and juvenile stages where non-invasive pathology methods can be applied for vertebral phenotype analysis, relevant for product quality in on-growing phases [10,13]. Postlarvae and juveniles proceeded from a random sample of Senegalese sole fertilized eggs from a commercial batch supplied by the company Stolt Sea Farm S.A. (Merexo, A Coru˜ na, Spain) and transferred to the ECIMAT facilities (Marine Research Center, Vigo, Spain; https://cim.uvigo.gal/en/), where the culture trial was performed following the European Directive 2010/63/EU. Eggs were incubated for 24 h in 10 L conic tanks at two thermal conditions (18 ◦C and 22 ◦C), which had been previously associated in this species with differential prevalence of skeletal anomalies, higher at 22 ◦C [25,26] and with long-term effects on gene expression [57,58]. Hatched larvae were distributed in triplicate tanks per egg incubation condition (18 ◦C and 22 ◦C), to ensure persistence of culture design and sampling collection for histological/genomic analysis. All the six tanks were maintained at the same temperature (mean 19.1 ◦C, SD 0.40) until the end of the trial. Standard ECIMAT rearing protocols were applied, in line with what was reported for industrial farm conditions by de Azevedo et al. [5]. Benthic larvae at 14 days after hatching (dah) were transferred to 100 L square tanks, keeping three replicates per egg-incubation condition, at a density of ~30 larvae/L under 13 h:11 h light:dark cycle, and dissolved oxygen of 6.9 mg/L (SD 0.4). Feeding was adjusted to the development stage. Larvae were fed with rotifers until 8 dah and Artemia was added from 6 dah to 31 dah. Fish were fed with commercial inert feed (Gemma Micro; Skretting) from 32 dah to the end of the trial. Random samples of 31 dah postlarvae and 123 dah juveniles were collected for vertebral phenotyping and RNA extraction. Fish were sacrificed using an overdose of Tricaine methanesulfonate (MS-222, Sigma-Aldrich, Madrid, Spain). Vertebral phenotyping consisted in assessing the vertebral body morphology, to detect normal and anomalous vertebrae and/or axis deviations along the spine in postlarvae and juveniles, considering three anatomical regions (abdominal and caudal vertebrae, and caudal complex) according to de Azevedo et al. [5] (Fig. 1A). The vertebral phenotype at 31 dah was explored and recorded using a binocular stereomicroscope model Olympus® (Barcelona, Spain) SZX16 coupled to a digital camera Olympus® DP72. Juvenile vertebral phenotypes were explored using a portable X-ray radiographic machine model Trophy 708-CCX (Trophy Radiologie, France), recorded with a Kodak® RVG 6100 Digital Radiography System (Carestream Health, USA). Further image analysis and processing was performed with ImageJ 1.50i software (National Institutes of Health, Bethesda, MD, USA). For each region, the vertebral and caudal complex elements were evaluated, focusing on vertebral bodies, as described in previous studies [9,10,13]. Fish were classified attending to the presence (and type) or absence of anomalies (Fig. 1B-C). Samples of the vertebral column from postlarvae and segments of the vertebral column from juveniles were collected and immersed in TRIzol Reagent (Life Technologies) for further analyses. The juvenile spinal segments contained the vertebrae at similar anatomical regions that exhibited some malformation or did not show anomalies (Supplementary Table S1). 2.2. RNA sequencing Postlarval and juvenile vertebral bone samples were sequenced to capture stage-specific gene expression profiles, as reported across different life-stages in fish [43,46,77,79]. Different spinal phenotypes and temperature conditions were included to increase transcriptome diversity according to previous reports [16,17,19,34,57,58]. Thus, three samples per life-stage and thermal condition were assayed using mRNA and small RNA (sRNA) sequencing, covering postlarval and juvenile C. Bouza et al. Genomics 116 (2024) 110802 3 stages after exposition to two egg-incubation temperatures (18 ◦C vs. 22 ◦C; 12 samples in total; Supplementary Table S1). Spinal anomalies were mostly found at 22 ◦C in the samples collected for our transcriptomic study, as reported [25,26]. Different spine phenotypes were included at each life-stage, half of the samples at 18 ◦C without anomalies or with minor preural fusions in the caudal complex, and the other half at 22 ◦C showing major vertebral anomalies, such as fusions, deformations or axis deviations in the caudal or abdominal regions [9,10,13]. Vertebral bone samples were dissected out and cleaned of adjacent tissues using a scalpel and a small brush under a binocular model Olympus® (Barcelona, Spain) SZX16. Each sample was stored and homogenized in TRIzol Reagent (Life Technologies), using a TissueLyser II (Qiagen, Germany) with one 3 mm zirconia bead for two cycles of 1 min at 30 Hz, and 1 min resting in-between. Total RNA was extracted using a TRIzol protocol, and RNA quantity and integrity evaluated with NanoDrop® ND −1000 spectrophotometer (NanoDrop® Technologies Inc.) and RNA 6000 Nano Bioanalyzer 2100 Assay (Agilent). The 12 RNA samples were delivered to Novogene-Europe, UK, to construct both poly-A enriched mRNA and sRNA libraries per sample (24 libraries in total) to be sequenced in 150 bp pair-end (PE) and 50 bp single-end (SE), respectively, on an Illumina NovaSeq 6000 S2 platform. 2.3. RNA-seq data processing and annotation Quality control of mRNA raw reads from 12 libraries of vertebral bone was assessed with FastQC v.0.11.7 (www.bioinformatics.babraha m.ac.uk/projects/fastqc/). The Harvard FAS Informatics pipeline was used for quality filtering [59], discarding read pairs for which one of the reads was deemed unfixable, removing residual adaptor sequences and erroneous k-mers. Clean reads were aligned against the Senegalese sole genome (NCBI Genome Assembly: GCA_919967415.2; [51]) with STAR 2.7.9 [60] with –quantMode to obtain GeneCounts per sample. Gene expression was recorded when at least two samples showed a TPM (Transcripts Per Kilobase Million) above 4. Annotation from the gff file was completed by searching FASTA sequences in Swissprot using Blastn [61]. GO functional annotations for each sequence were obtained by blasting FASTA sequences against the uniref90 database; the best homologous sequence and its name were used to extract the GO terms at the UniProt web portal (https://www.uniprot.org) using a custom Perl script. Transcriptome functional analysis was made using WEGO [62] and Enrichment Analysis performed in AgriGO v2.0 [63] using a 5% false discovery rate (FDR). KEGG pathways were analysed using KOBASi [64] and enrichment assessed using a Fisher’s Exact Test and FDR correction (Benjamini-Hochberg). FASTA clean reads from 12 sRNA libraries of vertebral bone were processed for miRNA analysis using Prost! V0.7.60 [65], aligning the unique set of sequences to the reference Senegalese sole genome (GCA_919967415.2 [51]). The stacked reads were annotated by aligning against the mature and hairpin sequences of Oryzias latipes, a closely related acanthopterygian species with high quality miRNA annotation using FishmiRNA database [47]. We configured Prost! to select for reads between 17 and 25 nucleotides long and retained only sequences with a minimum of five identical reads for the initial annotation, and further retaining only sequences with a minimum of 30 reads for differential expression analysis. Clusters of miRNAs were explored in the Senegalese sole genome, when located at <10 kb distance on the same strand [66]. 2.4. Differential gene expression between postlarval and juvenile stages Differentially expressed (DE) genes for mRNA (DEG) and DE miRNA data were identified using DESeq2 v.1.34.0 [67] in R v.4.1.2. Absolute gene counts from STAR (mRNA) and from Prost! (miRNA) were normalized by DESeq2 to obtain log2 fold changes (FC) and adjusted pvalues (Benjamini-Hochberg FDR). Significant DE transcripts (FDR < Fig. 1. Methodology for vertebral phenotyping and sampling collection in Solea senegalensis. A) Both developmental stages, postlarva and juvenile, are depicted in photographs and diagrams illustrating the skeletal structure. The three anatomical regions of the vertebral column are identified in all images: abdominal region (AR), caudal region (CR), and caudal complex (CCo). The black box highlights the spinal regions of interest where samples were collected. B) Micrographs showing representative samples of normal vertebral structures. C) various types of anomalies detected in the conditions studied, including fusions (F), deformations (D), and deviations (DV). The scale bars for the macrographs are provided with their corresponding measurements. The black scale bars for the micrographs of the postlarvae span 0.5 mm. The white scale bars for the micrographs of the juveniles span 5 mm. C. Bouza et al. Genomics 116 (2024) 110802 4 0.05) were hereinafter referred to as upregulated or downregulated in postlarvae taken as reference juveniles (positive or negative FC values, respectively), or in juveniles taking postlarvae as reference, depending on the genes or functions analysed. Hierarchical clustering and principal component analyses (PCA) were performed to assess the clustering of samples and identify putative outliers. The R packages “pheatmap” and “PCAtools” were used to plot heatmaps and PCA, respectively. Functional analyses for DE transcripts were performed as previously described. Searching was conducted to identify candidate coding genes and miRNA related to bone development, and osteogenic, chondrogenic and osteoclastogenic processes, based on GO and KEGG functional analysis, along with previous studies in fish and vertebrates. 2.5. miRNA target prediction Because UTR sequences in the sole genome were not available for many genes, a conservative identification of the 3 ′ UTR sequences was set up by taking 700 bp at the end of coding genes in the Senegalese sole transcriptome [51], according to the median size across fish species [68]. Potential miRNA targets were predicted with MiRanda3.5 [69] with energy threshold-25 and other parameters by default. MiRanda output was pruned and processed to extract relevant information with a custom Perl script. Functional enrichment analysis for miRNA target genes against the reference transcriptome was explored using WEGO and KOBAS-i. 3. Results and discussion 3.1. RNA sequencing In this study, we used a combined mRNA and miRNA sequencing approach to obtain large-scale transcriptomic resources from vertebral bone samples of the Senegalese sole (Supplementary Table S1A and B), a valuable farm flatfish with high incidence of spinal anomalies [1,5]. For mRNA sequencing, 387 million PE raw reads were generated for the 12 spinal samples studied, and 98% of them were retained after quality filtering (on average, 31 million per sample). For miRNA analysis in the same vertebral samples, 143 million of 50 bp SE raw reads were generated, and 121 million reads (85%) were retained after filtering (on average, ~10 million reads per sample). 3.2. mRNA transcriptome A total of 10,591 genes were identified in the vertebral transcriptome of Senegalese sole, 96.4% annotated using the reference genome [51] (Supplementary Table S2A). Based on the median TPM over all samples, 960 highly expressed genes were identified (>100 TPM). The 20 most expressed genes in postlarvae and juveniles included key candidates related to fish vertebrae, such as col1a1 (collagen type I alpha 1 chain), a major component of bone matrix, or pvalb (parvalbumin), a calciumbinding protein associated with spinal anomalies [33,70] (Supplementary Table S2B). Some ribosomal protein genes (rps27; rpl37a) were highly expressed, especially in postlarvae, suggesting very active protein synthesis during earlier development stages. Other highly expressed genes in both stages are putatively associated with remnant adjacent tissues (e.g., neural, haemal or skeletal muscle, such as myelin basic protein; haemoglobin proteins; creatine kinase M-type, glycolytic enzymes; Supplementary Table S2B), as expected for small spinal samples in fish [33]. The whole transcriptome catalogue of postlarval and juvenile stages revealed a clear differentiation in the PCA in the conditions studied (Supplementary Fig. S1A), allowing exploring stage-specific sets of expressed genes for comparative studies across life-stages [33,76–78]. GO terms were retrieved for 94.7% of annotated genes, revealing a wide range of functions. Prominent functions included: biological processes (BP), such as cellular and metabolic process, biological regulation, response to stimulus, signaling, biological adhesion, and growth; molecular functions (MF), such as binding and catalytic activity, transport, transducer activity and transcription regulation; and cellular components (CC), such as cell membrane and extracellular region (Supplementary Fig. S1B; Supplementary Table S2C). Functional enrichment of the vertebral bone mRNAome was inspected taking as reference the whole transcriptome of Senegalese sole (Supplementary Table S2C). GO terms related to potassium channel complex (CC), calcium and vitamin binding and catalytic activity (MF), and cellular and different metabolic processes (BP) were identified, pointing to noteworthy pathways, such as aromatic amino acids metabolism related to bone turnover and calcification [71]. 3.3. Differentially expressed coding genes (DEG) between stages A total of 4683 DEGs between postlarvae and juveniles were detected (Fig. 2A; Supplementary Table S3A), distributed across all Senegalese sole chromosomes (Chr) [51], although enriched at some of them (Chr12, Chr16 and Chr19) (Supplementary Fig. S1C). Mapping of vertebral mRNA transcripts and DEGs can provide useful positional candidates for comparative studies of skeletal development and growth in fish, further applied to dissect complex traits in aquaculture [17,27,44,72]. DEGs showed GO enrichment in terms such as regulation, transducer activity, response to stimulus, signaling and biological adhesion (Fig. 2B; Supplementary Table S3B). Among the top DEG, we could identify candidate genes related to the skeletal system and bone homeostasis, many of them upregulated in postlarvae, such as s100p, a S100 calcium-binding protein member involved in osteogenesis [73], while others were downregulated in postlarvae, such as col12a1, associated with osteoblasts response in mice [74] or casr (calcium sensing receptor), key for skeletal development in fish [75]. Functional enrichment of stricter subsets of DEGs between stages (−1 >FC >1; Supplementary Fig. S2A and B) revealed, among upregulated genes in postlarvae, MF terms related to compound binding, transport and regulator activity, and BP terms related to metabolic processes and growth, while among downregulated in postlarvae we identified response to stimulus, stress, and cell death. Enriched KEGG pathways including WNT, TGFβ or mTOR signaling related to growth and bone development were detected among upregulated genes in postlarvae, while actin cytoskeleton, focal adhesion, PPAR signaling or proteasome, more related to bone remodeling for downregulated in postlarvae, many of them previously associated with vertebral transcriptomes and spinal anomalies in fish [36,37]. These GO terms and KEGG pathways represent candidate functions on which to focus future studies on vertebral bone development and growth across lifecycle in flatfish. 3.4. Bone-related coding genes and pathways Candidate genes associated with skeletal system development, bone growth and homeostasis, involving bone cell functions were identified in the Senegalese sole vertebral transcriptome (Tables 1 and 2; Supplementary Table S2A), meeting to previous functional studies on fish spinal tissue [18,33,34,36,37,76]. Many DEGs between stages were included to be considered in future comparative studies of vertebral bone across lifecycle in Senegalese sole (Table 1; Supplementary Table S4). Most genes involved in morpho-skeletogenesis and development showed higher expression in postlarvae, whereas others related to extracellular matrix (ECM) increased their expression in juveniles, as reported in other fish [77]. Gene families associated with specific bonerelated functions were identified in the vertebral mRNAome of Senegalese sole, including ECM, development, regulatory processes, calcium binding, and other relevant pathways such as retinol, lipid metabolism or response to stimulus (Table 1; Fig. 3C and D; Supplementary Table S4). Among them, a baseline of selected genes and pathways was C. Bouza et al. Genomics 116 (2024) 110802 5 stablished for further studies of spinal bone development and pathology in Senegalese sole and other fish. ECM-related genes included conserved collagen members among vertebrates, such as col1a1 and col2a2, expressed in fish vertebrae centra [78], and col9a2 and col11a2, involved in spinal segmentation and intervertebral disk mineralization [79,80]. Some were related to control of bone formation and vertebral structure in farmed Atlantic salmon, such as col1a2 and col11a1 [81], or to early osteoblast maturation in medaka, such as col10a1 [82]. Also, ECM-related glycoproteins involved in skeletal maturation, such as thrombospondins (thbs3a, 4a, 4b), osteopontin (spp1, secreted phosphoprotein 1) and osteonectin (sparc, secreted protein acidic and cysteine rich), were detected, along with proteoglycans, such as aggrecan (acan), brevican (bcan), decorin (dcn) or lumican (lum), involved in matrix and collagen organization, inflammation, and cell signaling. Other ECM-related genes included fibronectin 1 (fn1), cartilage intermediate layer protein 2 (cilp2) and cartilage oligomeric matrix protein (comp), related to the chondrocyte pericellular matrix, bone development and homeostasis [33,83]. Matrix remodeling genes, including ADAM metallopeptidase with thrombospondin type 1 motif 14 (adamts14), related to bone formation and early development in zebrafish [84], and metalloproteinases related to bone homeostasis in early development and tissue remodeling, such as mmp2 or mmp13 [85], were also found. This latter involved in osteoclast activation and bone resorption was upregulated in juveniles in this study. Developmental pathways and regulatory functions in the Senegalese sole vertebrae included different nuclear receptors, such as nr4a1, regulator of osteoblast differentiation [76], other prominent regulators, such as parathyroid hormone and receptors (pthlhb, pth3r) involved in fish skeletal mineralization and phosphate regulation, and receptors for retinoic acid (rara, rarg) or vitamin D (vdrb), essential for bone formation and calcium homeostasis [76]. Receptors for thyroid hormone (thra, thrb), involved in bone growth [76] were upregulated in postlarvae, as also observed during early development and growth transitions in fish [86], while nfe2 like bzip transcription factor 2a (nfe2l2a), associated with oxidative stress in bone remodeling in fish vertebrae [87], was downregulated in postlarvae in the conditions assayed. Transcription factors involved in skeleton development during embryonic and larval stages were mostly upregulated in postlarvae, but also expressed in later juvenile stages, as reported in other fish [33,76,78]. These included sox and hox members involved in chondrogenesis (sox5, sox6, sox9) and osteoblastogenesis (sox8, sox11, hoxa10, hoxb13), and osteogenic-related genes, such as RUNX family 2 (runx2) and msh homeobox 2 (msx2), along with paired box genes (pax3, pax9) playing a role on axial skeleton development [76,88]. Among other interesting candidates, distal-less homeobox members (i.e., dlx15, dlx3b) expressed in fish chondroblasts and osteoblasts [82], T-box genes, such as tbx15, regulating skeleton development [76], twist BHLH members (twist1, twist2), involved in chondrogenesis, osteogenic inhibition and intermuscular bone development in fish [89,90], or the Fig. 2. Differentially expressed genes (DEG) from mRNA data between postlarva (plarva) and juvenile (juv) vertebrae of Solea senegalensis. A) Heatmap of significant DEGs (sample codes in Supplementary Table S1); B) Functional enrichment of DEGs (grey) respect to the vertebral bone transcriptome (red) (GO-term categories: CC cellular content, MF molecular function, BP biological process); C) and D) Candidate bone-related DEGs: Functional description and enrichment of KEGG pathways, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) C. Bouza et al. Genomics 116 (2024) 110802 6 osteoclastic-related hematopoietic regulator PU.1 (spi1), were identified. Conserved signaling crosstalk mechanisms were found in poslarvae and juvenile Senegalese sole vertebrae, including hedgehog signaling members (ihhb, hhip) involved in ossification and osteoblastogenesis, which have been related to hypertrophic chondrocytes and skeleton anomalies in fish [76]. Bone morphogenic proteins (BMP) and their receptors involved in mesenchymal cell proliferation, chondrogenesis and bone formation in fish were found (i.e., bmp1, bmp2, bmp4, bmp7b, bmpr1a, 1b, 2b) [33,91]), mostly upregulated in postlarvae in our study, along with BMP-activated MAD homeobox members (smad2 to smad5) and antagonists, such as chordin (chrd) and noggin 2 (nog2), critical for bone formation. We also identified members of TGFβ signaling pathway associated with skeletal development and bone homeostasis including tgfb3 and tgfbr2, involved in chondrogenesis and related to vertebral defects [35,92]. Fibroblast growth factors and receptors related to chondrogenesis, and ossification were also detected (e.g., fgf2, fgf8b, fgf10, fgf18; fgfr1 to fgfr3), including fgf23, a phosphate regulator related to the osteocytic function in cellular bones, as reported in other fish with acellular bone [33,76]. Other master growth regulators were found, such as igf1, involved in muscle and bone growth, osteoblastogenesis, bone turnover and bone response to water temperature in farmed fish [93]. In our study, igf1 and igf2 were down and upregulated in postlarvae; showing differential expression profiles during larval development of Senegalese sole [57] and other vertebrates [94]. In addition, members of the WNT signaling pathway, including major candidates related to osteoblastogenesis, chondrogenesis, bone metabolism and pathology [95], such as wnt11 [76], wnt3a and wnt4a or beta-catenin (ctnnb1) [33,35] were identified. We also found frizzle class receptor genes interacting with the WNT pathway including some candidates upregulated in postlarvae (fzd2, 4, 8, 9), related to osteogenic process, vertebrae development and mineralization [95]. Members of dickkopf WNT pathway inhibitors included dkk1, associated with osteoblastic and chondroblastic differentiation and deficient skeleton development, or dkk3, enhancing osteoclastogenesis [95], which was upregulated in juveniles in our study. Chemokine signaling members related to bone remodeling were also detected (e.g., ccl5, cxcl12, ccr1, cxcr4), including cxcr2 downregulated in postlarvae, and cxcr32, essential for bone resorption in medaka [96]. Key signaling mechanisms Fig. 3. Differentially expressed (DE) miRNA between postlarval (plarva) and juvenile (juv) vertebrae of Solea senegalensis. A) Heatmap of significant DE miRNA (sample codes in Supplementary Table S1); B) Functional enrichment for all predicted gene targets of the vertebral miRNAome (GO-terms categories: CC cellular content, MF molecular function, BP biological process); C) Predicted DEG targets of DE miRNA in the vertebral transcriptome (up and downregulated genes in postlarvae taken juveniles as reference are marked in red and blue colour, respectively). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) C. Bouza et al. Genomics 116 (2024) 110802 7 related to osteoclastogenesis included cytokine macrophage colonystimulating factor (csfb1) and members of the RANK/RANKL/OPG pathway regulating osteoclastic markers, such as tartrate-resistant acid phosphatase (acp5a) (Table 2). Other candidates involved in the retinol pathway, related to osteoblastogenesis, calcification and vertebral development in fish [5,76] were identified (e.g., alcohol dehydrogenase 1, adh1; aldehyde dehydrogenase 1 a2, aldh1a2; lecithin retinol acyltransferase, lrat; cytochrome p450 family 3 subfamily a40, cyp3a40), most of them downregulated in postlarvae in the conditions studied. Different heat shock protein genes, expressed in response to thermal and stressful conditions were also found, mostly downregulated in juvenile samples (e.g., hspa8, dnajc7, dnajc12, cryab, hsp27), interestingly related to bone metabolism and disease in mammals [97], and to long-term response to embryonic temperatures associated with growth in Senegalese sole [57]. Moreover, different upregulated genes in postlarvae for high-density lipoproteins (such as apoa2, apoa4, apoc1, apoda1, apodb) are related to enriched functions, such as cell adhesion, locomotion, and lipid metabolism, interacting with bone physiology and pathology [14,98]. Also, crucial members for calcium binding proteins were detected, as in other fish vertebrae [34] (e.g., pvalb4, s100a6, s100b, anxa2, anxa6), Table 1 Bone-related differentially expressed genes between postlarval and juvenile vertebrae of Solea senegalensis. Genes Coding gene description / miRNA (function) Chr: position (Mb) FC DEG (mRNA) bcan Brevican Chr10:1.6 -1.559 bgnb Byglican b Chr13:10.8 -1.603 col11a1b Collagen XI Alpha 1b Chr1:45.4 -0.777 col1a1a Collagen I Alpha 1a Chr15: 2.4 -1.227 col2a1a Collagen II Alpha 1a Chr13:23.5 -1.288 col9a1b Collagen IX Alpha 1b Chr18:14.8 -0.849 dcr Decorin Chr2:3.1 -1.562 lum Lumican Chr2:3.1 -1.367 mmp13b Matrix Metallopeptidase 13b Chr20:8.3 -0.437 mmp2 Matrix Metallopeptidase 2 Chr7:15.9 -0.710 bmp2b bone morphogenetic protein 2b Chr11:7.0 0.589 bmp3 bone morphogenetic protein 3 Chr6:0.2 0.663 bmpr1aa bone morphogenetic protein receptor, type Iaa Chr18:0.7 0.524 dkk3a Dickkopf WNT Signaling Pathway Inhibitor 3a Chr7:3.1 -0.969 dlx3b Dickkopf WNT Signaling Pathway Inhibitor 3b Chr8:1.4 -1.382 fgf23 fibroblast growth factor 23 Chr2:36.0 2.363 fgf8b fibroblast growth factor 8b Chr14:7.3 3.022 fgfr3 fibroblast growth factor receptor 3 Chr3:2.0 0.489 fzd2 frizzled class receptor 2 Chr15:17.5 0.933 fzd4 frizzled class receptor 4 Chr20:12.3 0.968 fzd9 frizzled class receptor 9b Chr20:11.6 0.796 hoxa10a Homeobox A10a Chr10:1.1 0.668 tbx15 T-box transcription factor 15 Chr4:3.2 1.301 wnt3a Wnt Family Member 3a Chr1:16.4 1.394 DE-miRNA mir-122 miR-122-5p (osteoblastogenesis) Chr 5:9.1 4.225 mir-138a miR-138a-3p (osteoblastogenesis) Chr 8:5.5 -4.137 mir-144 miR-144-3p (osteoclastogenesis) Chr 17:13.6 -4.064 mir-145 miR-145-5p (chondro/osteogenic) Chr 9:7.5 -0.997 mir-146a miR-146a-5p (osteoclasto/ chondrogenesis) Chr 18:16.3 -1.374 mir-182a miR-182a-3p (chondro/osteogenic) Chr 2:5.0 1.495 mir-194a miR-194a-5p (osteoblastogenesis) Chr 17:14.6 4.157 mir-210 miR-210-5p (osteoblastogenesis) Chr 8:20.4 -1.432 mir-223 miR-223-3p (osteoclastogenesis) Chr 17:17.0 -1.327 mir-3381b miR-338-1b-5p (osteoclastogenesis) Chr 15:13.1 1.082 Chr.: Chromosome number and location in Mb (NCBI Genome Assembly GCA_919967415.2); FC: log2 fold changes. Table 2 Candidate genes and miRNA targets related to bone cell processes identified in the vertebral transcriptome of Solea senegalensis. Gene Coding gene description / miRNA target Chr:position (Mb) Osteoblastogenesis alpl alkaline phosphatase, tissue non-specific Chr19:1.0 col1a1b Collagen I alpha 1b Chr12:1.1 ctnnb1 Catenin Beta 1 Chr13:0.8 dlx3 Distal-Less Homeobox 3 Chr5:2.6 dlx5a Distal-Less Homeobox 15 Chr10:1.1 hoxa10a Homeobox A10a Chr10:1.1 msx1 Msh Homeobox 1 Chr19:25 msx2 Msh Homeobox 2 Chr13:9.8 runx2 runt-related transcription factor 2 Chr11:6.4 sparc Secreted Protein Acidic And Cysteine Rich (ostenectin) Chr9:1.8 spp1 Secreted Phosphoprotein 1 (osteopontin) Chr5:0.003 jun-D transcription factor jun-D Chr21:14.0 twist2 twist-related protein 2 Chr19:19.5 mir-146a smad2, smad3 Chr18:16.3 mir-155 bmpr2, runx2 Chr7:3.3 mir-20a-1 bmp2, bmpr2 Chr7:3.3 mir-23a-2 hoxa10, runx2 Chr15:9.2 mir-24-2 runx2 Chr15:9.2 mir-27a runx2 Chr1:17.3 mir-217a dkk1 Chr18:17.8 Osteocytic function-related (cellular bones) fgf23 fibroblast growth factor 23 Chr2:36.0 phex Phosphate Regulating Endopeptidase X-Linked Chr1:3.5 sost sclerostin Chr12:1.2 mir-218a-1 sost Chr9:23.2 Chondrogenesis acan aggrecan Chr8:9.8 bmp2b bone morphogenetic protein 2b Chr11:7.0 bmp4 bone morphogenetic protein 4 Chr3:14.8 bmp7 bone morphogenetic protein 7 Chr13:24.7 col10a1a Collagen X alpha 1a Chr3:14.7 col2a1 Collagen II alpha 1a Chr13:23.5 fgf10 fibroblast growth factor 10 Chr5:13.6 fgfr2 fibroblast growth factor receptor 2 Chr18: 11.4 igf1 insulin-like growth factor I Chr2:33.3 ihh indian hedgehog B protein-like Chr4:36.1 sox9 Transcription Factor sox-9 Chr12:15.7 tgfb3 transforming growth factor, beta 3 Chr3:18.3 mir140 hdac4 Chr8:0.7 mir-145 sox9 Chr9:7.5 mir-29b-1a hdca4 Chr13:20.3 mir-206-1 col2a1 Chr11:10.6 Osteoclastogenesis acp5a (trap) tartrate-resistant acid phosphatase Chr14: 6.2 m-csf (csf1b) macrophage colony-stimulating factor 1b Chr13:25.9 mmp13a Matrix Metallopeptidase 13 a (collagenase 3) chr17:11.3 myo10 unconventional myosin X Chr1:40.5 spi1 (PU.1) transcription factor PU.1 Chr8:16.3 tfgb1 transforming growth factor, beta 1a Chr20:21.0 tnfrsf11a (rank) tumor necrosis factor receptor superfamily member 11 A Chr1:26.5 tnfrsf11b (opg) TNF Receptor Superfamily Member 11b Chr16:15.7 tnfsf11 (rankl) tumor necrosis factor superfamily member 11 Chr4:1.3 mir-148a tnfsf11 Chr4:16.9 mir-155 m-csf, spi1 Chr7:3.3 mir-187 tnfsf11 Chr10:13.1 mir-21a mmp2, cfos Chr17:12.6 Chr.: Chromosome number and location in Mb (NCBI Genome Assembly GCA_919967415.2). C. Bouza et al. Genomics 116 (2024) 110802 8 candidates related to osteoblast differentiation and mineralization in mammals [99]. 3.5. miRNA transcriptomic data The effect of miRNAs on bone homeostasis and disease has been studied in humans and model vertebrates [39,41,42,100]. In farmed fish, conserved miRNAs have been involved on productive traits, including growth, stress, immune system, or reproduction, but scarcely in relation to bone development and skeletal anomalies [44,101–103]. We identified a total of 317 miRNA, and 201 hairpins for annotated miRNA genes, with one or both annotated strands (85 and 116, respectively; among the former, 41 from 3p strand and 44 from 5p) (Supplementary Table S5A). Some miRNAs were described and mapped for the first time in Senegalese sole, identifying many miRNA clusters and new duplicated paralogous members, thus improving the annotation of small non-coding RNAs related to fish vertebrae in this species [51,57]. The whole vertebral miRNAome showed differentiated profiles between postlarval and juvenile stages in the conditions studied (Supplementary Fig. S3A), as reported across different life-stages and tissues in other fish [43,46]. The most abundant mature miRNAs included members related to bone cells and homeostasis in model fish and vertebrates [38–42], some of them only detected in postlarvae (mir-10a-5p, mir182a-5p) or juveniles (mir-101a-3p, mir-126b-3p), but mostly shared between stages (Supplementary Table S5B). Several bone-related candidate miRNAs were identified in the Senegalese sole vertebrae (Tables 1 and 2; Supplementary Table S5C), such as mir-20a, associated with osteogenesis via targeting bmp2 and bmpr2 in teleost and mammals, supporting the role of miRNA mediated regulation on fish skeletogenesis [104]. Other key members and families were annotated, such as mir-26a involved in osteogenic differentiation, let-7 and miR-140, related to hypertrophic chondrocytes and osteogenesis, or mir-148a regulating osteoclastogenesis and intermuscular bone in fish [105]. The physical mapping of miRNA genes showed quite uneven chromosomal distribution (from 2 at Chr11 and 12, to 22 at Chr13; Supplementary Fig. S3B), and allowed exploring their location respect to the coding transcriptome of Senegalese sole [51]. Among miRNA genes, 55.7% were found at intergenic regions and 42.3% at intragenic positions, plus a few overlapping intra-intergenic locations (2%) (Supplementary Table S5A), within the ranges reported in fish and vertebrates [106]. Thirty-five clusters of miRNA genes were identified across chromosomes, except Chr3, Chr5, Chr12 and Chr16 (Supplementary Fig. S3B and Supplementary Table S5A), the two largest clusters being found at Chr7 and Chr9. Comparative mapping revealed high syntenic conservation of miRNA genes and clusters between Senegalese sole and O. latypes, useful for functional and comparative genomics of the miRNA repertoires in fish [47,107,108]. 3.6. Differential expression of miRNA between stages A total of 114 DE-miRNAs were detected when comparing postlarval and juvenile vertebral bone, nearly half upregulated in postlarvae (Fig. 3A; Supplementary Table S6A), which included a large fraction (40%) of bone cells-related candidates in vertebrates (Supplementary Table S5C). Clustered miRNA genes showed concordant expression patterns, either upregulated in postlarvae (such as 216b-1/216a-1/217a at Chr18) or juveniles (23a-2/27c/24-2 at Chr15). Among upregulated in postlarvae, mir-190b-5p, mir-217a-3p and mir-216a-1-3p, have been also reported enriched at early stages in other fish [46]. Other candidates, such as mir-122-5p, mir-130c-5p and mir-206-3p, were related to thermal plasticity and myogenesis in Senegalese sole, through different pathways [58], including mTOR signaling, associated with skeletal development [109]. Similarly, mir-192 and mir-181a, associated with thermal imprinting during early development in fish [43], were upregulated in postlarvae. Candidate miRNAs associated with osteoblastic differentiation and bone formation were differentially expressed between stages, either upregulated in postlarvae (such as mir-122, mir-194, mir-206-1) or juveniles (such as mir-138a and mir-210); the same was observed for miRNAs related to chondroor osteoclastogenic functions (Table 1), revealing complex miRNA regulation of these processes (Supplementary Table S5C). Other upregulated miRNAs in juveniles included mir-146a, mir-148a, and the cluster mir-144/mir-451, involved in osteoclastogenic regulation in mammals [39,110]. This observation matches with DEG from mRNA data in our study, suggesting more active osteoclastic and bone remodeling mechanisms in juveniles in the conditions studied. 3.7. miRNA target prediction and miRNA-mRNA interaction A total of 473 potential targets of 135 miRNAs were identified widely scattered across chromosomes (Supplementary Table S7A). These results provide a catalogue of potential miRNA-mRNA interactions for target genes involving MF related to signaling and catalytic activity, and BP to cellular and metabolic processes, signaling, localization and, to a lesser extent, growth, immune response, cell proliferation and biological adhesion (Fig. 3B). Suggestive mRNA and miRNA interactions included key coding candidates and miRNA genes related to bone growth and development, such as runx2 bearing targets for mir-30e-2-3p and mir-75p, previously reported to be involved in inhibiting osteogenic differentiation in mammals [41,111]. A subset of gene targets (53) for DE-miRNA between stages were identified (Supplementary Table S7B), 16 of them targeting DEGs (Fig. 3C; Supplementary Table S7C), including bone-related candidate miRNA (e.g., mir-101b, mir-148a, mir-206-1, mir-222a, mir-23a, mir26b) and coding genes, such as fgf23, tbx6, cellular retinoic acid binding protein 1 (crabp1) or nuclear receptor interacting protein 1b (nrip1b). Among them, mir-23a and its predicted target in this study, rgs5 (regulator of G-protein signaling 5), were candidates related to chondrogenesis activation and osteoblastogenesis regulation via runx2 signaling [112,113]. Some DE-miRNAs and target DEGs between stages showed negative correlation, suggesting interactions between bonerelated candidates; for instance, mir-10d and mir-let-7c were upregulated in postlarvae respect to their respective predicted targets, tbx6 associated with spinal development and fgf23, a regulator of bone mineralization [114,115]. Inhibition of let-7c has been related to osteoblastic promotion in mammals [116]. Conserved coding genes and miRNAs retrieved functional networks and pathways (Supplementary Fig. S4 and Supplementary Table S8), which lay a groundwork for further exploration of regulatory mechanisms involved in bone development and growth. In addition, transcriptomic markers previously associated with osteoblastogenesis, chondrogenesis and osteoclastogenesis processes in vertebrates (Tables 1 and 2; Supplementary Tables S4 and S5C), constitute a useful reference for functional studies of the anosteocytic bone in Senegalese sole and other teleosts [33,117]. Candidate genes related to osteocytic functions found in our study included sclerotin (sost), a key inhibitor of the WNT pathway and bone formation almost exclusively expressed by mammalian osteocytes. Conversely, a variety of non-osteocytic cells (such as osteoblasts, chordoblasts and chondroblasts) found in the acellular bone of medaka, support an osteocyte-independent bone remodeling mechanism in fish [117]. Similarly, conserved miRNA members associated with osteocytic functions were identified in the Senegalese sole vertebrae, such as mir-145, mir-181b and mir-218a, this latter inhibitor of sost [42], interesting for functional studies of fish bones, in combination with histopathology analysis across stages in this species [11], as reported in other fish [78,81,91,117]. 4. Conclusions This study reports for the first time a combined mRNA and miRNA transcriptomic data from vertebral bone of Senegalese sole, covering a C. Bouza et al. Genomics 116 (2024) 110802 9 wide range of expression profiles at postlarval and juvenile stages for different spinal phenotypes and early thermal conditions, influencing fish skeletal development. This broad transcriptomic repertoire adds new genomic resources and insights into the genetic pathways and regulatory mechanisms involved in the vertebral development and growth for this flatfish with anosteocytic bone, useful for future comparative studies of fish skeletogenesis. Candidate coding genes and miRNAs related to bone homeostasis and disease were identified, involving a balance of osteogenic, chondrogenic and osteoclastogenic processes, which showed highly conserved molecular mechanisms in vertebrates and teleosts. Our data provide a baseline of transcriptome markers for future studies to compare normal and specific types of anomalous spines, and to explore cell lineages of vertebral tissues in the Senegalese sole using histological and transcriptomic approaches. The transcriptome here represented constitutes a reference for further comparisons on the functional mechanisms underlying spinal bone homeostasis and pathology in the Senegalese sole within fish and vertebrates. Supplementary data to this article can be found online at https://doi. org/10.1016/j.ygeno.2024.110802. CRediT authorship contribution statement Carmen Bouza: Writing – review & editing, Writing – original draft, Supervision, Resources, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. Ana P. Losada: Writing – review & editing, Writing – original draft, Methodology, Investigation, Formal analysis. Carlos Fern´ andez: Writing – review & editing, Writing – original draft, Investigation, Formal analysis. Jos´ e A. ´ Alvarez-Dios: Writing – review & editing, Investigation, Formal analysis. Ana Manuela de Azevedo: Writing – review & editing, Investigation, Formal analysis. Andr´ es Barreiro: Writing – review & editing, Resources, Methodology, Investigation. Dami´ an Costas: Writing – review & editing, Resources, Methodology. María Isabel Quiroga: Writing – review & editing, Resources, Methodology. Paulino Martínez: Writing – review & editing, Resources, Methodology, Investigation. Sonia V´ azquez: Writing – review & editing, Supervision, Resources, Methodology, Investigation, Formal analysis, Conceptualization. Declaration of competing interest The authors declare no conflict of interest. Data availability All RNA-seq data were submitted to NCBI-SRA database (BioProject ID: PRJNA820527, Biosample accessions: SAMN31285849SAMN31285789). Acknowledgments This study was supported by the Spanish Ministry of Science, Innovation and Universities, and European Regional Development Fund (RTI2018-097110-B-C21). AMdA has awarded a postdoctoral contract from Consellería de Cultura, Educaci´ on e Ordenaci´ on Universitaria, Xunta de Galicia, Spain. We thank Stolt Sea Farm S.A. for providing the egg samples, and the personal involved in fish culture and sampling at the Station of Marine Science of Toralla (ECIMAT, University of Vigo). We acknowledge the bioinformatic support of Andr´ es Hortas, ´ Oscar Aramburu and Galicia Supercomputing Center (CESGA). Thanks to Lucía ´ Insua and Sonia G´ omez for technical and administrative support. References [1] S. Morais, C. Arag˜ ao, E. Cabrita, L.E.C. Conceiç˜ ao, M. Constenla, B. Costas, J. Dias, N. Duncan, S. Engrola, A. Estevez, E. Gisbert, E. Ma˜ nan´ os, L.M.P. Valente, M. Yúfera, M.T. Dinis, New developments and biological insights into the farming of Solea senegalensis reinforcing its aquaculture potential, Rev. Aquac. 8 (2016) 227–263, https://doi.org/10.1111/raq.12091. [2] APROMAR, Aquaculture in Spain, The Spanish Marine Aquaculture Producers Association, 2022. https://apromar.es/wp-content/uploads/2022/10/Aquacult ure-in-Spain-2022_APROMAR.pdf. [3] I. Fern´ andez, M.S. Pimentel, J.B. Ortiz-Delgado, F. Hontoria, C. Sarasquete, A. Est´ evez, J.L. Zambonino-Infante, E. Gisbert, Effect of dietary vitamin a on Senegalese sole (Solea senegalensis) skeletogenesis and larval quality, Aquaculture 295 (2009) 250–265, https://doi.org/10.1016/j.aquaculture.2009.06.046. [4] P.J. Gavaia, S. Domingues, S. Engrola, P. Drake, C. Sarasquete, M.T. Dinis, M. L. Cancela, Comparing skeletal development of wild and hatchery-reared Senegalese sole (Solea senegalensis, Kaup 1858): evaluation in larval and postlarval stages, Aquac. Res. 40 (2009) 1585–1593, https://doi.org/10.1111/ j.1365-2109.2009.02258.X. [5] A.M. de Azevedo, A.P. Losada, I. Ferreiro, A. Riaza, V. Losada, T. Russo, C. Boglione, S. V´ azquez, M.I. Quiroga, Skeletal anomalies in Senegalese sole (Solea senegalensis, Kaup) fed with different commercial enriched artemia: a study in postlarvae and juveniles, Animals 11 (2021) 22, https://doi.org/10.3390/ ani11010022. [6] C. Boglione, P. Gavaia, G. Koumoundouros, E. Gisbert, M. Moren, S. Fontagn´ e, P. E. Witten, Skeletal anomalies in reared European fish larvae and juveniles. Part 1: normal and anomalous skeletogenic processes, Rev. Aquac. 5 (2013) S99–S120, https://doi.org/10.1111/raq.12015. [7] C. Boglione, P. Gavaia, G. Koumoundouros, E. Gisbert, M. Moren, S. Fontagn´ e, P. E. Witten, Skeletal anomalies in reared European fish larvae and juveniles. Part 2: main typologies, occurrences and causative factors, Rev. Aquac. 5 (2013) S121–S167, https://doi.org/10.1111/raq.12016. [8] A. Lorenzo-Felipe, H.S. Shin, S. Le´ on-Bernabeu, C. P´ erez-García, M.J. Zamorano, J. P´ erez-S´ anchez, J.M. Afonso-L´ opez, The effect of the deformity genetic background of the breeders on the spawning quality of gilthead seabream (Sparus aurata L.), Front. Mar. Sci. 8 (2021) 656901, https://doi.org/10.3389/ fmars.2021.656901. [9] A.P. Losada, A.M. de Azevedo, A. Barreiro, J.D. Barreiro, I. Ferreiro, A. Riaza, M. I. Quiroga, S. V´ azquez, Skeletal malformations in Senegalese sole (Solea senegalensis Kaup, 1858): gross morphology and radiographic correlation, J. Appl. Ichthyol. 30 (2014) 804–808, https://doi.org/10.1111/jai.12524. [10] A.M. de Azevedo, A.P. Losada, I. Ferreiro, A. Riaza, S. V´ azquez, M.I. Quiroga, New insight on vertebral anomalies in cultured Senegalese sole (Solea senegalensis, Kaup) at early stages of development, J. Fish Dis. 40 (2017) 987–1000, https:// doi.org/10.1111/jfd.12575. [11] A.M. de Azevedo, A.P. Losada, A. Barreiro, S. V´ azquez, M.I. Quiroga, Skeletal anomalies in Senegalese sole (Solea senegalensis), an anosteocytic boned flatfish species, Vet. Pathol. 56 (2019) 307–316, https://doi.org/10.1177/ 0300985818800027. [12] P.E. Witten, L. Gil-Martens, B.K. Hall, A. Huysseune, A. Obach, Compressed vertebrae in Atlantic salmon Salmo salar: evidence for metaplastic chondrogenesis as a skeletogenic response late in ontogeny, Dis. Aquat. Org. 64 (2005) 237–246, https://doi.org/10.3354/dao064237. [13] A.M. de Azevedo, A.P. Losada, A. Barreiro, J.D. Barreiro, I. Ferreiro, A. Riaza, S. V´ azquez, M.I. Quiroga, Skeletal anomalies in reared Senegalese sole Solea senegalensis juveniles: a radiographic approach, Dis. Aquat. Org. 124 (2017) 117–129, https://doi.org/10.3354/dao03110. [14] A. Printzi, D. Mazurais, P.E. Witten, L. Madec, A.A. Gonzalez, X. Mialhe, J. L. Zambonino-Infante, G. Koumoundouros, Juvenile zebrafish (Danio rerio) are able to recover from lordosis, Sci. Rep. 12 (2022) 21533, https://doi.org/ 10.1038/s41598-022-26112-2. [15] J. Castro, A. Pino-Querido, M. Hermida, D. Chavarrias, R. Romero, L.A. GarcíaCort´ es, M.A. Toro, P. Martínez, Heritability of skeleton abnormalities (lordosis, lack of operculum) in gilthead seabream (Sparus aurata) supported by microsatellite family data, Aquaculture 279 (2008) 18–22, https://doi.org/ 10.1016/j.aquaculture.2008.04.023. [16] D. Negrín-B´ aez, A. Navarro, I. Lee-Montero, M. Soula, J.M. Afonso, M. J. Zamorano, Inheritance of skeletal deformities in gilthead seabream (Sparus aur–ta) - lack of operculum, lordosis, vertebral fusion and LSK complex, J. Anim. Sci. 93 (2015) 53–61, https://doi.org/10.2527/jas.2014-7968. [17] D. Negrín-B´ aez, D. Negrín-B´ aez, S.T. Rodríguez-Ramilo, J.M. Afonso, M. J. Zamorano, Identification of quantitative trait loci associated with the skeletal deformity LSK complex in Gilthead seabream (Sparus aurata L.), Mar. Biotechnol. 18 (2016) 98–106, https://doi.org/10.1007/s10126-015-9671-7. [18] F. Lavajoo, M. Perell´ o-Amor´ os, E.J. V´ elez, A. S´ anchez-Moya, S. Balbuena-Pecino, N. RieraHeredia, J. Fern´ andez-Borr` as, J. Blasco, I. Navarro, E. Capilla, J. Guti´ errez, Regulatory mechanisms involved in muscle and bone remodeling during refeeding in gilthead sea bream, Sci. Rep. 10 (2020) 184, https://doi.org/ 10.1038/s41598-019-57013-6. [19] E. Ytteborg, J. Torgersen, G. Baeverfjord, H. Takle, Morphological and molecular characterization of developing vertebral fusions using a teleost model, BMC Physiol. 10 (2010) 13, https://doi.org/10.1186/1472-6793-10-13. [20] T.W. Fraser, T. Hansen, M.S. Fleming, P.G. Fjelldal, The prevalence of vertebral deformities is increased with higher egg incubation temperatures and triploidy in Atlantic salmon Salmo salar L, J. Fish Dis. 38 (2015) 75–89, https://doi.org/ 10.1111/jfd.12206. [21] J.S.F. Lee, M.A. Cook, J.A. Luckenbach, B.A. Berejikian, C.A. Simchick, S. M. Oden, F.W. Goetz, Investigation of long-term effects of larval rearing temperature on growth, deformities, flesh quality, and phenotypic sex of cultured C. Bouza et al.