Estimation of genetic parameters for growth and carcass traits in turbot (Scophthalmus maximus)

Abstract Information on phenotypic and genetic (co)variance for production traits in turbot is required to improve breeding programs. So far, information on morphometric growth traits is sparse and completely lacking on quality carcass traits like fillet weight or fillet yield for turbot. As part of a long-term study we explored the phenotypic and genetic (co)variance of 16 biometrical and carcass traits of three different European turbot strains. Fish were reared under commercial grow-out conditions, including size grading. We used molecular relatedness (MR) methods based on genotyping with 96 microsatellite markers and animal models. We included an adapted condition factor for Pleuronectiformes (FCIPLN) and average daily weight gain (ADG) between the ages of 300 and 500 d post-hatch (dph) for their potential correlation with body weight at harvest. Heritability estimates for all traits were low to medium (0.04–0.29) when strains were jointly analyzed. Separate analysis of strains yielded higher heritability estimates (0.12–0.43). Genetic correlations between weight-related traits were highly positive (0.70–0.99), while runs with yield and ratio traits often resulted in unreliable estimates of genetic correlation due to high standard errors. Body weight (h2=0.19), fillet yield (h2=0.15), and dressing percentage (h2=0.17) are particularly promising selection traits for turbot breeding.


Introduction
Aquaculture production of turbot (Scophthalmus maximus) has been steadily increasing over the last years, especially in the PR China (64 000 t in 2013), which is the biggest producer worldwide (FAO, 2012). In Europe, turbot is mainly valued as a high-priced food fish on the upscale gourmet market, with the vast majority of production taking place in Spain and Portugal. Smaller aquaculture facilities for turbot are located in France, Norway, Denmark, and the Netherlands. Turbot is primarily sold whole, with a smaller portion of aquaculture production being sold either gutted or filleted (Turbot Species Leaflet -European Commission, 2016). Regardless of these different marketing strategies, the improvement and harmonization of growth traits remains one of the biggest and economically most relevant challenges of turbot breeding programs. Despite a recent rise in the development of genetic and genomic resources for turbot, reliable estimates of genetic parameters are still needed to improve existing breeding programs. Studies on heritability and genetic correlations in turbot have mostly been focused exclusively on body weight as the trait of the highest economic value (e.g., Gjerde et al., 1997;Guan et al., 2015) and have revealed moderate levels of genetic variation and therefore potential for selection. However, even when fish are marketed whole, traits like fillet weight, fillet yield, and dressing percentage are of great importance with regard to customer preferences and financial profits (Rutten et al., 2005) and have therefore been incorporated as additional selection criteria in many of the established breeding programs for other important aquaculture species (Neira et al., 2004). Additionally, carcass traits like the hepatosomatic index (HSI) and morphometric traits like length, width, or Fulton's condition index (FCI) can be used as readily obtainable indicators for health and fitness conditions. Possible interactions on a genetic level between traits of interest have to be quantified to avoid unfavorable correlated selection response. As far as we know, this is the first time that a quantitative genetic analysis of carcass traits in turbot has been done. As part of a long-term study we explored the genetic and phenotypic (co)variance of 16 biometrical and carcass traits of different turbot strains under commercial grow-out conditions using molecular relatedness (MR) and restricted maximum likelihood (REML) methods. In addition to standard traits like weight and length, we included FCI and average daily weight gain (ADG) between the ages of 300 and 500 d post-hatch (dph) to explore their potential correlation with body weight at harvest. Fulton proposed his condition index to measure fish health in 1904, assuming that the standard weight of the fish is proportional to the cube of its length. Since this assumption is based on the morphometrical traits of roundfish, we additionally included an adapted condition factor for Pleuronectiformes (FCI PLN ) in our analysis. The adapted condition factor relates the weight of the fish to the volume of a hypothetical cube, with its edges defined as √ length × width, and has been demonstrated by Arfsten et al. (2010) to be much better suited for use in flatfish than the original one.

Experimental design
Nearly 8000 turbot of four different European strains were raised from juveniles to approximately market age in two experiments to collect individual weight and morphometric growth data. During the first trial in 2010-2011 fish originating from Norway (strain A, N = 3423) and Iceland (strain B, N = 2563) were reared as part of a larger project on growth modeling (Oesau et al., 2013). Each fish was tagged intra-abdominally with passive integrated transponder (PIT) tags (Hallprint, PTY Ltd., Hindmarsh Valley, Australia). Wet body weight (BW, rounded to the nearest gram (g)) and a digital photograph were recorded for each fish individually using a semiautomatic camera setup (Fishback et al., 2002). The experiment lasted for 15 months and measurements were taken in 42 d intervals. Standard length (L), width (W), and area (A) were determined in centimeters from the photographs manually with the ImageJ software (Schneider et al., 2012). Length was defined as length from tail fin to snout, W was measured at the widest point of the fish excluding the fins, and A was defined as the entire body surface area of the fish excluding the fins. Measurement protocols were established in a previous experiment and are described in detail in Oesau et al. (2013). Starting in 2013, animals from France (strain C, N = 940) and Norway (strain D, N = 1026) were subsequently raised and measured in 30 d intervals for a duration of 18 months, utilizing the same procedures. Both experiments took place in two identical recirculating aquaculture systems (RASs) containing 10 round tanks with a total water volume of 40 m 3 each at the Gesellschaft für marine Aquakultur mbH (GMA) in Büsum, Germany. Water temperature was kept constant at 17 • C (±1). Fish were fed a commercial diet (Aller Aqua, Denmark). Further details on feeding and water parameters are described by Lugert (2015) and Oesau (2012). For management reasons fish were graded based on body weight between three and six times during the experiment. This is common practice under commercial grow-out conditions and serves to establish tanks with same-size cohorts to minimize social hierarchies and facilitate feeding management. Each animal was assigned a new number reflecting their contemporary group (CG) after each grading event. CG numbers were unique between the two experiments. CG thus reflects common environment factors such as year, size cohort, and stocking density and has been included as an adjustment in modeling. Gonadal sex and carcass traits were recorded postmortem on an individual basis. After a fasting period of 5 d fish were killed adhering to animal welfare standards using a mixture of ice water and CO 2 , then carved and filleted manually. Since this process was time-consuming and took place on several consecutive days, the age at slaughter differs within strains. Age at slaughter also differs between strains due to nonsynchronized European production cycles, which resulted in different delivery dates and ages of fry. All further analyses have therefore been adjusted to include age at slaughter as a fixed effect. Four filets were cut from each fish, two each on the dorsal and ventral sides. A total of 26 different persons were involved with filleting, and preliminary tests showed large inter-individual variation in fillet yields due to personal experience levels and susceptibility to fatigue. The filleting person (FP) has therefore been recorded for each fish and included as a fixed effect in the models. Individual weights have been collected for fillet without skin (FW), head (HE), liver (LI), gut and stomach (GU), and the rest of the carcass (CA) including bones, skin, and gonads. Total weight of intestines (IW), gutted weight (GW), hepatosomatic index (HSI), fillet yield (FY%), adapted condition index for Pleuronectiformes (FCI PLN ), and dressing percentage (DP%) have been calculated according to the following formulas.
average daily gain: ADG * = BW 2 − BW 1 dph 2 − dph 1 (8) BW 1 is body weight at previous measurement, BW 2 is current body weight, dph 1 is age in days post-hatching at previous measurement, and dph 2 is age at current measurement. Combining all records from both parts of the experiment resulted in a phenotypic dataset of 7952 turbot, with the larger part originating from the 2010-2011 time period (strain A: N = 3423, strain B: N = 2593). Since strains A and B were originally part of a larger project with a different research objective the availability of fish for genotyping differed between experiments. In this study genotyping efforts thus had to be focused on animals of strains C (N = 940) and D (N = 1026). Animals were removed from the dataset if the number of growth recordings per fish was fewer than three (N = 69), if gender could not be reliably determined (N = 146), or if the animal showed obvious signs of illness based on visual examination during slaughter (N = 79); 534 fish of strains A and B had been randomly selected in 2010, resulting in a total number of 2500 genotyped animals. A further 344 individuals had to be excluded later due to poor genotyping quality, leading to a total dataset of 2156 turbot with ages at slaughter ranging from 496 to 758 dph (Table 1). For strain-wise analysis, we excluded stain B due to an overall low contribution to the dataset and combined strains A and D since they originated from the same breeding facility in Norway.

Markers and genotyping
Lysates with a high concentration of genomic DNA were produced from fin samples according to a modified protocol by Miller et al. (1988). Fish were genotyped using a panel of 94 polymorphic microsatellite markers ( Table S1 in Supplement) that were tested in a small subset of fish and chosen based on their high polymorphism information content. Markers were selected to be evenly distributed across the genome (Bouza et al., 2012), covering 22 of 24 linkage groups with an average inter-marker distance of 10.5 cm. Primer pairs were designed with the Primer 3 software (Untergasser et al., 2012). Primer sequences and marker information are presented in Table S1. Multiplex PCR amplifications were performed using fluorescence labeled primer pairs, the Qiagen Multiplex PCR kit (Qiagen GmbH, Hilden, Germany), and an MJ Research PTC-200 thermal cycler (Global Medical Instrumentation Inc., Minnesota, USA) under the following cycling protocol: initial denaturation at 95 • C for 15 min, 40 cycles each of 94 • C for 30 s, 56-64 • C for 1 h 30 min (depending of the optimal annealing temperature for each multiplex set), 72 • C for 1 h 30 min, and final extension at 72 • C for 20 min. Capillary electrophoresis and fragment length analysis were carried out with an ABI 3130xl Genetic Analyzer and the Genemapper 4.0 software (Thermo Fisher Scientific Inc., USA).

Molecular relatedness (MR)
All data handling and statistical computing were done in R (R Development Core Team, 2011) unless otherwise specified.
As neither pedigree information nor DNA samples of possible parents were available for strains C and D, we estimated molecular relatedness for the total dataset based on co-ancestry using the similarity index described by Eding et al. (2001): where S xy,l is the similarity index for individuals x and y at locus l, and I ij is an indicator variable that is 1 when allele i at locus l in the first individual equals allele j at the same locus in the second individual; otherwise, it is 0. The values were averaged across loci and according to Hayes and Goddard (2008), and each element of the resulting matrix was transformed as where min is the minimum similarity in the matrix. This is an estimate of co-ancestry, and molecular relatedness r xy was calculated as 2S xy (Lynch, 1988;Blonk et al., 2010). The resulting co-ancestry matrix (r xy ) was then used to calculate genetic parameters instead of conventional pedigree (PED).

Genetic analysis
For the estimation of heritability we used a univariate additive animal model: where y ij klm is the mth observation of BW, L, A, or W, µ is the general mean, ST i is the fixed effect of strain i (i = 1 to 4), CG j is the fixed effect of contemporary housing group j (j = 1 to 64), S k is the fixed effect of gonadal sex k (k = 1 (male) or 2 (female)), b n (dph) is the fixed regression coefficients on days post- , and x ij klm5 (dph) = 1/8(63 dph 5 − 70 dph 3 + 15 dph), animal l is the random additive genetic effect of animal l (l = 1 to 2156), and e ij klm is the random error term.
The model includes the fixed effects of strain, contemporary housing group, and gonadal sex. The filleting person was included in the model only for traits that were affected. The age curve was modeled by applying Legendre polynomials of degree 5. Polynomials of this order were chosen since we experienced convergence issues in some traits with lowerorder polynomials. Polynomials of degree 5 yielded consistent results across all traits. We utilized the restricted maximum likelihood method as implemented in the WOMBAT software package (Meyer, 2007) to calculate variance components. The random effect of "animal" is distributed proportionally to the "arbitrary" relationship matrix, which has been supplied to WOMBAT as the inverse. Molecular relationship information was supplied to WOMBAT as a positive-definite, general inverse matrix file. Fixed effects were tested for significance for each trait by stepwise forwards selection comparing Akaikes information criterion and Bayes information criterion (AIC and BIC) values and by chi-square statistics (significance level: p = 0.05) in the model adding one trait at a time. Estimates of correlations were obtained using the same model in a bivariate setting, since we experienced convergence issues with multivariate settings. Heritability estimates (ratio of additive genetic variance to phenotypic variance h 2 ) and genetic correlations including standard errors could be directly retrieved from WOMBAT output files.

Descriptive statistics
Mean weights and ages for all strains at the start of the experiment are given in Table 1: fish of strain A (N = 261) were the oldest and thus had the highest mean starting weight of 108.1 g. Although strains B (N = 250) and D (N = 862) were of comparable ages, Icelandic fish had a higher mean weight of 57.1 g. Starting weights of strains A and B, which originated from the first trial, had markedly higher standard deviations than fish in the second trial. While fish of strain D were around 70 d older than fish of strain C, both strains had similar mean weights. Table 2 shows means, standard deviations, range, and the coefficient of variation (CV) for all traits at the time of slaughter. In the complete dataset, mean BW was 765.5 g with a high standard deviation of 281.9 g. Mean FW was 227.5 g, which translates to a mean fillet yield (FY%) of 29.7 %. The adapted Fulton's condition index (FCI PLN ) had an average of 8.0, which was higher than the conventional FCI with 5.2. ADG between the ages of 300 and 500 dph was 1.9 g. All weight-based traits showed high variation (CV: 31.7-57.3), whereas coefficients of variation for morphometric traits and calculated ratios were remarkably lower (CV: 0.8-22.9). For strain-wise comparison BW, FW, FY%, and morphometric traits were in similar ranges. Differences in viscera traits and ratios were mainly based on a higher mean liver weight (+3.6 g) in Norwegian fish. Standard deviations, and accordingly CVs, were generally higher for all traits in Norwegian fish due to higher age ranges within this strain.

Heritability estimates of growth and carcass traits
Heritability estimates in this study were low to moderate, ranging from 0.04 to 0.29 in univariate runs (Table 3). Heritability estimates for the economically most important traits BW, FW, and FY% were low with 0.19, 0.18, and 0.15, respectively. Estimates of heritability for individual carcass elements GU, CA, IW, and HE were also low between 0.14 and 0.19, with the exception of LI, which had a moderate heritability of 0.24. For morphological production traits heritability estimates were within the same range (L: 0.18 and A: 0.12), apart from W with a heritability of 0.08. A comparison of FCI (0.25) and FCI PLN (0.04) reveals a lower estimate for the latter condition index. ADG between the ages of 300 and 500 dph showed the highest heritability (0.29) of all traits we examined. Table 2.
Descriptive statistics of all growth and carcass traits at slaughter for the complete dataset and separately for strains "Norway" and "France". BW: body weight at slaughter (g); L: body length (cm); W: body width (cm); A: body area (cm 2 ); FW: fillet weight (g), FY%: fillet yield (%), HE: head weight (g), LI: liver weight (g); GU: gut weight (g), HSI: hepatosomatic index; FCI PLN : adapted Fulton's condition index; DP%: dressing percentage (%); GW: gutted weight (g); IW: weight of intestines (g); CA: weight of carcass discard (g); ADG: average daily gain between the ages of 300 and 500 dph; FCI: Fulton's condition index; Min: minimum, SD: standard deviation; Max: maximum; CV: coefficient of variation. Separate analysis of strains using the same models generally resulted in higher heritability estimates, particularly for BW (0.27-0.31), FW (0.25-0.29), and FY% (0.14-0.23). Standard errors of all estimates were slightly higher than in the complete dataset. Estimates for biometrical traits were also higher in both separate strains, although trends and patterns between those traits remained the same with W (0.12-0.16) showing the lowest estimates and L (0.24-0.35), ADG (0.16-0.36), and DP% (0.12-0.43) showing the highest heritabilities. While estimates for FCI were similar between strains, heritability for FCI PLN was virtually zero in Norwegian fish and 0.12 in French fish, indicating that low and unreliable estimates for this trait in the complete dataset are based on one strain.

Genetic and phenotypic correlations between traits
Genetic and phenotypic correlations between BW (the target trait for selection) and all other traits are given in Table 4, while a complete table of all correlations for all traits can be found in the Supplement (Table S2). Phenotypic correlations between BW and other weight-based, as well as morphometric, traits were generally high (0.70-0.94). FW, HE, and CA had phenotypic correlation estimates higher than 0.9. Genetic correlations for these combinations were also above 0.9. Of the morphometric traits A and W had the highest positive correlations (> 0.9) with BW. While phenotypic correlations of both condition indices with BW were weak, we found moderate positive genetic correlation for BW-FCI PLN . Genetic correlations between FCI PLN and other traits were always higher than the ones between FCI and other traits. Genetic correlation between BW and FW was also high (0.95), while it was very low with large standard errors for fillet yield (0.11 ± 0.17) and dressing percentage (0.15 ± 0.15). Ratio traits expressed in percentages and both condition index traits often produced estimates difficult to interpret for both genetic and phenotypic correlations, since standard errors were often higher than the estimate itself. Additionally, negative estimates solely occurred for correlations involving those traits. The three biometric traits L, W, and A were highly genetically correlated with FW (r g : 0.86-0.94) and GW (r g : 0.89-0.99). Both genetic and phenotypic correlations for BW-AVG were high (r p : 0.79 and r g : 0.87). Due to model convergence problems correlations between traits BW-GW, FW-AVG, HE-CA, HE-W, and CA-AVG could not be estimated. Table 4. Estimates of phenotypic (r p ) and genetic (r g ) correlations for body weight at slaughter (BW) and other traits. BW: body weight at slaughter (g); L: body length (cm); W: body width (cm); A: body area (cm 2 ); FW: fillet weight (g), FY%: fillet yield (%), HE: head weight (g), LI: liver weight (g); GU: gut weight (g), HSI: hepatosomatic index; FCI PLN : adapted Fulton's condition index; DP%: dressing percentage (%); GW: gutted weight (g); IW: weight of intestines (g); CA: weight of carcass discard (g); ADG: average daily gain between the ages of 300 and 500 dph; FCI: Fulton's condition index; standard errors are given in brackets; NC: no convergence.

Body weight
Phenotypic Genetic (g) and the following correlations: r p correlations: r g

Discussion
The heritability estimates we obtained in this study are low to moderate, indicating potential for improvement via selective breeding programs. A number of studies on genetic parameters for body weight at harvest exist in turbot, which generally report slightly higher estimates. Sánchez-Molano et al. (2011) reported heritability estimates of 0.45 for body weight in their study. Xu et al. (2015) estimated heritability estimates ranging from 0.2 to 0.4 for body weight under different temperature regimes, whereas Gjerde (1997) reported higher estimates ranging from 0.4 to 0.7. We suspect that our results might be slightly underestimated due to both the use of molecular relatedness and the joint analysis of multiple strains. So far, there are only two publications using MR to estimate genetic parameters in flatfish: Blonk et al. (2010) calculated heritability for body weight in common sole (0.11-0.13), while Guan et al. (2015) compared the MR method to the classical pedigree-based approach and found heritability estimates for body weight at 100 dph of 0.19 (MR) and 0.66 (pedigree). Lower values for estimated variance components based on molecular relatedness have also been reported in a number of other species (e.g., Thomas et al., 2002;van Kleunen, 2005) and have been potentially linked to the denseness of the relationship matrix. Powell et al. (2008), among others, discussed a lack of comparabil-ity of different populations within the same species, which might partly be due to genetic differences in founder animals. Although the calculation of the genomic relationship matrix includes a standardization step to minimize founder effects, we received higher estimates in all traits when we analyzed strains separately, an effect that can be expected for within-family variance estimation. Since most turbot breeding programs in the different European countries treat their populations as separate from each other, these results could be considered more realistic. To our knowledge, this is the first time that genetic parameters for carcass and fillet traits are being reported for this species. Despite the obvious difficulties in direct comparability with other aquaculture species, particularly roundfish, it can be stated that our results seem to be within a reasonable range. Saillant et al. (2009) reported heritability estimates of 0.25 for FY% in sea bass, while Navarro et al. (2009) presented heritability of 0.15 for FW and 0.12 for FY% in gilthead sea bream. Neira et al. (2004) found very similar estimates (FW: 0.18, FY%: 0.11) in coho salmon, whereas heritability estimates in other salmonids are usually much higher (e.g., Powell et al., 2008; Atlantic salmon FW: 0.52). Many of the above studies reported problems in the analysis of mathematically derived ratio traits, which we experienced as well, especially regarding correlations involving those traits. Powell et al. (2008) and Rutten et al. (2005), in particular, argued the usefulness of yield traits as selection criteria since ratios are calculated from closely proportional traits. Consequently, we tried an alternative, indirect method to calculate heritability for fillet yield by including body weight as an extra covariable in the model for fillet weight, but we received close to identical results from both methods (results not shown here). Alternative methods to measure and quantify ratio traits are needed to make them more useful as selection target traits. Efforts have been made in other livestock species to develop direct, in vivo methods for fillet or carcass yield determination based on X-rays, computer tomography, or ultrasound (e.g., Sather et al., 1996;Milisits et al., 1999).
Genetic correlations between weight-based traits that were directly measured were all positive and high (> 0.7) (for example, BW at harvest, GW, or FW). We also found high positive correlations for the target traits BW and FW with morphometric traits. The highest correlations were found for A, which makes it a very suitable alternative selection trait. The use of morphometric traits such as length, width, or area as alternative selection criteria is particularly interesting since selection candidates do not have to be sacrificed to determine these measurements. With the advance of noninvasive, camera-based measurement techniques handling stress could be further minimized as well.
Another interesting selection trait could be ADG, since it is comparatively highly heritable (0.3) and strongly genetically correlated with body weight at harvest (0.88), but it can be measured nonlethally at significantly earlier life stages (300-500 dph in this case).
While a prior study (Arfsten et al., 2010) suggested the adapted condition index for Pleuronectiformes FCI PLN as a useful selection trait based on high phenotypic correlations with BW, we could not confirm this finding. When compared to the classical Fulton's condition index, the adapted FCI PLN showed higher genetic correlations with other traits. This resulted in lower h 2 estimates since a large part of the variance is already explained by the variance of the underlying traits. Like with other types of yield and ratio traits this makes them less suitable for direct selection.

Conclusions
This study presents heritability estimates and genetic correlations for weight and carcass traits in turbot. Estimates of heritability for weight-related traits were moderate, thus indicating that there is inter-individual variation and selection on these traits is beneficial. We were able to demonstrate the usefulness of biometrical body traits as alternatives for selection and found no potential for inadvertent negative selection response between body weight and carcass quality traits.
Data availability. Access to genetic and phenotypic data used in this study will be provided upon request to Jens Tetens (jens.tetens@uni-goettingen.de).
Author contributions. JT developed the research topic and acted as study principal investigator (PI). GT acted as co-PI. KS and VL collected phenotype data. KS collected genetic data, performed statistical analysis, and wrote the manuscript with support from JT, NK, and GT. NK provided valuable statistical expertise in the estimation of genetic parameters. CS provided useful expertise on animal rearing and husbandry.