Body-shape trajectories and their genetic variance component in Gilthead seabream (Sparus aurata L.)

The phenotype of juvenile fish is closely associated with the adult phenotype, thus consisting an important quality trait for reared fish stocks. In this study, we estimated the correlation between the juvenile and adult body-shape in Gilthead seabream, and examined the genetic basis of the ontogenetic trajectories. The body shape of 959 pit-tagged fish was periodically examined during the juvenile-to-adult period. Individual shape ontogenetic trajectories were studied in respect to the initial (juvenile) and final (adult) phenotypes, as well as to the rate that adult phenotype is attained (phenotypic integration rate). We found that the juvenile body-shape presented a rapid change up to 192.7 ± 1.9 mm standard length, followed by a phenotypically stable period (plateau). Depending on the shape component considered, body-shape correlations between juvenile and adult stages ranged from 0.22 to 0.76. Heritability estimates (h2) of the final phenotype ranged from 0.370 ± 0.077 to 0.511 ± 0.089, whereas h2 for the phenotypic integration rate was 0.173 ± 0.062. To our knowledge, this is the first study demonstrating that the variance of the ontogenetic trajectories has a substantial additive genetic component. Results are discussed in respect to their potential use in selective breeding programs of Gilthead seabream.

Body-shape is a significant quality trait of reared fish, especially in the case of species that are marketed as a whole [1][2][3][4] . Except of its normal intra-species variation, under rearing conditions, body-shape is frequently subjected to the effects of skeletal abnormalities, which develop mostly during the hatchery phase 5 . Reasonably, juvenile phenotype has been considered as a good predictor of fish morphological quality at harvesting, with commercial hatcheries adopting certain practices to sort out the abnormal individuals, before they are transferred in sea cages 5,6 . After their transfer in cage farms, normal juveniles have rarely recorded to develop severe skeletal abnormalities 7 , or to recover from an existing skeletal abnormality 8 .
As every phenotypic trait, fish body-shape is determined by the action of both the environment and genotype. Abiotic parameters and nutrition can exert a direct effect on fish body-shape at given developmental periods [9][10][11][12] , or modify the ontogenetic rate of body-shape [13][14][15] , without necessarily leading to long-lasting phenotypic changes. In a similar way, genotype has been proven as a significant source of body-shape variation in reared fish, with the majority of existing studies targeting the near harvesting stages to improve the final product 1,2,4,[16][17][18] .
In most finfish species, body-shape ontogeny does not follow a uniform pattern, but presents different periods of rapid allometric changes, usually ending to a period of rather isometric growth [19][20][21] . Despite our interest in predicting the final phenotype of reared fish, to our knowledge, no studies exist on the correlation of fish bodyshape between different developmental stages, at the individual level. Relevant literature on other organisms (e.g., humans) demonstrates significant phenotypic correlations between different ontogenetic stages and a substantial variability in individual body-shape trajectories 22 . Furthermore, other studies show that the ontogenetic trajectories of various phenotypic traits (pigmentation pattern in snakes, herbivore defensive traits in plants, timing of developmental events in pond snails) may be heritable [23][24][25] .
In the present study, we examined whether juvenile body-shape is a good predictor of the adult body-shape in Gilthead seabream (Sparus aurata L.). Individual ontogenetic trajectories of body-shape were quantified for 959 fish and their genetic parameters were estimated. The studied organism is an important species of European Body shape ontogeny during the on-growing period. As it was demonstrated by the graph of Procrustes distances on SL, there was a significant size-effect on seabream body shape during the on-growing period (Fig. 4). In the SL range between ca 70 and 192.7 ± 1.9 (± SE) mm SL the PDs increased with the growth of the fish, whereas in the following ontogenetic period it was independent of SL ( Fig. 4, Table S1). Following relative warp analysis, size-effects on body shape were evident for RW1 (43.4% of the total variance explained, Fig. 5A). Simlarly to what was observed for PD, RW1 changed with fish growth up to 202.0 ± 1.7 (± SE) mm SL, whereas over this size RW1 was independent of SL ( Fig. 5A). Ontogenetic shape variation between juvenile and adult samples mainly concerned the anterior body parts. The transition from the juvenile to adult phenotype was char-  (B) Pearson's correlation coefficients between the Procrustes distances at the end of on-growing (434 dpt) with the Procrustes distances of the previous sampling points, up to the beginning of on-growing period (1 dpt). The average standard length (SL) of each sample is given. All estimated correlations were significant (P < 0.05). www.nature.com/scientificreports/ acterised by posterior shift of the snout and the anterior margin of the eye, as well as by an anterior transposition of the gill-cover, pelvic and pectoral-fin bases (Fig. 5A).
In opposite to RW1, no abrupt changes during fish growth were observed in the case of RW2 (12.4% of the total variance explained, Fig. 5B). Shape variation across RW2 axis mainly concerned the proximal shift of the dorsal profile (landmark 13), of the pectoral-fin bases (landmark 8) and of the anterior anal-fin base (landmark 7), as well as the ventral shift of the snoot (Fig. 5B). Similarly, to RW2, no size effects were evident in the case of the rest four relative warps (RW3-RW6, Fig. S1) that cumulatively explained 24.3% of the total body-shape variance.
Quantification of body-shape trajectories. To visualize the body-shape trajectory of each individual in respect to the consensus trajectory, shape data were categorized into twelve SL classes of 20 mm each (Table S2). Consensus allometric trajectories were plotted on the 5, 25, 50, 75 and 95 percentiles of each SL class (Fig. 6, Table S2). Interestingly, different body-shape (RW1) trajectories could be observed in different individual fish (Fig. 6).
To quantify the body-shape trajectory which was followed by each individual fish, the linear-regression parameters (slope, intercept) of SL-RW1 and SL-PD relationships were estimated separately for each specimen during the period of the first three samples (phase 1); i.e., when body-shape changes rapidly with the growth of SL (1-282 dpt, Fig. 6). Estimated coefficients of determination (r 2 ) were high for all the fish, both in the case of RW1 (0.34-1) and PD (0.13-1) (Table S3). For the period of the last three samples (phase 2), when RW1 and PD were not affected by SL (282-434 dpt, Fig. 6), the mean RW1 and PD were calculated separately for each individual, as a body-shape estimate at the trajectory plateau (pL-RW1, pL-PD, Table S3). Following the lack of  www.nature.com/scientificreports/ abrupt changes of RW2 with SL growth (Fig. 5B), the mean RW2 of each individual (av-RW2) was calculated to describe the plateau of SL-RW2 trajectory throughout the entire studied period (1-434 dpt).
Genotyping and parental assignment. From the 959 animals that were used to collect data, 913 (95%) were assigned to either one single pair or one single parent eligible to be used to quantitative further analysis. From these, the 84% of the offspring was assigned to a single pair whereas 11% to one parent but multiple pairs.  Genetic parameters of body-shape trajectories. Heritability estimates of all PD-SL trajectory parameters were statistically significant, ranging from 0.098 (Y-intercept, a-PD), to 0.173 (slope, b-PD) and 0.390 (plateau, pL-PD) ( Table 1). Similar results were observed for the RW1-SL trajectory parameters, with however insignificant h 2 estimate for the slope (b-RW1) ( Table 1). A high heritability estimate was present for av-RW2 (0.511, Table 1). Finally, h 2 estimates for fish standard length (SL) increased with fish age, ranging from 0.159 (1 dpt) to 0.291 (434 dpt, Table 1). Phenotypic and genetic correlations were estimated for all the examined traits, including those with insignificant heritabilities (Table 2). Generally, significant genotypic correlations were substantially fewer than the significant phenotypic. Negative high genetic (-0.717 to -0.800) and phenotypic (-0.897 to -0.905) correlations were found between the slope and Y-intercept of both SL-RW1 and SL-PD ontogenetic relationships examined. Phenotypic and genetic correlations between fish SL at different sampling ages increased as age difference Table 1. Estimates of heritabilities (h 2 ) for body-shape trajectory traits (RW1 and PD) and fish standard length (SL). P values indicate the significance of the estimates' difference from zero, ns not significant difference, dpt days post-tagging, PD Procrustes distance, RW1, RW2 scores of the 1st and 2nd relative warp axis respectively. 1, mean RW1 (282-434 dpt). 2, mean PD (282-434 dpt). www.nature.com/scientificreports/ decreased, ranging between 0.179-0.888 (phenotypic, r p ) and 0.273-0.948 (genetic, r g ). In the most of the rest trait pairs examined, phenotypic correlations were significant but low ( Table 2). The significant genetic correlations between body-shape traits and SL were low, with a maximum value of 0.461 (pL-PD with SL 282 ). No significant genetic correlations were observed between av-RW2 and SLs. Interestingly, significant genetic correlations found between av-RW2 and the plateau phenotypes of SL-PD (a-PD, 0.518) and SL-RW1 (a-RW1, -0.644) trajectories (Table 2). Concerning the rest trajectory parameters, a significant medium genetic correlation was observed between pl-PD and b-PD (0.672), as well as between pl-RW1 and a-RW1 (0.652) ( Table 2).

Discussion
The present study showed a significant but low correlation of body-shape between the juvenile and adult stage of Gilthead seabream. Body-shape ontogeny during the juvenile-to-adult period presented a substantial intrapopulation variation, in respect to the final (adult) phenotype (pl-PD), as well as the rate that this phenotype is attained (b-PD, slope or phenotypic integration rate). Significant heritability estimates for both of these traits (Table 1) showed that variation in shape ontogenetic trajectory has a significant genetic component. High heritability estimates for body-shape have been reported in a variety of fish species 1,26-28 . To our knowledge however, this is the first study documenting that variation in phenotypic integration rate has additive genetic component.
To follow the ontogeny of overall body-shape, in this study, Procrustes distances (PD) were estimated as an overall measure of multivariate shape components (i.e. relative warps, RWs). In agreement with previous studies 29 , PD successfully described the overall shape ontogeny with fish size (i.e., SL), additionally allowing the estimation of genetic parameters for body-shape trajectory. Compared with the relative warps however, PD showed smaller phenotypic correlation between juvenile and adult stage, and thus inferior capacity in predicting the adult from the juvenile body-shape. Relative warps analysis, analogous to Principal Components Analysis, replaces original shape variables with new ones (RWs) that are independent of each other, thus partitioning the overall shape variation into different components 30,31 . In the present study, significant differences were observed between the different RWs, not only in respect to the explained shape variation, but also in respect to their dependence on fish growth (allometry), correlation between the juvenile and adult stages and heritability estimates for final phenotypic scores (i.e. plateau phenotypes). Size-dependent shape variation (allometry) was expressed along the RW1 axis, whereas the rest RWs accounted mostly for size-independent shape variation. Such a partitioning is typical of relative warp analysis, especially when shape is studied over a relative wide size range 1,10,19,21 .
Body-shape correlations between different stages of fish ontogeny have been studied in the past, mostly with respect to the phenotypic evolution of abnormalities-related variation [32][33][34] . In Gilthead seabream, it was demonstrated that body shape is not altered after ca 70 mm total length, thus strengthening the hypothesis on the high predictability value of juvenile body shape 21 . In spite of their significance, in our study, shape correlations between different ages were rather low (0.22-0.76), even in the case of size-independent shape variables (i.e. RW2). Taking into account the environmental variation during fish growth in sea cages (e.g., water temperature 8 ), this result suggests that seabream body-shape is subjected to the effects of genetics and environment throughout the entire juvenile-to-adult period. In agreement to this hypothesis, previous studies revealed that on-growing environment (i.e., tanks vs sea cages) has a significant effect on seabream body-shape at harvesting 3 . In agreement to other studies 21 , our results revealed a clear breakpoint during the ontogeny of Gilthead seabream body-shape, at however a substantially bigger fish length (193-202 mm SL). The difference in the estimated breakpoints between the two studies, might be attributed to the inclusion of larval period and to its effects on the overall shape trajectory 21 . The breakpoint observed in the present study might be associated with the process of sexual maturation, which in Gilthead seabream takes place after ca 20 cm TL 35 .
In the present study, moderate to high heritabilities (0.370-0.511) were estimated for Gilthead seabream shape-related traits at harvest (i.e. pl-PD, pl-RW1, av-RW2), clearly suggesting that body-shape can be exploited in a selective breeding program. Among the examined traits, interest should focus especially on RW2 because it explains a part of morphological variation (i.e., proximal/distal shift of the dorsal, anal and pelvic fins) which is related with the consumers' preferences for the body-shape of reared Gilthead seabream 3 . This trait (RW2) has a high heritability (0.511) indicating that the accuracy and expected gain would be noticeable in a breeding program. Heritability of SL increases with age and this has been reported earlier for other species [36][37][38] since the larger the distance from the early stages the less the impact of non-additive effects (genetic or environmental) 38,39 . On the contrary, the correlation estimates decrease with the age indicating that the distant ages have lower correlations, either genetic or phenotypic, since we are referring to different ontogenetic stages and developmental mechanisms. The RW2 is slightly correlated (either phenotypically or genetically) with the other shape traits (RW, PD and SLs) indicating that selection on one trait will slightly affect the other ones. The same stands for most of the genetic or phenotypic associations among SL, PD and av-RW traits. There are no reported correlations estimates for such trait associations (SL with PD or RW) in the literature. However, the reported genetic or/and phenotypic correlations of body weight with several traits of shape coincide very well with our estimates since they range between 0.2-0.5 1,2,4,[16][17][18] .
During their ontogeny, animals follow rapid allometric changes to attain the adult phenotype. In organisms, like finfish, with external fertilization and planktonic early life stages, species-specific allometric growth patterns usually present multiple breakpoints between periods of different ontogenetic rates 21,40,41 . Environmentally-driven plasticity of ontogenetic trajectories may significantly alter the rate of shape ontogeny and the breakpoints (ontogenetic scaling), frequently leading to changes of body-shape at specific stages (i.e. juveniles, adults, etc.) 3,15,42 . Following the results of the present study, the rate of body-shape ontogeny is also subject to a genetically driven variation. In the future, it would be interesting to examine whether other organisms, besides fish, also present a strong genetic component with respect to their ontogenetic trajectory. www.nature.com/scientificreports/

Methods
Fish origin and samplings. During this study no experimentation with alive animals was performed. The examined biological material consisted exclusively of digital photographs, which were derived from our previous study on lordosis recovery 8 . Examined fish group consisted of 959 seabream juveniles (86 ± 7 mm standard length, SL) with normal external morphology, which were tagged electronically (FDX-B, Trovan Ltd, USA) and examined periodically for their body-shape at tagging (1 dpt, days post-tagging, 155 dph, days post-hatching), and at 77 (232 dph), 282 (437 dph), 371 (526 dph) and 434 dpt (end of on-growing) (589 dph) 8 . During samplings, all specimens were anaesthetised (ethyleneglycol-monophenylether, Merck, 0.2-0.5 mL·L −1 ), scanned for ID recognition, photographed on the left side and returned to the sea-cage. The digital photographs of the specimen were taken with a Canon PowerShot G9 camera, mounted on a tripod 8 . The examined fish population originated from a common larval population and egg batch, and reared according to the standard methodology for the species, in a land-based hatchery up to the juvenile stage and in a sea-cage farm for on-growing. The eggs were spontaneously spawned from a breeding population of 60 males and 59 females, which were kept under controlled photoperiod and temperature conditions. Fish rearing was performed at Andromeda S.A., a registered company (registration number GGN 5,200,700,699,992) for aquaculture production in Greece. The company is certified with GLOBAL G.A.P quality certification, with certified Veterinary doctor that verified the health and welfare of the fish. Animal sampling followed routine procedures and samples were collected by a qualified staff member from a standard production cycle. The legislation and measures implemented by the commercial producer complied with existing Greek (PD 56/2013) and EU (Directive 63/2010) legislation (protection of animals kept for farming). Production and sampling, by an experienced staff member, were optimised to avoid unnecessary pain, suffering or injury and to maximise fish survival.
Morphometric analysis. The ontogeny of body shape throughout the on-growing period was studied by geometric morphometry. On the digital photographs of the examined fish, thirteen homologous landmark measurements were taken by means of tpsDig2 software 43 (Fig. 7). A generalized least square method was applied 44 to adjust landmark configurations for centroid size and remove any irrelevant to shape variation. To compute the body-shape ontogenetic trajectories, the Thin-plate Spline algorithm was then applied and relative warp analysis (a principal component analysis of the partial warp scores 45 ) was carried out on the resultant weight matrix (TpsRelw software 45 ). To visualize the body-shape changes during the on-growing period, vector diagrams were obtained after the regression of shape components on the relative warp scores (tpsRegr software 46 ).
As a proper metric for shape dissimilarity in the Kendall shape space, the Procrustes Distances (PDs) were used to estimate the overall body-shape changes of each specimen during growth 31,47 . For PD calculation, a subtotal of the thirty smaller juveniles (71.1 ± 2.1 mm SL, 1 dpt) was chosen as reference point. The Procrustes distance (PD) of each individual from the reference point was calculated by the formula: where RWi is the individual score for each one of the estimated relative warps (RW1-RW22) and rRWi is the mean RWi for the reference group of fish. To examine whether the juvenile body-shape is a good predictor of fish body-shape at later stages, Pearson's correlation coefficient between the individual PD values of the first (1 dpt) and each one of the next four sampling points (77, 282, 371 and 434 dpt) was estimated. This analysis was also repeated for the RW scores, each of which explained a specific part of body shape variation.
The SL at which body shape trajectories presented a breakpoint was estimated by a piecewise linear regression (software SPSS was used 48 ), fitted with a non-linear estimation procedure  www.nature.com/scientificreports/ where Y is the studied body-shape variable (PD or RWs), SL is the standard length, b 0 is the y-intercept, b 1 is the slope of the Y-SL relationship during the first trajectory phase, b 2 is the change in b 1 that results in the slope of Y-SL relationship at the second trajectory phase, and L is the SL at the breakpoint 49 . To quantify the body-shape trajectory which was followed by each individual fish, the linear-regression parameters (slope, intercept) of the SL to shape-variables were estimated separately for each specimen and trajectory phase.
kit (Machery-Nagel, Germany), according to the manufacturer. DNA of the offspring was extracted using the Chelex-100 resin 50 . According to that, the tissue sample (of approximately 1 mm 2 ) was added to a 96-well plate and mixed with 100 μl of 10% Chelex-100 resin solution and 15 μl of proteinase K (10 mg/ml, Boehringer Mannheim). The plate was then incubated at 55 °C for one hour (shaken from time to time) followed by a 30-min incubation at 100 °C to deactivate the proteinase K and extent of protein denaturation. The DNA extracts were stored at 4 °C or at -20 °C. Before every use, the mixture was vortexed and centrifuged at 10,000 rpm for 10 min to separate the surface layer in which the DNA can be found and the lower layer which contains the Chelex-100 resin, the denatured proteins and other elements. 1-2 μl of supernatant are used for each 10 μl of final volume of the PCR reaction mixture. All parents and offspring were genotyped using a multiplex-PCR of nine microsatellite markers [51][52][53][54] . 10 μl volumes containing 0.4 unit of Taq polymerase (KAPA Biosystems), 1 × Taq buffer, 0.2 mM dNTPs mix, 1.5 mM MgCl 2 , 0.35 μM of forward and 5′-fluorescently labelled reverse primer and approximately 20 ng of template DNA were used to perform the multiplex PCRs. An initial three-minute 95 °C denaturation step was followed by 34 cycles of 30 s at 95 °C, 30 s at 53 °C and 30 s at 72 °C, with a final extension at 72 °C for 20 min. An ABI PRISM 3500 DNA Analyzer (Applied Biosystems), was used to separate fluorescently labeled PCR products, with a GeneScan 500 LIZ Size Standard (ABI) internal size standard. GeneMapper (Applied Biosystems) software was used to size alleles and genotype individuals. The presence of genotyping errors or null alleles was checked with the software Micro-Checker 55 . Parentage assignment was performed by VITASSIGN software allowing two mismatch allele 56 .

Data analysis.
Heritabilities and phenotypic correlations were calculated using phenotypic data collected on 911 animals. Heritabilities were analyzed for all data using WOMBAT 57 . An animal model was fitted: where Y is the vector of observations, β is the vector of fixed effects (overall mean), u is the vector of random additive genetic effects and e is the vector of random residual effects. X, Z are known incidence matrices. The genetic correlations were estimated for the traits applying a multivariate model.
Ethical statement. During this study no experimentation with alive animals was performed. The examined biological material consisted exclusively of digital photographs, which were derived from our previous study on lordosis recovery 8 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.