Received date: February 05, 2015; Accepted date: April 22, 2015; Published date: April 24, 2015
Citation: Le Luyer J, Deschamps MH, Proulx E, Poirier Stewart N, Droit A, et al. (2015) RNA-Seq Transcriptome Analysis of Pronounced Biconcave Vertebrae: A Common Abnormality in Rainbow Trout (Oncorhynchus mykiss, Walbaum) Fed a Low-Phosphorus Diet. Next Generat Sequenc & Applic 2:112. doi:10.4172/2469-9853.1000112
Copyright: © 2015 Le Luyer, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Next Generation Sequencing & Applications
The prevalence of bone deformities, particularly linked with mineral deficiency, is an important issue for fish production. Juvenile triploid rainbow trout (Oncorhynchus mykiss) were fed a low-phosphorus (P) diet for 27 weeks (60 to 630 g body mass). At study termination, 24.9% of the fish fed the low-P diet displayed homogeneous biconcave vertebrae (deformed vertebrae phenotype), while 5.5% displayed normal vertebral phenotypes for the entire experiment. The aim of our study was to characterize the deformed phenotype and identify the putative genes involved in the appearance of P deficiency-induced deformities. Both P status and biomechanical measurements showed that deformed vertebrae were significantly less mineralized (55.0 ± 0.4 and 59.4 ± 0.5,% ash DM, for deformed and normal vertebrae, respectively) resulting in a lower stiffness (80.3 ± 9.0 and 140.2 ± 6.3 N/mm, for deformed and normal phenotypes, respectively). The bone profiles based on μCT observations showed no difference in the osteoclastic resorption while no difference in matrix production was observed between deformed (total bone area 5442.0 ± 110.1 μm2) and normal vertebrae (total bone area 5931.2 ± 249.8 μm2) in this study. Consequently, the lower P content rather results from a reduced degree of mineralization in the deformed phenotype. Finally, we quantified differential gene expression between deformed vertebrae (pronounced biconcave) and normal phenotype by employing deep RNA-sequencing and mapping against a reference bone transcriptome for rainbow trout. In total, 1289 genes were differentially expressed. Among them, in deformed fish we observed that BGLAP, MGP and NOG, an inhibitor of BMP signalling pathway, were up-regulated while COL11a1 was down-regulated. These genes are central actors involved in the reduced degree of mineralization triggering vertebral deformities. These results will further the understanding of P deficiency-induced deformities; hence providing new tools for improved P management in production settings
RNA-sequencing; Deformities; Rainbow trout; Phosphorus
BGLAP: Osteocalcin; BMPR2: Bone Morphogenetic Protein Receptor, type II; CEBPB: CCAAT/Enhancer Binding Protein (C/EBP) Β; COL11a1: Collagen type XI alpha 1; COL12a1: Collagen type X2 alpha 1, CREBBP; CREB Binding Protein; DCHS1: Dachsous Cadherin-Related 1; FAT4: FAT Atypical Cadherin 4; FBN2: Fibrillin; FGFR4: Fibroblast Growth Factor Receptor 4; FLNA: Filamin A, alpha; GCG2: Glucagon-2; GPC1: Glypican 1; HIVEP3: Human Immunodeficient Virus Type 1 Enhancer Binding Protein 3; HSPG2: Heparan Sulfate Proteoglycan 2; JUNB: Jun B proto-oncogene; JUND: Jun D proto-oncogene; LRP6: Low density lipoprotein receptor-Related Protein 6; M3K1: Mitogen-Activated Protein Kinase Kinase Kinase 1; MEF2C: Myocyte-Specific Enhancer Factor 2C; MGP: Matrix Gla Protein; MMP10: Matrix Metallopeptidase 10; MYCB2: Myc Binding Protein 2; MYCPP: C-myc Promoter-Binding Protein; MYO10: Unconventional Myosin-X; NF1: Neurofibromin 1; NFIX: Nuclear Factor 1/X (CCAAT-binding transcription factor); NIK: Mitogen- Activated Protein Kinase Kinase Kinase 4; NKX3.2: Homeobox Protein Nkx3.2; NOG: Noggin; OPG: Osteoprotegerin; SFRP1: Secreted Frizzled-Related Protein 1; SHH: Sonic Hedgehog; SOS1: Son of Sevenless Homolog 1; SOX9: Sex Determining Region Y-Box 9; SPARC: Osteonectin; SPP1: Osteopontin; SRCAP: nf2-Related CREBBP Activator Protein; TNR5: CD40 molecule TNF Receptor Family Member 5; TRAP: Tartrate Resistant Alkaline Phosphatase; TMEFF1: Tomoregulin-1; TWIST2: Twist Family BHLH Transcription Factor 2; USP9X: Probable Ubiquitin Carboxyl-Terminal Hydrolase FAF-X; VDR: Vitamin D Receptor; VGFR1: Vascular Endothelial Growth Factor Receptor 1
Extensive work regarding deformity development in teleost fish have revealed that biotic and abiotic factors affect this condition including genetics, ploidy as well as nutrition [1-3]. Among the main nutrients affecting skeletal development in fish, phosphorus (P) particularly has drawn particular attention [4-6]. Bone tissue mineral (hydroxyapatite) constitutes the main reservoir of P, comprising up to 85% of the total body content [7,8]. To date, many studies in salmonids have attempted to identify and characterize the response of the fish skeleton to low-P diets. Under extreme conditions such as reproduction, migration or starvation, teleosts mobilize P from scales and bone tissues [9-13]. During a prolonged P deficiency, fish may also display lower growth performance , lower mineral status of bone, scales and flesh [3,15,16] and, in some severe cases, increased occurrence of specific phenotypes of vertebral deformities [17-20].
The study of vertebral abnormalities is particularly difficult in regard to the complexity and the role of the axial skeleton. In trout, this system is comprised of 64 vertebrae, composed of several tissues and cell types, interconnected by intervertebral tissue and undergoing compression forces during swimming. As in other vertebrates, under non-pathological conditions, a tight balance in bone tissue is maintained between bone apposition by osteoblasts and bone resorption, resulting in a continuous turnover. In trout and in other species that belong to basal teleostean lineages, osteoblasts are trapped in the matrix and become osteocytes, a condition that defines cellular bone . The osteocytes also play a role in bone apposition and bone resorption around their lacunae . In vertebrae, chondrocytes are also present at the base of neural and hemal arches (basalia). Their progressive endochondral mineralization enables the complete fusion of the arches with the vertebral body during juvenile growth, at least in the caudal portion of the axial skeleton . Thus, the occurrence of vertebral anomalies reflects a disturbance of the homeostatic mechanisms and well-coordinated actions of several cell types.
The molecular pathways involved in P-deficiency-induced deformities are far from being well understood. In Atlantic salmon (Salmo salar) presenting vertebral fusions [23,24], more than 20 actors including transcription factors and genes encoding molecules involved in bone cell activation and differentiation, and extracellular matrix (ECM) formation were detected. Similarly, a large-scale study on the impact of photoperiod on vertebral bone formation was conducted in Atlantic salmon  and highlighted the structural role of COL11a1 in the occurrence of vertebral deformities. Recently, characterization of vertebrae transcriptomes in sea bream (Sparus auratus)  and rainbow trout  suggested that homologs of key genes involved in bone metabolism are conserved in vertebrates. As bone metabolism in fish is complex, transcriptomic-wide analyses may be the best approach to identify the principal affected pathways during P deficiency.
A recent study  showed that the most common deformities observed in P-deficient all-female triploid rainbow trout are pronounced biconcave and homogeneously-compressed vertebrae, biconcave vertebrae being considered an early stage of vertebral compression. Most of the fish undergoing prolonged P deficiency develop vertebral abnormalities, however, some individuals are able to maintain a normal phenotype. Given this phenotypic plasticity, it is of interest to better understand the cellular mechanisms underlying the appearance of vertebral abnormalities.
The aim of this study was to identify molecular changes that are involved in the appearance of pronounced biconcave vertebrae (Figure 1) when trout are fed a low-P diet. We compared gene expression in normal and biconcave vertebrae from P-deficient all-female triploid rainbow trout using next-generation sequencing (NGS) on the Illumina HiSeq2000 platform. To support the transcriptomic results, in each individual we characterized the P status and the quality of bone tissue based on chemical, histological and biomechanical determinations. We discuss our findings with regard to transcription regulation and the formation, mineralization and resorption of bone tissue.
Figure 1: Illustrations (left) and X-rays (right) of the reported vertebral phenotypes in this study (normal and pronounced biconcave, respectively) and associated with mineral deficiency in teleost. Illustrations were modified after Witten et al.  For further details regarding the phenotypes see Poirier Stewart et al. .
All-female triploid rainbow trout (N= 1,680, initial mass 60.8 ± 1.6 g) were transferred from St-Alexis-des-Monts Inc., Canada to experimental facilities at the Laboratoire de recherche des sciences aquatiques (LARSA), Université Laval in Québec City (Canada). An acclimation period of two weeks took place, during which fish were fed a commercial feed (Corey Optimum 3 mm) in accordance with manufacturer’s tables. Thereafter, fish (N=840) were fed with practical P-deficient diet (available P: 0.29%) using either apparent satiation or pair-fed feeding regimes. The P-deficient diet was found to induce severe deficiency in rainbow trout as already reported [16,19,20].
The fish were reared for 27 weeks in six circular tanks (2000 L) in a recirculating aquaculture system (6 replicates, n=140/tank, 8.5 to 52.7 kg m-3, 12 ± 0.3°C), with near oxygen saturation and at constant, long (LD, 18:6) photoperiod. During experimental manipulations, fish were anaesthetised in a MS-222 bath (75 mg L-1; Syndel International Inc., Vancouver, BC, Canada). Experiments took place in compliance to the guidelines of the Canadian Council on Animal Care (1984) and supervised by the Animal Protection Committee of Université Laval.
Deformities assessment and tissue sampling
At week 5, all trout were PIT-tagged with a microchip enabling individual survey (weight, length and condition factor: K=105 × g/ mm3) for each sample. X-rays permitted individual monitoring and was performed at week 5, 15 and 24. X-rays were performed on the trunk-caudal region of all fish at 60 kV (2s, 15 mA, distance of 40 cm). Monitoring of vertebral abnormalities was assessed directly on developed films and abnormalities along the vertebral axis were sorted according to the classification proposed for Atlantic salmon . Fish were associated to a phenotype when all caudal vertebrae (V31-V44) showed a homogenous pattern of abnormalities (see details in ). At the termination of the study (week 27), P-deficient fish displaying normal (n=3 fish) and pronounced biconcave (n=3 fish) vertebrae were randomly sampled. This abnormality type was chosen as it was by far the most represented type observed in the current study . Moreover, pronounced biconcave vertebrae are related to early stages of mineral deficiency and might be a transition stage to vertebral compression, which are of major concern for fish growth . For each fish, blood was collected from the caudal aorta, centrifuged, and the plasma stored at -20 C. Scales were scraped from tail to head and were collected in a 70° ethanol solution. Caudal vertebrae were removed for bone histology (V32-33 and V38-39), bone genomics (V35-37), bone structure (V40- 42), biomechanics and bone chemistry (V31-32 and V43-44). All of these analyses were conducted on vertebrae from the same individuals.
Inorganic phosphorus in plasma
Colorimetric assays (phosphomolybdate; K003912, JAS diagnostics, Inc., USA) were performed to measure inorganic phosphates (Pi=PO43-) in individual samples of plasma (150 μl) using an automatic spectrophotometer (Trilogy Multi-purpose diagnostic analyzer system, Drew Scientific, USA).
To quantify the mineralization of basalias, two attached vertebrae (V38-39) were fixed in 4% formaldehyde and decalcified in 25% decal solution for 24 hrs. Vertebrae were then dehydrated in graded series of ethanol prior to being embedded in paraffin. Up to 15 adjacent longitudinal cross-sections (5 m) from either side of the mid-region of the vertebrae (i.e., where the notochord canal is present at the narrowest) were obtained using a microtome (Microm HM 330, Heidelberg, Germany). The sections were stained using Sirius red techniques . Two other attached vertebrae V33-34 (fixed in cold acetone) were embedded in glycol methacrylate. One longitudinal section was obtained by sanding one side of the block to the mid-region of the vertebrae and by cutting the other side of the block to a thickness of 1.5 mm using a band saw and a cutting guide. The sections were polished before being stained using Von Kossa techniques .
Morphometrical results at the vertebral body endplates (VBE) were estimated from the Von Kossa sections. For each individual, tiff images from the longitudinal section (4x magnification) were obtained using a digital camera (EXI Blue Q35618, Q Imaging Inc., BC, Canada) mounted on a microscope (Nikon eclipse E600) and captured with Q-capture Pro 7 software (Q Imaging Inc., BC, Canada). Once transferred to ImageJ v1.40 (U. S. National Institutes of Health, USA, http://rsb.info.nih.gov./ij/) software (with Java, 32 bits), the freehand selection tool was used to manually delineate (graphic tablet Intruos3, model PTZ-631W, Wacom Co. Ltd, China) areas and calculate a number of pixels for the mineralized trabecular area (tb) and total bone area, comprising tb and the un-mineralized osteoid at the apposition front (os), for each vertebral body endplate (VBE) (n=8/individual). The mean proportion of os (%) at the VBE was calculated as follows: os/(os + tb) x 100. The length (LVBE) and thickness (TVBE) of the vertebral body endplates were also measured. LVBE included the total distance from the notochord canal to the apposition front. TVBE was defined as the perpendicular distance from the point of attachment of the notochord tissue to the inner limit of the compact bone layer. A schematic representation including the different vertebral body parts mentioned in this section is proposed in Figure 2.
Figure 2: Scheme of a longitudinal section of the mid-region of two adjacent vertebrae including basalias and notochordal tissues. With b = basalia; LVBE = length of the vertebral body endplate; N = notochordal tissue; os = osteoid; T = trabecular bone; V = vacuolated tissue; VBE = vertebral body endplates; TVBE = thickness of the vertebral body endplate.
The relative mineralization of the basalias was estimated from digital pictures of Sirius red sections at 200x magnification (N=5 images/basalia/individual section) using a digital camera (Infinity2- 3C; Lumenera 184.108.40.206) mounted on a binocular microscope (Olympus BX53) and captured with Image-Pro Plus software (version 220.127.116.111). The basalias were cut and cell nuclei removed using the selection tools of GIMP v2.8.10 software (ww.gimp.org). The number of pixels for two RGB colours (unmineralized collagen=red (222, 333, 136); mineralized collagen=grey (142, 59, 142) were analysed using ImageJ and colour deconvolution Plug-In as published elsewhere . The proportion of mineralized and unmineralized areas of the basalias were calculated as the number of pixels for a specific colour divided by the total number of all defined pixels × 100.
Bone density and structure of vertebrae V40-42 were estimated using the approach of the modelling of bone area profiles according to previously-published methods  and presented in Figure 3. Briefly, total bone area (Tt-B.Ar.) and mean bone area according to the relative distance from the notochordal canal were measured using Bone profiler software  and modeled using 12 parameters in four predefined regions of the vertebral body (notochord, transition, middle and periphery). Measurements were performed on 125 μm-thick transverse sections of the mid-region of the vertebrae, in which the notochord canal is the narrowest and bone tissue the largest. Virtual cross-sections were obtained from μCT images (National Museum of Natural History, France).
Figure 3: Vertebral bone profile for normal and deformed groups. Bone area (expressed as a ratio from 0 to 1) as a function of the relative distance to the notochord. nt = notochord area; tr = transition area; mi = middle area; pe = periphery area. Asterisk (*) shows a slight significant difference for the lost of bone in the middle area between groups (ANOVA P < 0.05; n=3 fish/ condition).
Four vertebrae (V31, 32, 43 and 44) were used for mechanical testing prior mineral content analyses. Neural and haemal arches were removed, and the amphicoelus centra were compressed along the cranial–caudal axis at a steady rate of deformation (0.1 mm sec–1) using a texture analyzer (TA-XT2 Texture Analyzer, Stable Micro Systems, UK) and the Texture Expert Exponent 32 acquisition software (TE32 version 18.104.22.168, Texture Technologies, USA). For each vertebra, load (N) and deformation (mm) data were continuously recorded and the stiffness (g.mm–1) was calculated according previously-published methods .
The mineral content was determined for scales and vertebrae (V31, 32, 43 and 44). The samples were dehydrated in a graded series of ethanol (70°, 90°, 100°; 24 h/bath), and lipids were extracted in acetone (two baths of 24 h), then in trichloroethylene (two baths of 24 h). Diets and tissues were analyzed for dry matter (drying in a vacuum oven for 18 h at 105°C) and ash (incineration in a muffle furnace for 18 h at 550°C), weighed to the nearest 0.1 mg according to AOAC guidelines (1990). P content was determined by ion chromatography (ICS-3000, Dionex Corporation, Sunnyvale, CA) following ash digestion in nitric acid solution and filtering (Whatman paper No. 1, rinsed 3 times). Separation was carried out on an IonPac AS 17 column (Dionex) with an IonPac AG 17 guard column. For sample loading, a 50 μl sample loop (PEEK) was used. Multi-step gradient concentrations ranging from 9 to 30 mM allowed for the complete separation of phosphorus anions in aqueous media in 16 min at a constant flow rate of 1.0 ml.min−1.
Three connected vertebrae (V35-37, comprised of ligaments and intervertebral tissue) were rapidly dissected, arches removed and the remaining vertebral body quickly cleaned of flesh by rinsing and brushing with PBS (1X) to remove all muscle tissue. The vertebrae from each individual was then immersed in liquid nitrogen, crushed into a fine powder, and total RNA extracted using TRIZOL Reagent (Life Technologies) at a ratio of 1 ml per 100 mg tissue according to manufacturer’s recommendations. Total RNA was subjected to DNAse treatment on a column (ArcturusPicoPure RNA isolation kit, Life Technologies), and integrity and purity assessed using a Bioanalyzer 2100 (Agilent Technologies). RNA was also quantified using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc.). For each sample, 1 μg of total RNA was used to create cDNA libraries with TrueSeq Sample Prep kit v2 (Illumina, San Diego, CA, USA) following manufacturer’s instructions. Resulting libraries were quantified using a Bioanalyzer 2100 (Agilent Technologies). Samples were multiplexed (6 samples per lane) and sequenced at McGill genomic platform with HiSeq2000 sequencer and 100 paired-end (PE) technology.
Differentially expression study and annotation
Reads from HiSeq2000 Illumina were first processed with Trimmomatic v0.30  to remove low quality (trailing: 20, lowest quality: 30) and short reads (<60 bp). Trimming with Trimmomatic also included removal of Illumina adapters together with the most common contamination vectors from UniVec database (https://www.ncbi.nlm.nih.gov/tools/vecscreen/univec/). Read abundance and differentially-expressed genes were estimated using Trinity pipeline published recently for non-model species  by mapping reads on a reference transcriptome for rainbow trout bone tissue . Briefly, gene abundance estimations have been assessed individually using RSEM tools  and bowtie v1.0.0 . The differentially-expressed gene assessment based on RPKM values, was conducted using both EdgeR v3.8.5  and DESeq v1.18.0  packages with correction of false discovery rate (FDR) of 5%. Only common genes detected as differentially expressed with both analyses have been used for further analyses and functional annotation. Blast-x searches (cut-off e-value<10-8) were performed on Uniprot-Swissprot database. Gene ontology (GO) association was realised using GOminer software . To estimate the proportion of long non-coding RNAs (lncRNAs), remaining sequences non-matching against Uniprot database were analyzed using PLEK software  (length threshold >200 bp). Visualization of biological functions was conducted using IPA core analysis (Ingenuity Systems, www.ingenuity.com). The Ingenuity platform allows novel protein sequence association with networks, canonical pathways and biological functions relative to human and mouse systems. The genes specifically related to bone metabolism were manually searched by compiling relevant metabolic pathways (Table 1). The analysis pipeline for gene expression quantification is described in Figure 4.
|Gene header||UniprotAcc.||Gene name||Log2FC||FC|
|Down-regulated genes in deformed vertebra|
|comp149051_c0||A4IFW2||Protein tyrosine phosphatase receptor type F||1.64||3.12|
|comp138459_c1||Q64487||Protein tyrosine phosphatase receptor type D||1.62||3.08|
|comp144178_c2||Q6ZRS2||Snf2-related CREBBP activator protein||1.49||2.8|
|comp149014_c5||Q8BNA6||FAT atypicalcadherin 3||1.37||2.59|
|comp144153_c0||O95644||Nuclear factor of activated T-cells, calcineurin-dependent||1.31||2.49|
|comp133481_c0||Q9Y3S1||Serine/threonine-protein kinase WNK2||1.31||2.49|
|comp148107_c0||Q8WYB5||K (Lysine) acetyltranferase 6B||1.24||2.35|
|comp143062_c0||Q13332||Protein tyrosine phosphatase receptor type S||1.19||2.28|
|comp126856_c0||P17948||Vascular endothelial growth factor receptor 1||1.17||2.24|
|comp134803_c0||Q14517||FAT atypicalcadherin 1||1.14||2.21|
|comp143303_c7||Q8VI24||SATB homeobox 2||1.12||2.18|
|comp145538_c0||Q99715||Collagen type X2 alpha 1||1.12||2.18|
|comp146160_c1||Q9Z180||SET bindingprotein 1||1.1||2.14|
|comp149571_c0||Q9P2D1||Chromodomain helicase DNA binding protein||1.09||2.13|
|comp140581_c1||Q92794||K (Lysine) acetyltranferase 6A||1.05||2.07|
|comp128868_c0||O75581||Low density lipoprotein receptor-related protein 6||1.04||2.05|
|comp142199_c3||O70496||Chloridechannel voltage-sesitive 7||1.01||2.01|
|comp147156_c2||O89026||Roundabout axon guidance receptor homolog 1||0.97||1.96|
|comp135768_c0||P27658||Collagen type VIII Alpha 1||0.97||1.97|
|comp149750_c0||Q6V0I7||FAT atypicalcadherin 4||0.97||1.96|
|comp134251_c0||Q8BTM8||Filamin A, alpha||0.97||1.96|
|comp144389_c0||A7XYQ1||Jackson circlerprotein 1||0.94||1.92|
|comp67973_c1||Q6NSM8||SIK family kinase 3||0.93||1.91|
|comp135751_c0||Q9UBS9||Sun domain-containing ossification factor||0.93||1.9|
|comp135751_c0||Q9QYP0||Multiple EGF-like-domain 8||0.91||1.88|
|comp145151_c0||P26012||Integrin β 8||0.9||1.87|
|comp136915_c0||Q13233||Mitogen-activated protein kinase kinasekinase 1||0.89||1.85|
|comp149523_c0||O88491||Nuclear receptor binding SET domain protein 1||0.87||1.83|
|comp149485_c0||Q9UHF7||Trichorhinophalangeal syndrome I||0.87||1.83|
|comp149389_c0||O35607||Bone morphogenetic protein receptor, type II||0.85||1.81|
|comp146007_c0||Q14938||Nuclear factor 1/X (CCAAT-binding transcription factor)||0.85||1.8|
|comp146819_c2||Q90413||Fibroblastgrowth factor receptor 4||0.81||1.75|
|comp146502_c4||P20909||Collagen type XI alpha 1||0.8||1.74|
|comp149588_c0||A2A884||Human immunodeficient virus type 1 enhancer binding protein 3||0.79||1.73|
|comp147279_c4||P98160||Heparan sulfate proteoglycan 2||0.78||1.72|
|comp144268_c0||Q09472||E1A bindingprotein P300||0.78||1.71|
|comp148958_c0||Q9ULT8||HECT domaincontaining E3 ubiquitinprotein ligase 1||0.78||1.71|
|comp149573_c0||P10587||Myosin heavy chain 11 smooth muscle||0.77||1.71|
|comp149223_c0||P79114||Unconventional myosin-X, role n podosom belt osteoclast||0.77||1.7|
|comp138258_c6||Q07889||Son of sevenlesshomolog 1||0.69||1.61|
|comp147381_c0||Q91YM2||Rho GTPaseactivatingprotein 35||0.68||1.6|
|Up-regulated genes in deformed vertebra|
|comp142082_c0||P57102||Heart and neural crest derivatives expressed 2||-1.84||3.59|
|comp148652_c4||Q9JK00||Sodium channel voltage gated type III β subunit||-1.63||3.09|
|comp149418_c1||P05549||Activating Enhancer-Binding Protein 2-Alpha||-1.51||2.84|
|comp142563_c0||P09238||Matrix metallopeptidase 10||-1.46||2.76|
|comp88970_c0||P79703||Jun B proto-oncogene||-1.42||2.67|
|comp144447_c2||Q05826||CCAAT/Enhancerbindingprotein (C/EBP) Β||-1.27||2.42|
|comp144817_c0||O75911||Dehydrogenase/Reductase (SDR family member 3)||-0.93||1.9|
|comp136344_c0||P52909||Jun D proto-oncogene||-0.85||1.8|
|comp148539_c0||P70061||NK3 homeobox 2||-0.85||1.81|
|comp143712_c0||P25942||CD40 molecule TNF receptor family member 5||-0.78||1.72|
|comp98891_c0||O55052||Macrophage migration inhibitory factor||-0.7||1.62|
|comp148620_c0||P18203||FK506 binding protein 1A, 12kDa||-0.7||1.63|
Table 1: Selection of differentially-expressed genes putatively related to bone metabolism. The homologies were revealed with blast-x search against Uniprot-Swissprot database (cut-off evalue<1e-8). Differentially expressed genes were conserved if they display a unique and FDR<5%.
We used RT-qPCR to assess differential gene expression for keys genes involved in each cellular types: MEF2C, TRAP and BGLAP, for chondrocytes, osteoclasts and osteoblasts, respectively. For each sample, total RNA was converted into complementary DNA (cDNA) using Omniscript reverse transcriptase kit (Qiagen). Specific primers (Table 2) were designed using Primer Quest from DNA Technologies (IDT, Coralville, IA) for Real-Time quantitative PCR (480 Roche lightcycler 2.0, Roche Diagnostics, Canada) and SYBR green I staining . Quality and specificity of amplicons were visually confirmed on 2% gel electrophoresis and by DNA sequencing. A normalization factor for each sample was assessed using GeNorm software  based on two housekeeping gene expression levels, β-actin and Heat Shock Protein 90B (HSP90B). Expression of genes assessed by RT-PCR is presented in Figure 5.
|Gene name||Amplicon (bp)||Sense||Sequence (5'to3')||Length (bp)||Tm|
|Heat shock protein 90 B||232||F||CGACTTCATCAAGCTGTTTGTC||22||59.9|
Table 2: Primer sequences, product sizes and annealing temperature (Tm) of selected genes for RT-qPCR.
The histological results and fish performance are reported as mean ± standard error of the mean (sem). After ensuring normality (Shapiro-Wilk, Skewness and Kurtosis tests) and homoscedasticity (Bartlett test), data were analyzed using one-way ANOVAs. For the length of trabecular bone observed with Von Kossa, data were Log10 transformed. The results were assumed to be significant at p-values below 0.05. All the tests were conducted with R software version 3.0.2  using Moments package .
This experiment employed all-female triploid individuals fed with a single P-deficient diet in order to overcome possible factors increasing phenotypic variation such as differences in somatic growth, sex and sexual maturation. After 27 weeks of feeding low-P diets, mortality was similar and negligible (3.2%) in all tanks and no significant difference in growth performance was observed in control and experimental fish. Moreover, no difference in length, weight and condition factor (Table 3) were observed between experimental individuals (p=0.06 and 0.09, respectively) () and external visualisation could not detect any defects of conformation (gross deformity).
|Phenotypes||Sirius red (basalias)||Von Kossa (trabecular bone)||mCT|
|Non mineralized area (%)||Intermediary mineralized area (%)||Mineralized area (%)||Total area (mm2)||Mineralized area (mm2)||Osteoid (mm2)||Tt.B.Ar (%)|
|Normal||59.4 ± 11.5||32.0 ± 9.6||8.5 ± 1.9||5931.2 ± 249.8||5711.0 ± 233.1||220.1 ± 16.7||39.3 ± 3.2|
|Deformed||68.6 ± 7.8||17.7 ± 2.0||13.7 ± 6.0||5442.0 ± 110.1||5166.1 ± 99.8||275.9 ± 10.7||43.3 ± 3.6|
|Deformed||68.6 ± 7.8||17.7 ± 2.0||13.7 ± 6.0||5442.0 ± 110.1||5166.1 ± 99.8||275.9 ± 10.7||43.3 ± 3.6|
|Phenotypes||Fork length (cm)||Weight||Condition||Stiff.||P Vt.||Ash Vt.||P Sc.||Ash Sc.||P plasma|
|Normal||35.1± 0.6||657.3± 33.2||1.5 ± 0.0||140.2 ±6.3*||9.3 ± 0.3||59.4 ± 0.5*||5.5 ± 0.1||29.8 ±0.7*||11.5 ± 0.4*|
|Deformed||32.6 ± 0.8||526.9 ± 49.6||1,6 ± 0.1||80.3 ± 9.0*||8.8± 0.3||55.0 ± 0.4*||4.4 ± 0.5||26.3 ±1.1*||7.6 ± 0.7*|
Table 3: Growth, condition factors, stiffness, P mineralization status (dry tissues) calculated as ash or P content and results of histological observations by Sirius red and Von Kossa. Values are expressed as means (n=3) ± sem. P: Phosphorus; Sc: Scales; Vt: Vertebrae; Tt.B.Ar : Total Bone Area. Asterisks (*) represent significant difference (P<0.05).
Here, we focused on pronounced biconcave phenotype because of its high frequency in salmonids prone to mineral deprivation and because it is considered as an initial development stage of more severe compression deformities [20,28]. Pronounced biconcave vertebrae after 27 weeks were observed in 24.9% (N=113) of the fish. Only 5.5% (N=25) of individuals displayed normal vertebrae (healthy control fish) all along experiments (i.e., from week 9 to week 27). Other fish showing normal phenotype at week 27 (31.1%; N=141) were not considered because they were known to recover from an abnormal phenotype at week 18. Moreover, other types of deformities such as homogeneous compressions (29.4%; N=133), small and widely spaced vertebrae (5.5%; N=25) and other types of deformities (3.5%; N=16) were also observed but were not considered for this study.
The P status of deformed fish displaying pronounced biconcave vertebra was negatively affected as shown by the lower ash content in their scales and vertebra, and the lower plasma Pi concentration (Table 3). Our results are consistent with previous observations that showed a significant correlation between low mineral content in scales and higher occurrence of vertebral deformities in rainbow trout fed low-P diets . Whether the altered mineral status resulted from differences in individual P requirements or from difference in feed intake due to limited food supply inherent to feeding approaches and/ or competition between congeners deserves to be further studied. The lower mineralization of biconcave vertebrae resulted in lower biomechanical competence when submitted to a compression force (stiffness, Table 3), a finding which is coherent with similar results reported in Atlantic salmon [47,48]. Lower P content and mechanical strength are recognised as initial indicators of deformity occurrence, particularly for the various compression types . Phenotypically, pronounced biconcave vertebrae differ from normal vertebrae by their apparent "transparent" vertebral body endplates upon X-ray evaluation (Figure 1); the shape of the vertebral body appears as a vertical apple core and the X-shaped centrum of the vertebrae display more acute angles [20,28]. Therefore, our results may suggest that the lower biomechanical competence of the vertebrae promotes the alteration of the vertebral body endplates (apposition site of newly formed bone) as well as in binding the X-structure of the vertebral body (more acute angles). If this hypothesis is true, the mechanisms of vertebral compression would not only affect intervertebral tissues at the vertebral body endplates as suggested by some authors [23,28], but also the entire structure of the vertebral body.
A recent study from our group reported the first reference transcriptome for rainbow trout vertebral bone and its characterization revealed that key regulators of bone metabolism have been well conserved in vertebrates and are present in trout . Here, we compared the overall gene expression by mapping paired-end reads on the transcriptomes in pronounced biconcave (deformed) vertebra with normal vertebral phenotype after 27 weeks of feeding a P deficient diet (Table 4). Using edgeR and DESeq package (FDR<5%; P<0.05), a total of 1289 genes were found differentially expressed. Of these, 700 were down- and 589 up-regulated in biconcave vertebrae. A blast-x search (cut-off 1e-8) was conducted to identify these sequences in Uniprot-Swissprot database reference and only genes with unique match were conserved. A total of 793 sequences (61.5% of the total genes) corresponding to 485 down- (61.9%) and 308 up-regulated (38.8%) genes were associated with known proteins. Furthermore, from the 1289 genes DE, 23.1%(N=298) were associated with putative lncRNAs. The unidentified sequences were not retained in our analyses although we believe that the study of these non-referenced sequences may improve our understanding of the process involved in deformity appearance.
|Ave. number reads||44.6 M||46.1 M|
|Ave. number trimmed reads||35.3 M||37.8 M|
|Ave. mapping (%)*||71.7||71.6|
|DE analysis||Def. Vs. Norm.|
|Up-regulated in Def.||700|
|Down-regulated in Def.||589|
FDR: 5%, triplicates. Unique gene match in Uniprot-Swissprot database (cut-off evalue<1e-8). Analyses on genes levels expression with both edgeR and DESeq packages.
*Trimmed reads with at least one reported alignment using RSEM output.
Table 4: Sequencing results of biconcave (Deformed, Def.) and normal (Norm.) vertebra. DE: Differentially Expressed.
Gene function analysis using the Ingenuity platform revealed that the top three networks differentially represented in this study were “cellular morphology, embryo development, hair and skin development and function”, “organismal development, skeletal and muscular system development and function in cancer” and “hereditary disorder, neurological diseases, organismal injuries and abnormalities”. They comprised 30, 28 and 28 genes, respectively. In addition, pathway analysis showed that a wide range of canonical pathways were differentially represented (p<0.05). They included, among others, protein kinase A signalling (p=1.2E-4), actin cytoskeleton signalling (6.5E-4), epithelial adherens junction signalling (8.5E-4), axon guidance signalling (1.7E- 3), paxillin signalling (5.4E-3), notch signalling (6.9E-3), inhibition of metalloproteinases (7.8E-3) and FGF signalling (1.9E-2). Enrichment analysis using GOminer highlighted 123 clusters with enrichment >2 and FDR<5%. Within these clusters are included relevant biological processes, cellular components and molecular functions involved in bone regulation: “GO:0008147 structural constituent of bone”: 23.1; “GO:0003308 negative regulation of Wnt receptor signalling pathway involved in heart development”: 15.4; “GO:0051019 mitogen-activated protein kinase binding”: 13.86; “GO:0030520 estrogen receptor signalling pathway”: 6.08; and “GO:0051216 cartilage development”: 2.1. Integration of functional analysis coupled to manual screening pointed to putative relevant genes involved in bone metabolism regulation (Table 1).
To understand the mechanisms involved in the depletion of the mineral status linked to the apparition of biconcave vertebrae in P-deficient fish, we combined RNA-seq experiments with recently developed morphometrical techniques. Our first hypothesis was that the lower mineral content and stiffness of deformed vertebrae resulted from higher osteoclast activity (bone resorbing cells). In human osteoporosis, the loss of bone and mineral is commonly linked to higher bone resorption by osteoclasts . Similarly, RANKL-induced ectopic osteoclast activity led to osteoporotic phenotype in medaka (Oryzias latipes), a teleost fish . Mature osteoblasts produce RANKL, which stimulate osteoclast proliferation, differentiation and activation, but they also synthesize osteoprotegerin (OPG), which is the major inhibitor of RANKL/RANK ligation . Here, no differences were observed between both OPG and RANKL gene expression. However, RNA-seq results showed that Tumor Necrosis Factor Receptor Superfamily, Member 5 (TNR5/CD40) was up-regulated in the deformed group. The ligation to TNR5 receptor induces RANKL pathway activation and bone resorption . To verify whether the lower mineralization in deformed trout vertebrae resulted from higher bone resorption, we verified Tartrate Resistant Acidic Phosphatase (TRAP) gene expression using RT-PCR. Expression of this marker of osteoclast activity was not significantly affected in deformed fish compared to normal individuals (Figure 5). However, in the former, unconventional Myosin 10 was down-regulated. Myosin 10 is a crucial protein for podosome belt formation allowing osteoclast attachment . These results are consistent with bone profiling assessment using μCT, which did not indicate major differences in bone resorption when comparing deformed vertebrae to those from normal individuals (Figure 3). Therefore, the lower mineral status in fish with deformed vertebrae did not likely result from higher osteoclast activity. On the contrary, it appeared that low levels of unconventional myosin 10 might impair osteoclast activity after 27 weeks of feeding low-P diet. Thus, our study suggests that impaired vertebral mineralization resulted from either alteration of matrix production and/or mineralization processes.
The mechanisms by which osteoblasts and osteoclasts control bone formation are not fully understood. However, various metabolic pathways affecting bone and cartilage formation have been identified in mammals and recent whole-transcriptome annotation provided evidence that these pathways were also present in teleosts, which suggests high conservation during osteichthyan evolution [26,27]. Genome-wide experiments revealed that most of the pathways (e.g. TGFβ/BMP, WNT and ERK/MAPK pathways) are involved, and may putatively cross-talk, during transcriptional regulation of osteoblast differentiation, proliferation and maturation [54-58]. Our results suggest that the pathways involved in the regulation of osteoblast differentiation and proliferation are down-regulated in deformed fish. Indeed, the following genes were down-regulated in deformed vertebrae: co-receptor activators or receptors (TGFβ/BMP: BMPR2; WNT: LRP6 and MACF1) and downstream actors (ERK/MAPK: NIK, MEKK1, SOS1, MYCB2 and MYCPP; TGFβ/BMP: USP9X). Conversely, inhibitors of ligand binding receptor were up-regulated (TGFβ/BMP: NOG, TEFF1 and FBN2; WNT: SFRP1). For instance, NOG expression was substantially up-regulated (fold change>9) in deformed compared to normal fish. In transgenic mice, over-expression of NOG triggered osteopenia, i.e. the abnormally low mineral status of bone tissue [59,60]. Finally, most of these signal pathway are either directly or indirectly controlled by a master regulator of osteoblast differentiation: RUNX2 . The expression level of RUNX2 was similar in deformed and normal vertebrae, but conflicting results were obtained regarding RUNX2 regulation. Indeed, in deformed fish we observed a lower expression level of HIVEP3, a gene negatively regulating osteoblast differentiation by stimulating degradation of RUNX2 and acting as a damper for signalling pathways such as WNT, TGFβ and BMP2 [62-66]. TWIST2, a RUNX2 inhibitor, was up-regulated in deformed fish while activators, such as STAB2, were down-regulated . In vitro experiments have shown that osteoblastogenesis is a physiological process that spends several weeks, during which specific actors are expressed at different times/stages . Yet, in deformed fish RUNX2 expression was not significantly affected, a finding that may be consistent with the higher expression of osteocalcin (also known as bone Gla protein, BGLAP), which normally occurs several weeks after a peak of RUNX2 .
BGLAP expression is a marker of mature osteoblasts. In deformed fish, the higher expression of BGLAP is concomitant with the higher expression of JUND, an AP-1 transcription factor complex, responsible for binding to BGLAP promoter and enhancing BGLAP expression in mammalian osteoblasts . In vitro experiments showed that JUND expression was relatively low during osteoblast proliferation and then increased during matrix maturation to reach a peak during matrix mineralization .
Taken together our analyses of biconcave vertebrae of deformed fish suggest that the processes involved in osteoblast differentiation and proliferation are reduced while the high level of BGLAP expression indicates the presence of numerous mature osteoblasts. We observed consistently in these deformed fish a low expression of COL12a1, a collagen type highly expressed by differentiating osteoblasts but not in mature osteoblasts, in which its expression is reduced . Furthermore, previous studies report that COL12a1 knockout mice displayed altered vertebral structure with decreased bone strength , which is coherent with the lower bone stiffness observed in deformed fish. Interestingly, the same authors reported no difference in TRAP expression or in other markers of bone resorption in COL12a1 -/- mice vs. wild-type control, indicating that decreased bone mass resulted from abnormal bone formation rather than from changes in bone resorption . The expression of COL11a1, a gene coding for a structural component of cartilage matrix, trabecular and compact bone as well as notochord tissue was also down-regulated in deformed fish compared to normal individuals. In mammals, COL11a1 ensures cohesiveness of collagen matrix and its altered expression results in thicker fibers and increased mineralization . In fish, COL11a1 is mostly present in notochord tissue, and is also expressed by osteoblasts in trabecular bone . Interestingly, in Atlantic salmon a long photoperiod resulted in a lower amount of COL11a1 in compact bone of compressed vertebrae . COL11a1 expression was also reduced in Atlantic salmon subjected to a combination of low P-dietary treatment and long photoperiod .
Nevertheless, at the histological level our study demonstrated that both trabecular bone area and the surface of mineralized bone areas were not significantly different (p=0.39 and 0.32, respectively) in pronounced biconcave compared to normal vertebrae (Table 3). We conclude that if the capability of osteoblasts to form new osteoid was not differentially affected in P-deficient fish with either biconcave or normal vertebrae, then the observed differences in mineral content might thus rather be related to the alteration of mineralization degree of the bone matrix.
Bone mineralization corresponds to the precipitation of calcium phosphate as hydroxyapatite (HA) crystals within a collagen-rich matrix, although how mineral formation is controlled is not fully understood. The combination of a favourable environment along with removal of mineralization inhibitors are known to be key factors regulating crystal formation . Similarly, the interaction of HA crystals with the organic matrix and their growth and expansion have not been fully elucidated.
However, we know that bone matrix, including non-collagenous proteins such as BGLAP, matrix Gla protein (MGP), osteopontin (also known as secreted phosphoprotein 1, SPP1) and SPARC, also control mineral deposition in bone tissue [77,78]. As said above BGLAP expression was higher in biconcave than in normal vertebra with regard to RNA-seq experiment as well as by confirmation with RT-qPCR (Figure 5). In mammals BGLAP is known to be expressed exclusively by mature osteoblasts and osteocytes [68,70] and constitutes the major non-collagenous protein in bone tissue. The role of circulating BGLAP is particularly relevant for energy metabolism and glycemia regulation . Interestingly, we observed a substantial up-regulation (FC>34) of GCG2, involved in lipid and glycogen hydrolysis and higher calcium intake . However, despite extensive studies, the direct role of BGLAP in bone tissue is still subject to debate. For instance, BGLAP deficient mice displayed no detectable impact on bone mineralization using Von Kossa staining , but infrared microspectrometry showed its positive role in the maturation of hydroxyapatite crystals . In juvenile Atlantic salmon with high growth rate BGLAP expression was negatively affected and the activity of mature osteoblasts and mineralization were reduced . Low BGLAP expression could result in a default in mature osteoblast function, hence triggering default in mineralization status. However, by controlling the direction of crystal growth or even inhibiting it, high BGLAP expression may also have negative impacts on mineralization [77,83,84]. The up-regulation of BGLAP is also consistent with a higher expression of C/EBP-β, a putative positive regulator of BGLAP, a finding which was also reported in rat . Interestingly, vitamin D is known to be a regulator of C/ EBP-β expression and has also been involved in regulation of BGLAP expression [70,86]. Similarly, vitamin D and expression of its receptor, VDR are affected by P-level in the diet [48,75,87,88]. In rainbow trout vitamin D induces NaPi-II expression, which is responsible for P absorption and retention . However, conflicting responses of P dynamics to vitamin D treatments have been observed in fish, indicating differences between teleosts and mammals [87,89]. MGP, encoding a protein involved in the regulation of the mineralization process, was up-regulated in deformed fish. MGP is a potent inhibitor of bone and cartilage mineralization [90,91]. In contrast to BGLAP, MGP is expressed in various tissues and is not limited to bone tissue in Atlantic salmon . The expression of MGP in deformed fish is negatively correlated to bone ash content and mechanical quality (stiffness) in deformed fish.
The ossification by matrix proteins was related to expression of DCHS1 and FAT4, both being lower in deformed than in normal fish. Interestingly, studies in transgenic mice indicated that the two proteins act as ligand-receptor complexes; Mice in which DCHS1 and FAT4 were invalidated show defects in vertebral column development and abnormal patterns of skeletal ossification . The mechanisms by which these proteins play a role in skeletal tissue are not fully understood. However, according to the authors, DCHS1 and FAT4 could play a role in interacting with planar cell polarity effector, a pathway that was significantly impacted in the current study. Taken together, these results highlight various actors, including matrix proteins that could be involved in the default of degree of mineralization observed in deformed fish and thus putatively triggering the development of deformities.
Finally, under normal conditions, the teleost skeleton displays a wide spectrum of intermediate skeletal tissues (from cartilage-like to bone-like tissues), which are complex to study [94,95]. In teleost vertebrae, more or less mineralized healthy cartilage as well as ectopic cartilage can be observed, the former being related to a metaplastic shift of bone and co-transcription of chondrocytic and osteogenic markers [23,24,96,97]. During development, the cartilage-type tissue localized in the basalia is replaced by bone-type tissue by a process called endochondral ossification. This process requires chondrocyte hypertrophy and apoptosis, blood vessel invasion and ECM degradation. It was thus important to assess whether molecular activity or overall mineral status might be biased by the high cellular activity within the basalias. Both osteoblasts and chondrocytes differentiate from mesenchymal stem cells and thus share numerous key events involved in their differentiation and regulation . Here, we observed that pathways involved in chondrocyte differentiation were affected in deformed fish (down-regulation of SHH and FGFR4) but not the expression level of a key mediator of cartilage formation SOX9  and of MEF2C, a marker of chondrogenesis [23,99], which remained constant. Expression of MEF2C was further confirmed by RT-qPCR (Figure 5). Furthermore, NKX3.2, a transducer of SOX9 signals was upregulated in deformed fish [67,100]. Our results also suggest that the invasion of blood vessels in the cartilage matrix, its degradation and replacement by osteoblasts might be negatively impacted in deformed fish, even if a distinctive phenotype is lacking. Indeed, deformed fish displayed impairment of the VEGF signalling pathway as indicated by lower expression of HSPG2, GPC1 and VEGFR1. The VEGF pathway is involved in angiogenesis (invasion of blood vessels) and ossification [101-103]. HSPG2 -/- mice displayed strong impairment of cartilage formation, ultimately leading to diminished osteoblast replacement and endochondral bone formation . Interestingly, in deformed fish we did not observe differences in the mineralization status of basalia. Consistent with observations of HSPG2 -/- mice, MMP10 expression was increased in deformed fish. MMP10 is expressed in chondrocytes degrading the cartilage matrix as well as in some, but not all osteoclasts and is also responsible for activation of several other MMPs including MMP9, that is involved in the vascularization process of the growth plate [105-107]. Here however, we observed an over-representation of canonical pathways for the inhibition of MMPs, hence the complexity to interpret the role of MMPs in deformed vertebrae. Thus, our results suggest that basalia deformed vertebrae might experience a subsequent endochondral ossification but with no impact of the calcification of the basalia by the time of sampling. To confirm the gene expression results, we used Sirius red staining coupled with image processing to extract the basalias and quantitatively assess the presence of mature collagen fibers (putatively mineralized) and that of less mature collagen fibers. In both groups, we observed that the vertebral arches were not fully ossified, but no significant differences were observed in their mineralization status (Table 3). Furthermore, no obvious presence of ectopic cartilage or remodelling of the notochord was observed in either deformed or normal fish, as previously reported in severe compression type deformities in other salmonids [24,96,108].
Interpreting data from RNA-seq experiments, particularly in a complex tissue like vertebra, remains challenging and is limited to the gene expression levels. This approach however, provides an unprecedented overview regarding molecular changes that are involved in the occurrence of vertebral deformities in P-deficient trout. The aim of our study was not to describe the impact of P-deficiency in trout but rather to characterize discrete differences between fish fed low-P diets presenting either deformed or normal vertebral phenotypes. We observed that the initial appearance of malformations in rainbow trout is associated with changes in numerous signalling pathways similar to those described in mammals, confirming their evolutionary conservation. Our results also suggested that after 27 weeks feeding low P diets, deformed vertebrae showed a decrease in osteoblast differentiation and proliferation while they display a higher expression of BGLAP, a marker of mature osteoblast activity. Together with histological observations, this suggests that during the weeks prior to sampling, deformed fish were still differentiating osteoblasts and producing matrix. However, higher expression of MGP and BGLAP as well as lower circulating P significantly impaired mineralization. Furthermore, impaired cohesiveness together with low mineralization reduced vertebrae stiffness. Consequently, reduced vertebral integrity might increase vertebrae susceptibility to deform under normal physical forces. The conserved osteoclast activity together with the production of BGLAP and GCG2 suggests that deformed vertebrae are undergoing significant bone remodelling with high bone resorption process not yet observable in our study. This process may lead to new vertebral phenotypes. As suggested in previous studies [20,28], pronounced biconcave vertebrae may be a transitory step toward a more dramatic abnormal vertebral phenotype (homogeneous compressed). Similar analysis of vertebral gene expression over time will highlight molecules involved in the occurrence of vertebral deformities. Our study provides consistent evidence of the importance of the matrix proteins and other markers involved in reduced mineralization and mechanical strength of vertebra, recognised as initial indicators of deformity development when feeding low-P diets to fish . Previous family-based studies in fish showed that genetics has an impact on the prevalence of deformities [109-111]. With the aim of reducing P excretion to promote sustainable freshwater aquaculture production, the indicators revealed in our study may ultimately serve as markers for genetic selection of fish strains displaying normal growth performance with a low prevalence of deformities, when feeding low-P diets.
We greatly thank the technical assistance of the Laboratoire des sciences aquatiques de l'Université Laval and Dr. Brian Boyle for his assistance with the Illumina library preparation at the Institut de Biologie Intégrative des Systèmes (IBIS), Québec. We thank as well Éric Normandeau for his help in bioinformatics analyses. Computations were carried out on the supercomputer Colosse, Université Laval, managed by Calcul Québec and Compute Canada. The assistance of Benjamin Bourdon regarding biomechanical analyses of vertebrae is acknowledged. This project was supported by the Ministère du Développement économique, Innovation et Exportation du Québec, the Ressources Aquatiques Québec (RAQ), the Société de recherche et de développement en aquaculture continentale (SORDAC) Inc. (PSR-SIIRI-443) and the Aquaculture Collaborative Research and Development Program, Fisheries and Oceans Canada (Q-08-01-001).