Differential gene expression and SNP association between fast- and slow-growing turbot (Scophthalmus maximus)

Growth is among the most important traits for animal breeding. Understanding the mechanisms underlying growth differences between individuals can contribute to improving growth rates through more efficient breeding schemes. Here, we report a transcriptomic study in muscle and brain of fast- and slow-growing turbot (Scophthalmus maximus), a relevant flatfish in European and Asian aquaculture. Gene expression and allelic association between the two groups were explored. Up-regulation of the anaerobic glycolytic pathway in the muscle of fast-growing fish was observed, indicating a higher metabolic rate of white muscle. Brain expression differences were smaller and not associated with major growth-related genes, but with regulation of feeding-related sensory pathways. Further, SNP variants showing frequency differences between fast- and slow-growing fish pointed to genomic regions likely involved in growth regulation, and three of them were individually validated through SNP typing. Although different mechanisms appear to explain growth differences among families, general mechanisms seem also to be involved, and thus, results provide a set of useful candidate genes and markers to be evaluated for more efficient growth breeding programs and to perform comparative genomic studies of growth in fish and vertebrates.

related to growth, identification of allelic variants of moderate to large effect could be useful in marker-assisted selection (MAS), provided a cost-benefit balance. Finally, genome-editing techniques may allow recovering or even introducing variants not present in the population of interest, increasing genetic gain 9 . As such, different studies have been focused on understanding growth differences in different fish species in the last decade [10][11][12][13][14][15] .
Turbot is the flatfish with highest aquaculture production, mainly held by Spain and China. Females reach sexual maturity at ~3 years (~50 cM length) and can spawn between 5 and 10 million eggs per season. Reproduction can be achieved along the whole year under a variable controlled photoperiod 16 . Thus, production can benefit from the high fecundity of a single selected couple over a relatively long period. Several studies have been carried out in this species to identify quantitative trait loci (QTL) associated with growth [17][18][19] , but the information provided, though useful, is still of limited application. Knowing the genes or even the allelic variants behind growth differences is important, not only for a more efficient selection, but also to guide future studies in turbot and other fish species. However, functional studies on differential growth have not been performed in turbot yet.
To gain knowledge on the key genes involved in growth, we carried out a transcriptomic study in fast-and slow-growing turbot to identify differential expression profiles and genetic markers between extreme growth phenotypes. Samples of muscle and brain tissues were taken from family-matched fast-and slow-growth fish selected from eight unrelated full-sib families, and their transcriptomes compared by RNA-Seq. These organs play a key role in the regulation of somatic growth, and muscle also represents the main edible part of fish 5,20 . In order to evaluate a high number of individuals and to avoid family bias, a pooling strategy was followed. Differential gene expression and marker association between high-and low-growth turbot were obtained. The detected candidate genes and markers were mapped using the recently assembled genome of the species and their positions were compared with previously reported QTLs 21 . The most interesting genes and markers were validated using individual Real-Time PCR (qPCR) and SNP genotyping to disentangle general and family-specific factors involved in turbot growth.

RNA sequencing output.
A total of ~183 million paired-end (PE) reads were generated, for an average of ~15 million reads for each of the twelve samples analyzed. After filtering, ~171 million reads were retained (93.44%; >13 million reads per sample). On average, 88.4% of the PE filtered reads aligned to the turbot genome at unique (83.9%) or multiple (4.5%) genomic positions, whereas the remaining 11.6% did not match. The resulting muscle and brain transcriptomes revealed a total of 30,674 and 49,034 transcripts, respectively, belonging to 17,173 and 23,104 genes (Table 1), which included some non-annotated genes in the turbot genome 21 . Annotation against NCBI's non redundant protein database was successful for 87.2% and 86.7% of the muscle and brain transcriptomes, respectively. Differential gene expression. A total of 17 and 4 differentially expressed (DE) genes were detected in muscle and brain, respectively, when comparing fast-(FG) and slow-(SG) growing pooled samples using the log2 FG/SG fold change (FC) ( Table 2). FC values were moderate, ranging from 1.18 to −0.86 FC values, positive FC corresponding to higher expression in FG compared to SG (up-regulated onwards), and negative FC to lower expression in FG compared to SG (down-regulated onwards). Since the number of DE genes was low, an extended list of those genes with FDR < 0.3 was also explored (Supplementary Dataset 1). Overall expression patterns were consistent across the three biological pools tested (Fig. 1), especially considering that FG and SG pools belong to eight different families and that the mechanisms driving growth differences are not necessarily the same for all the families. No significant sex-ratio differences between pools and families were found, supporting that sex was not a confounding factor in this study. KEGG enrichment analysis revealed the prevalence of certain pathways among the extended list of DE genes. The pathways consistently enriched were "Glycolysis/Gluconeogenesis" and "Starch and sucrose metabolism" (log 2 fold enrichment of 5.3 and 4.8, respectively) (Supplementary Figure 1), involving a set of 16 up-regulated genes coding for enzymes of the glycolytic pathway.
Among other interesting up-regulated genes in muscle, we detected six myosin heavy chain genes (FDR p-values ranging from 0.016 to 0.295 and FC from 1.00 to 0.48), strongly related to muscle growth; and the elongation factor 1-alpha (p-val = 0.099, FC = 0,70), which suggests enhanced protein synthesis for muscle growth. Interestingly, growth hormone receptor 2 (GHR2) was down-regulated (p-val = 0.063, FC = −1.01), showing an opposite pattern to that expected according to the function of the growth hormone receptor in mammals. Furthermore, the cytosolic creatine kinase (CKs) was up-regulated (p-val = 0.081, FC = 0.59), while its mitochondrial counterpart (CKm) was down-regulated (p-val = 0.170, FC = −0.67), in agreement with the higher production of cytosolic ATP due to the activation of the glycolysis pathway in FG fish. Accordingly, the mitochondrial oxidative phosphorylation showed a lower activity in FG turbot as suggested by down-regulation of ADP-ATP translocase 2 (p-val = 0.051, FC = −0.64), which codes for a membrane protein responsible for transporting ATP from the mitochondria to the cytosol. The most interesting genes detected in the family pools were evaluated in individual samples by qPCR: CKs, CKm, GHR2, lactate dehydrogenase (LDH), muscle glycogen phosphorylase (PYGM) and transforming growth factor β1 (TGFβ1). Although a notable concordance was observed between individual qPCR and pooled RNA-Seq data when comparing the expression pattern of FG and SG fish, some discrepancies were detected (Fig. 2), suggesting that some mechanisms driving growth might be different among families. Specifically CKm and TGFβ1 showed important variation among families. In fact, when individual qPCR data were pooled following the same distribution of individuals as for RNA-Seq, correlation was positive and significant between both datasets, except for CKm and TGFβ1 (Supplementary Table 1).
A lower number of DE genes were detected in the brain and not so straightforwardly related to growth as in the case of muscle. Among these, we observed down-regulation of phosphorilase kinase alpha 2 (p-val = 0.172, FC = −0.49) encoding for a subunit of the phosphorylase b enzyme, which is expressed in liver and brain. This gene modulates glycogen catabolism and plays an important role in supplying cell energy 22 . Also, we detected several up-regulated genes associated with sensory control, such as t-box brain protein 1 (p-val = 0.088, FC = 0.61), a transcription factor required for normal brain development and expressed in the olfactory bulb among other brain regions 23 ; synaptoporin (p-val = 0.039, FC = 0.66), a vesicular transport protein involved in synaptic vesicles and plasticity 24 ; and olfactomedin (p-val = 0.190, FC = 0.46), a secreted glycoprotein involved in regulating chemosensory cilia in olfactory neurons and optic nerve extension 25 . KEGG pathway and GO enrichment did not reveal any significant term.

Differential SNP allelic frequencies (DAfreq).
A total of 146,217 SNPs were called using both muscle and brain reads. The sequences containing these SNPs were annotated according to the turbot genome information 21 . A total of 45,348 SNPs (31%) were found in coding regions, and among these, 16,133 corresponded to non-synonymous mutations, representing an overall Ka/Ks ratio of 0.55. Allelic frequencies were estimated for FG and SG fish and the significance of allelic divergence tested (DAfreq onwards). Significant DAfreq between FG and SG were detected for 75 SNPs in muscle and for 358 in brain (Supplementary Dataset 2). These 435 DAfreq SNPs identified were located in 326 genes: i) 43% in untranslated regions (3′ and 5′ UTR); ii) 49.2% in coding regions of annotated genes, including mostly synonymous (74%), but also missense variants (24%), and 2% located in splicing regions; iii) 3.4% in introns, and iv) 4.4% in non-annotated genes. The genomic position of these SNPs in the turbot genome 21 and their co-localization with previously described growth-associated markers 17,18 was evaluated (Fig. 3). DAfreq SNPs were widespread throughout the genome, although some specific regions at particular linkage groups (LG) of the turbot genetic map showed a higher concentration (i.e. LG02, LG05 or LG13), and some of them co-localized with previously reported growth associated markers. We also  Table 2. Differentially expressed genes between fast-and slow-growing turbot. FDR p-val: false discovery rate corrected p-value; FC: log2 fold change; LG: linkage group; UN: gene located in a scaffold not assigned to a LG of the turbot genetic map.
checked for the correspondence between DAfreq SNPs and DE genes in the turbot map ( Fig. 3). A selection of the most suggestive SNPs is shown in Table 3. Several DAfreq SNPs co-localized in a narrow region of LG2, with two glycolysis up-regulated genes in muscle, glyceraldehyde-3-phosphate dehydrogenase and triosephosphate isomerase, including a synonymous SNP variant of the gene glycogen synthase kinase 3, which is involved in the negative regulation of skeletal muscle growth (Table 3; Supplementary Dataset 2). Highly significant SNPs showing a large DAfreq (>0.4) were also observed in two expressed genes in brain coding for sensory-related vesicle proteins, synaptotagmin 1 and synaptophysin, mapping at LG11 and LG16, respectively. Further, the muscle up-regulated deoxynucleoside triphosphate SAMHD1 gene showed a missense mutation located at LG10. Except for the SNP corresponding to fast myosin heavy chain protein, all the SNPs assayed were technically feasible and could be characterized using the SNaPshot methodology (8 out of 9). These SNPs were successfully genotyped in all FG and SG samples, and significant frequency differences for three novel candidate genes were individually validated (Table 3; Supplementary Table 2). Worth noting is the DAfreq observed for SAMHD1, which also presented differential expression between extreme growth phenotypes in this study.

Discussion
Transcriptomes of fast-(FG) and slow-(SG) growing turbot involving eight different families from a breeding population of Atlantic origin were compared in this study. Several genes showing differential gene expression (DE) and/or divergent allelic frequencies at SNP markers (DAfreq) were found when comparing both phenotypic groups. The use of pools of families represents a cost-effective strategy for a preliminary exploration of the most important factors involved in differential growth. This approach enables identifying candidate genes, pathways and genomic regions associated with differential growth rates, as recently reported in other aquaculture species [26][27][28][29] . This strategy was here applied for the first time in the turbot, the most important farmed flatfish worldwide 30 . We identified relevant candidate genes and genomic regions associated with differential growth in   this species, essential for understanding the mechanisms underlying phenotypic variation of this complex trait and useful for genetic assistance in breeding programs. The transcriptome profiles at 270 dpf here obtained also contribute key information to be further compared with other life-cycle stages with different weight of hyperplasia and hypertrophy into myogenesis 5 . Some candidate genes were individually validated using qPCR in FG and SG fish in turbot families. Aggregated individual qPCR results were mostly in agreement with RNA-Seq pools; however, some genes showed inter-family transcriptomic differences, suggesting that functional strategies associated with growth rate in turbot may be diverse among families, as previously reported for other complex traits in teleosts 31 . Our results are consistent with previous growth-association studies in turbot and other fish, where some genomic regions showed association across families, while others were family specific 17,18,32 . Furthermore, SNaPshot assays validated the association of some growth-related SNP variants identified through differential allele frequencies between FG and SG pools in the transcriptomic analysis. These gene markers could be useful for studies in other families and/or populations, or even for genome editing approaches for those potentially functional variants.
Up-regulation of all glycolytic enzymes along with lactate dehydrogenase was consistently observed in our study, which supports a higher oxidation rate of glucose by fermentation in muscle to increase growth rate. Unlike other vertebrates, fish muscle is separated in discrete layers of different fibre types. Vertebrates have two main types of striated muscle fibres, red and white, which in fish specialize in low speed swimming and bursts of maximum speed, respectively 33 . Red muscle fibres (less than 10% of the myotomal musculature) mainly depend mainly on aerobic, mitochondria driven, energy metabolism 34 . On the contrary, the metabolism of white muscle fibres (never less than 70% of fish skeletal muscle) relies mainly on anaerobic glycolysis, aided by the intracellular energy shuttle of cytosolic creatine kinase (CKs) 35 . The higher expression of glycolytic enzymes, lactate dehydrogenase (LDH) and CKs in fast-growing turbot strongly suggests a higher white muscle metabolic activity, which would presumably explain their higher growth. Glycolysis has also been found up-regulated in the muscle of domestic rainbow trout compared to their wild counterparts 36 , and also in faster-growing fish of the same species; this up-regulation has been associated with an increased muscle energy demand in fast-growing animals 29 . However, we cannot discriminate cause and consequence from our results. Either high growth rate demands an increasing anaerobic energy metabolism or, conversely, higher activity of glycolytic enzymes, LDH and CKs, results in higher growth rate. Further work is needed to clarify this issue.
Two suggestive candidates for follow-up experiments are the glycolytic enzymes triosephosphate isomerase 1 (TPI) and glyceraldehyde-3-phosphate dehydrogenase (G3PDH). Both genes were up-regulated in turbot muscle and they are closely linked in LG2 (181 kb; ~0.3 cM 21 ). Further, these two genes lie within the confidence interval of a growth-related QTL in turbot 17,19 , in the vicinity (<500 kb; ~1 cM 18 ) of several DAfreq SNPs. TPI markers have been found associated with length variation in other fish 37 , and its expression has also been correlated with body size, muscle fat content and growth potential in bulls 38 , as also observed for G3PDH in fast-growing rainbow trout 29 . Other clusters of a few DE genes in short genomic stretches (<1 Mb) were observed in our study, suggesting positional clustering of co-expressed genes of related pathways as reported for Arctic char 31 . The integration of structural and functional data provides clues for understanding the genetic architecture underlying growth in turbot and fish in general. The co-localization of several DE genes with DAfreq SNPs and growth-related QTL suggests that multiple genes could be underneath these growth-related genomic regions in turbot, as observed for complex traits in other fish 31 .
A higher activity of the glycolytic pathway was also supported by the up-and down-regulation of cytosolic (CKs) and mitochondrial (CKm) creatine kinase, respectively, in the muscle of turbot. Creatine is used as an energy reservoir in its phosphorylated form, and CK is responsible for creatine phosphorylation, capturing LG16 9,519,616 C/A 0.39 (7.9) 0.05 synaptophysin protein 1 5′UTR LG24 2,244,489 A/G 0.37 (7.8) 0.07 glucose-6-phosphatase 3 Synon Table 3. Selection of SNPs showing divergent allelic frequencies between high and low growth turbot. For each SNP the chromosome (Chr), the position in the chromosome (Pos), the major variant (V1), the minor variant (V2), the allelic frequency difference (DAfreq) and Fisher's test −log10 p-value (−log pval) obtained from the RNA-seq, the DAfreq and significance of genotypic/allelic divergence (G-tests; s: p < 0.10, *p < 0.05, **p < 0.01) obtained from SNaPshot, the gene the SNP belongs to (Gene) and the SNP annotation regarding the gene (Annotation; Synon: synonimous) are shown. available cellular energy; but it can also catalyse the reverse reaction, using phosphocreatine to regenerate ATP, acting as a metabolic regulator 39 . While CKs activates glycolysis by removing the ATP generated, CKm is coupled to the mitochondrial oxidative phosphorylation capturing the ATP generated 40,41 . In this context, the up-regulation of CKs observed in FG turbot is likely consistent with the higher anaerobic metabolism in the white muscle. Both CKs and CKm genes have been reported to be highly regulated in muscle between large and small rainbow trout, associated with reactive oxygen stress imbalance 29 . Growth hormone receptor 2 (GHR2), a fish-specific paralogous gene, was found down-regulated in the muscle of six out of eight turbot families, a similar expression pattern to that observed in fast-growing rainbow trout 29,41 . Growth hormone pathways are crucial for vertebrate growth, and specifically the variation at GHR2 has been also associated with growth in different tilapia species 42 . GHR2 is involved in fish myogenesis and compensatory growth 43 , although further functional analyses are required to understand its signal transduction mechanisms. Interestingly, GHR2 is located within the confidence intervals of a QTL explaining up to 18% of the phenotypic variance for body weight and length at LG14 in turbot 17 .
Unlike muscle, only a few DE genes were detected in brain when comparing FG and SG turbot. Detection of DE genes could be hindered by sampling the whole brain instead of dissecting specifically the hypothalamus or hypophysis, a very difficult task before the onset of sex growth dimorphism in turbot. Even so, to our knowledge, this is the first transcriptomic study showing brain differences between FG and SG in turbot. Interestingly, some up-regulated genes had been previously related to sensory control, including t-box brain protein 1, synaptoporin and olfactomedin, the latter two associated with domestication behaviour and feeding preference in other fish, and related to higher growth rate 44,45 . DAfreq SNPs were also found in two relevant genes coding for vesicle proteins, synaptophysin, associated with behavioural variation in domesticated zebrafish 44 , and synaptotagmin 1, related to memory formation and retention, and feeding response in fish 46 . Moreover, these SNPs mapped within the intervals of growth QTL explaining more than 10% of observed variance 17 , where major candidate genes are also located (lumican at LG16 and myogenin at LG11, respectively) 19 . Synaptotagmin 1 SNP was validated using SNaPshot, showing significant and suggestive differences between FG and SG fish at genotypic and allelic level, respectively. All these data suggest that sensory-related genes and pathways regulating learning and memory in brain might explain differences in feeding behaviour and, thus, be involved in turbot growth rate.
Numerous DAfreq SNPs were detected in muscle and brain transcriptomes distributed across all chromosomes of the turbot genome, even in linkage groups such as LG9 and LG22, where growth QTLs have not been previously reported 21 . Although many of them may reflect linkage disequilibrium with causal polymorphisms or be false positives, others may explain growth differences in turbot, especially if they co-localize with previously reported QTL, or better, if they are within candidate genes showing differential expression. Thus, this list represents a catalog of candidate markers affecting gene expression and/or protein function with small to moderate effects associated with growth. Some differential SNPs concentrated in genomic regions co-localizing with previous growth markers reported in turbot 21 , particularly focused on those found across different families, which may represent the QTL variation in the wild Atlantic population of origin 17,18 . For example, a cluster of differential variants were detected in LG12, close to 3/9CA15 explaining up to 9% of length variance; or in LG16, close to Sma-USC-E11, explaining up to 25% of Fulton's factor variance 17 , and apparently conserved in different fish species 21 . The largest accumulation of DAfreq SNPs is observed in LG2 where, as previously mentioned, there are also two glycolitic enzymes showing differential expression. One of the SNPs in these regions is located in glycogen synthase kinase 3 (GSK3). This SNP has been validated individually and shows significant genotypic frequency differences between FG and SG fish. GSK3 is involved in the negative regulation of muscle growth and differentiation 47 . A decrease in the activity of this protein is positively related to muscle growth. Further molecular and biochemical analysis are needed to determine if the identified SNP has an impact on GSK3 and consequently on the differential growth rates.
DAfreq SNPs within DE genes represent the most robust evidence for causality. In this study, DAfreq SNPs were often found in the vicinity of DE genes, but rarely mapped within DE genes. Two noteworthy exceptions were missense variants in the gene SAM and HD domain containing deoxynucleoside triphosphate (SAMHD1) located in linkage group LG10 and the myosin heavy fast skeletal gene located at turbot LG6. SAMHD1 was found to be up-regulated and its non-synonymous SNP individually validated, being significant at allelic (p = 0.034) and genotypic (p = 0.008) level. This gene encodes a major regulator of dNTP reservoir in the cell, playing an essential role in cell-cycle progression and cell proliferation, with impact on viral replication and cellular DNA polymerization 48,49 . SAMHD1 dNTP reduction has been associated with a reduction of the number of mitochondria per cell by impeding mtDNA replication 49 . This is consistent with down-regulation of CKm in fast-growing turbot. Furthermore, this gene mapped within the confidence interval of a suggestive QTL for growth in turbot 18 . On the other side, the myosin heavy fast skeletal gene was in the vicinity (<500 kb; ~1 cM) of a growth QTL marker 18 and was up-regulated together with other myofibrillar component genes, as also reported in fast-growing rainbow trout 29 . These observations support their interest as growth candidates for future studies.

Conclusion
The transcriptome comparison of fast-and slow-growing turbot from eight full-sib unrelated families enabled the identification of genes showing differential expression and SNP association with growth across the turbot genome. Sound candidate genes and genomic regions emerged as potential biomarkers of differential growth in turbot. Up-regulation of the glycolytic pathway in fast-growing turbot was consistent with a lower mitochondrial activity presumably related to higher proportion of white muscle fibres. A handful of sensory-related genes showing differential brain expression or divergent allelic frequencies suggest differential behaviour related to feeding between growth phenotypes. Nevertheless, family transcriptome variability of particular genes suggests that heterogeneous mechanisms might also be involved in differential growth rate, in accordance with the polygenic nature of this complex trait. Future studies should aim to expand these findings in other developmental stages to better comprehend growth differences in turbot. Nonetheless, all genomic resources gathered here have been useful for validating differential genetic variation between growth phenotypes and will aid in further comparative studies of growth in fish, providing a set of markers with an effect in turbot growth to be evaluated in future studies.

Methods
Families and sampling. Turbot from eight unrelated full-sib families coming from a pedigreed breeding strain of Atlantic origin were reared in tanks at 18 °C in the facilities of CETGA (Aquaculture cluster of Galicia; Ribeira, Spain). Fish were weighed at 270 days post fertilization (dpf), before significant differences in growth rate between males and females occured 50 . The 3-4 largest and smallest fishes from each family were sampled, rendering a total of 26 individuals for each phenotype: Fast-(FG; average weight = 402 g, SD = 51.38) and slow (SG; average weight = 180 g, SD = 35.65) growing turbot ( Fig. 4; Supplementary Dataset 3). We genotyped all animals using a sex-linked molecular tool 51,52 to identify putative sex-associated growth differences using contingency chi-square tests (p < 0.05). Muscle and brain tissues were collected and individually conserved in RNAlater RNA sequencing. RNA extraction was performed using the RNeasy mini kit (Qiagen) with DNase treatment following manufacturer's instructions. RNA quality and quantity were evaluated in a Bioanalyzer (Agilent Technologies, USA) and in a NanoDrop ® ND-1000 spectrophotometer (Nanodrop Technologies, USA), respectively. Three libraries were constructed for each growth group (FG and SG) and tissue (muscle and brain), totaling 12 libraries. Eight to nine individuals from 2-3 different families were pooled in each library: 3 FG and 3 SG libraries for each tissue (Supplementary Dataset 3), every family being present in both FG and SG groups. Pooled samples were barcoded and prepared for sequencing at the Wellcome Trust Centre for Human Genetics, Oxford, where 100 bp paired-end (PE) reads were obtained on an Illumina HiSeq. 2000. The quality of the sequencing output was assessed using FastQC v.0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Quality filtering and removal of residual adaptor sequences was conducted on read pairs using Trimmomatic v.0.32 53 . Specifically, Illumina adaptors were clipped from the reads, leading and trailing bases with a Phred score <15 were removed and the read trimmed if a sliding window average Phred score over four bases was less than 20. Only reads where both PE reads encompassed a length greater than 32 bp post-filtering were retained. Filtered reads were mapped to the recently assembled turbot genome 18 60 . Genes with Benjamini-Hochberg false discovery rate (FDR) corrected p-values < 0.05 were considered differentially expressed (DE) genes. Due to the low amount of significant DE genes, those showing FDR corrected p-values < 0.3 (standard p-values < 0.005 and 0.0005 for muscle and brain, respectively) were also explored to identify a wider list of candidate genes for investigating enriched functions and pathways. Log2 library-size normalized counts were used for generating heatmaps using "aheatmap" function in R/NMF 61 . Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was assessed using KOBAS 2.0 62 and Gene Ontology (GO) terms using Blast2GO 63 , the DE genes were compared to the whole muscle or brain transcriptomes and those terms or pathways showing FDR corrected P-values < 0.05 were considered enriched. Differential allelic frequencies. Pileup files describing the base-pair information at each genomic position were generated from the alignment files using the mpileup function of Samtools 64 discarding those aligned reads with a mapping quality < 20 and those bases with a Phred score < 20. SNP calling required a minimum read depth of 50 and a minimum of 5 reads for the least frequent allele. Differential allele frequency (DAfreq) between FG and SG samples were obtained using PoPoolation2 65 and the significance of DAfreq between the two groups was tested by Fisher's exact test. Those SNPs showing Bonferroni corrected p-values < 0.05 and further DAfreq > 0.25 were considered significant. All the SNPs were functionally annotated using SnpEff v.4.2 66 and the turbot genome protein annotation 21 . SNP physical location was compared with that of previously reported growth-associated QTLs in turbot 17,18 , based on other farmed families sharing an Atlantic origin with the CETGA families here studied. All these families come from broodstocks derived from the turbot Atlantic area that shows very low genetic differentiation 16 . Real-Time PCR. The same RNA samples used for RNA-Seq were analysed individually by qPCR. Primers for selected genes (LDH, PYGM, GHR2, CKs, CKm, and TGFβ1) were designed using the Primer 3 software 67 (Supplementary Table 3). To determine DE between conditions the ΔΔCT method was used with the slow-growth (SG) condition set as normalizer, and the ribosomal protein S4 (RPS4) and ubiquitin (UBQ) as housekeeping genes. These two genes were previously validated for gene expression determination by qPCR in turbot 68 . Reactions were performed using a qPCR Master Mix Plus for SYBR Green I No ROX (Eurogenetec) following the manufacturer instructions, and qPCR was carried out on a MX3005 P (Agilent Technologies). Analyses were performed using the MxPro software. A t-test was used to determine significant differences between conditions. Two technical replicates of each sample were included. SNP validation. PCR and SNaPshot primers for SNPs included in selected candidate genes were designed using the NCBI primer designing tool and blast checked against the nr database 69 . Primer length was modified to allow multiplexing for nine SNPs showing DAfreq (Supplementary Table 4). Total DNA was purified after protein precipitation (5 M NaCl) with freezing cold absolute ethanol (1 mL). Motility rate for each potential single-base extension product was calculated using SnaPshot Primer Focus kit (Thermo Fisher Scientific) and samples were divided into two multiplex reactions and one singleplex (Supplementary Table 4). A sample of DNA from each individual was amplified and the SNaPshot minisequencing reaction was carried out using ABI Prism SNaPshot Multiplex kit (Applied Biosystems) in an ABI Prism 3730 DNA sequencer according to manufacturer's instructions. The genotyping analysis was performed using GeneMapper ID software (Applied Biosystems). Exact G-tests were used to determine the significance of genotypic and genic differentiation between FG and SG population groups using Genepop 4.2 70 .