Multivariate random regression analysis for body weight and main morphological traits in genetically improved farmed tilapia (Oreochromis niloticus)

Background Because of their high economic importance, growth traits in fish are under continuous improvement. For growth traits that are recorded at multiple time-points in life, the use of univariate and multivariate animal models is limited because of the variable and irregular timing of these measures. Thus, the univariate random regression model (RRM) was introduced for the genetic analysis of dynamic growth traits in fish breeding. Methods We used a multivariate random regression model (MRRM) to analyze genetic changes in growth traits recorded at multiple time-point of genetically-improved farmed tilapia. Legendre polynomials of different orders were applied to characterize the influences of fixed and random effects on growth trajectories. The final MRRM was determined by optimizing the univariate RRM for the analyzed traits separately via penalizing adaptively the likelihood statistical criterion, which is superior to both the Akaike information criterion and the Bayesian information criterion. Conclusions In the selected MRRM, the additive genetic effects were modeled by Legendre polynomials of three orders for body weight (BWE) and body length (BL) and of two orders for body depth (BD). By using the covariance functions of the MRRM, estimated heritabilities were between 0.086 and 0.628 for BWE, 0.155 and 0.556 for BL, and 0.056 and 0.607 for BD. Only heritabilities for BD measured from 60 to 140 days of age were consistently higher than those estimated by the univariate RRM. All genetic correlations between growth time-points exceeded 0.5 for either single or pairwise time-points. Moreover, correlations between early and late growth time-points were lower. Thus, for phenotypes that are measured repeatedly in aquaculture, an MRRM can enhance the efficiency of the comprehensive selection for BWE and the main morphological traits. Electronic supplementary material The online version of this article (10.1186/s12711-017-0357-7) contains supplementary material, which is available to authorized users.


Background
From an economical point of view, growth and developmental characters are the most important traits in farmed fish species, and persistent efforts are made to genetically improve these traits in fish breeding. Growth and developmental traits, such as body weight (BWE) and morphological traits, are measured at different times. Faster BWE growth shortens time to market and selection on morphological traits allows fish shape and size to be standardized. Growth and developmental traits are dynamic quantitative traits because they vary spatially and temporally [1]. They are also called infinitedimensional traits, since they are expressed on a continuous time and space scale [2]. In tilapia breeding, growth and developmental traits at a specific age (in day, week, or year), such as the time to market are the main targets of genetic improvement. As such, genetic analysis of the growth and developmental process can increase the efficiency of selection compared to that of traits measured at specific ages [3]. He et al. Genet Sel Evol (2017) 49:80 Considerable attention has been paid to the genetic analysis of BWE and morphological traits at a specific age in breeding tilapia. In earlier studies, growth traits at specific ages were considered as separate traits and analyzed by using univariate animal models [4][5][6][7][8][9][10][11][12][13]. Subsequently, multivariate animal models were applied to estimate both the heritabilities of growth traits at specific ages and the genetic correlations between traits measured at different specific ages [8,11,12,[14][15][16][17][18][19][20][21][22]. Although multivariate genetic analysis of growth traits has improved the accuracy of parameter estimation by using more records, at multiple ages, such an analysis is strongly limited by the variable and irregular timing of these measurements [23].
For growth and developmental traits that are recorded at multiple time-points during growth, the patterns of growth curves have been shown to be heritable [2,3,24]. The genetic analysis of dynamic quantitative traits was initially conducted by first fitting individual growth curves and then analyzing their estimated parameters within the framework of a multivariate animal model. However, with such genetic analyses, it was not possible to determine whether some individuals had too few records to fit the growth curves. Since then, random regression models (RRM) [3,25] were developed to dissect dynamic phenotypes into different functions that describe the effects of various genetic and environmental factors. Only three papers on the application of RRM for the genetic analysis of growth curves have been published, in rainbow trout [26] and tilapia [9,27]. The use of RRM allows heritabilities of growth traits at any age to be estimated, as well as genetic correlations between pairwise traits measured at different ages. Moreover, multiple trait RRM (MRRM) can estimate genetic correlations between pairs of traits measured at different ages and improve heritability estimates for each trait. This not only facilitates the genetic analysis of dynamic traits, but also improves the prediction of breeding values [3,28].
The objective of this study was to construct and implement a MRRM to simultaneously model the genetic changes in growth traits and estimate the heritability of each growth trait, as well as the genetic correlations between pairs of growth traits measured at specific ages, such as the time to market for a population of genetically improved farmed tilapia in order to formulate a selection criterion for simultaneously selecting BWE and morphological traits. For this purpose, BWE, body length (BL), and body depth (BD) were measured six times on 1451 fish from 45 mixed families of full and half-sibs.

Experimental population
The base population of 1800 one-month old tilapia fingerlings (sex ratio of 1:1) was imported from WorldFish, Malaysia, in 2006, and was derived from 60 families with complete pedigree. The fish were systematically selected for seven generations at the experimental station of the Freshwater Fisheries Research Center in Wuxi, China. Based on breeding values for BWE at 120 days of age estimated using an animal model, 40 broodstocks per family were chosen at each generation, with a sex ratio of 1:1. The selection intensity in males and females was ~ 5%. Broodstocks from different families were randomly mated with each other and 100 to 125 families were retained in each generation. There was one generation of selection per year. In the seventh generation, in 2014, 120 males and 120 females from different families (one male and one female from each family) were selected as experimental parents. In May of that year, each male and two females ready to spawn were maintained in a 1-m 3 fiberglass tank for one week. Next, 120 females with fertilized eggs in their mouths were separately placed into multiple hapas (1 m × 1 m × 1 m) in a large concrete pond (50 m × 7 m × 1 m) for a one-week incubation period and were then isolated from their progenies. At 50 days of age, 50 families, each consisting of no less than 1000 surviving progenies above 30 g, were used to construct the experimental population. At the same time, 32 progenies that were randomly chosen from each family were tagged with passive integrated transponder tags, and 1600 tagged fish were mixed in a larger concrete pond (50 m × 20 m × 1.8 m). Because of this, we did not include any effect of the rearing facilities on the measurements in the statistical model. In the subsequent experiment, all fish were fed on a standard commercial diet (crude protein: 28%, crude lipids: 4%, crude fiber: 15%, ash: 18%, total phosphorus: 1%, lysine: 1.2%) manufactured by the Feedstuff Incorporated company (Ningbo, China). Dissolved oxygen was maintained at 8.23 to 8.67 mg/L by air pumps and the water temperature naturally fluctuated between 24.3 and 26.7 °C during the entire experiment.
The tagged fish were weighed on electronic scales to measure BWE, while BL and BD were measured with calipers. Before each measurement, the fish were anesthetized with clove oil at a density of 100 mg/L. The fish were measured once every 15 days during the experiment, with the first recording at the time of tagging. A maximum of six records for each fish and each trait were available, since 4.6% of the individuals died before the end of the experiment. Although experimental fish were measured only six times, the traits were recorded on 20 different days of age because of different spawning times. A total of 7560 records were extracted for BWE, BL and BD.
To edit the dataset, fish with missing phenotypic records for any trait and no information on the parents were excluded and abnormal phenotypic values that were more than three standard from the phenotypic mean of the trait were removed. After editing, 1451 individuals with 7235 records remained, with each individual having at least two records per trait. Table 1 presents the descriptive statistics for the dataset for BWE, BL, and BD for the 20 age groups.
Pedigree data were collected by tracing back three generations, thus including 1604 individuals from 77 sires and 88 dams. Although maternal effects on growth traits may be large in the early growth period for tilapia, they were not considered in the analyses, because observations from only a single generation and from such a small number of dams without records render the maternal effects unidentifiable.

Random regression model
In the same setting of cultivation programs, the two sexes and six test days were considered as fixed effects, while family, additive genetic, and permanent environmental effects were considered as random effects. For the growth-related traits that were recorded at multiple time-points, the phenotype of the trait at t days of age was described by the following single-trait animal model: where y ijkl (t) is an observation at t days of age, µ j is the j th test day effect, b k (t) is the kth sex effect, f l (t) is the lth family effect, a i (t) is the additive genetic effect for the ith individual, pe i (t) is the permanent environmental effect for the ith individual, and e ijkl (t) is the residual error following a normal distribution with expectation 0 and variance σ 2 e . To analyze the growth curves of phenotypes that are measured repeatedly, Legendre polynomials are generally used to fit the changes in the fixed and random effects in Model (1) and are represented by: where m is the order, β k is the kth regression coefficient, and ψ k (t) is the kth covariate of the Legendre polynomial [2], defined by: By modeling changes with age in the effects of sex, family, additive genetics, and permanent environment of the Legendre polynomials with orders m 1 , m 2 , m 3 , and m 4 , respectively, Model (1) becomes the following RRM: In matrix notation, Model (2) can be written as: where

The model satisfied the following:
where A is the numerator relationship matrix and F , G, and P are the family, additive genetic, and permanent environmental covariance matrices for the random regression coefficients of the Legendre polynomials, respectively.
To estimate the genetic correlations between growth traits measured at multiple time-points, the RRM for those traits must be solved simultaneously by a multivariate genetic analysis. The covariance matrices in the MRRM were estimated by using restricted maximum likelihood (REML), as implemented in DMU Version 6 [29]. The starting values were set to 0 for each fixed effect, to the identity matrix for each random effect covariance matrix and to 1 for the residual variance. The convergence criterion for REML was set to 10 −6 . After estimating the covariance matrices, heritabilities for the traits at any age and the genetic correlations between traits measured at different ages were estimated by the covariance functions [30]. The covariances for a trait between the ith and jth days of age can be calculated as z i Fz T j for the family effect, w i Aw T j for the additive genetic effect, and s i PEs T j for the permanent environmental effect. Therefore, genetic correlations between the ith and jth days of age for a trait were estimated by , phenotypic correlations by , and heritability at the ith day of age by

Choice of model
When choosing a model for a longitudinal analysis, Stoica and Babu [31] observed that neither the Akaike information criterion nor the Bayesian information criterion had optimum properties in terms of consistency and efficiency, and therefore the authors introduced a new criterion, i.e. penalizing adaptively the likelihood (PAL) [31], that overcomes these issues. More recently, Corrales et al. [32] applied

Choice of MRRM
The best MRRM for BWE, BL, and BD was selected by separately optimizing the univariate RRM for growth traits. The Legendre polynomial for the effect of sex was generally modeled by the population's growth curve, chosen as the Legendre polynomial of three orders according to the highest goodness-of-fit for the three analyzed traits. The three random effects, i.e. family, additive genetic, and permanent environmental effects, was characterized by using Legendre polynomials of different orders that ranged from 0 to 3 [3]. The models were designated as LPm 1 m 2 m 3 , e.g., LP121 is a model with a Legendre polynomial of order 1 for the family effect, order 2 for the additive genetic effect, and order 1 for the permanent environmental effect. A total of 64 RRM for each analyzed trait were compared based on the PAL criterion to determine the best combination of orders of the Legendre polynomials for the three random effects.
When fitting RRM with constant permanent environmental effects, the variance estimates of these effects were statistically inferred as not significant, using student t statistics (estimate/standard error). In addition, RRM with dynamic permanent environmental effects converged poorly in DMU. Therefore, permanent environmental effects were not considered in competing models. In the remainder, for the nine competing models, model LP00, which had constant family and additive effects, was taken as the reference, because the −2 log (ML) for the RRM with dynamic additive genetic and family effects was consistently lower than for model LP00. Based on Eq. (4), Table 2 presents the calculated PAL values for the nine competing RRM. Among these, LP23 was chosen for BWE with Legendre polynomials of order 2 for the family effect and order 3 for the additive genetic effect, while LP13 and LP22 were selected for BL and BD, respectively.

Genetic parameters for growth traits
Based on the model choices as described above, the final MRRM was constructed to jointly analyze BWE, BL, and BD. Based on Model (2), the MRRM, in which Legendre polynomials of different orders were nested, was expressed as: where y 1 (t), y 2 (t), and y 3 (t) are the phenotypes of BWE, BL, and BD, respectively.

Covariance estimates
Estimates of the (co)variances (standard errors) of the random regression coefficients were for the family effects and the genetic additive effects as shown in Figs. 1 and 2, respectively.  Estimates of the (co)variances (standard errors) of the random regression coefficients were for the random errors as follows: For each estimated covariance matrix of the random regression coefficients (including β 0 ), the diagonal block matrices are the covariance matrices for BWE, BL, and BD, respectively, while the off-diagonal blocks are the covariance matrices between the three growth traits. All covariances of the random regression coefficients were inferred as significantly different from 0. Estimates of the genetic correlations between most of the random regression coefficients for the additive genetic effects were positive, and the few negative genetic correlations were mainly between the intercept and the cubic regression coefficients or between the linear and cubic regression coefficients. In addition, estimates of genetic covariances for the same order random regression coefficients between the pairwise growth traits were all positive and decreased with increasing order. of age, using the univariate and multivariate RRM. The estimated ratios of family variances to phenotypic variances for BWE ranged from 0.259 and 0.365 and were consistently lower than these ratios for BL and BD for the period between 60 and 140 days of age. The estimated heritability obtained with the MRRM ranged from 0.080 to 0.635 for BWE, from 0.173 to 0.562 for BL, and from 0.098 to 0.610 for BD for the period between 60 and 140 days of age. Heritability estimates obtained from the MRRM were slightly lower than those obtained from univariate RRM for BWE and BL, but consistently higher for BD. Figure 3 shows that the heritabilities estimated with the MRRM increased with age for all three growth traits. Heritabilities changed most rapidly for BD and most slowly for BL. Heritabilities for BL were consistently higher than those for the other two traits from 60 to 75 days of age and those for BWE were higher than those for the other two traits from 75 days of age onward.

Heritability estimates
Considering 120 days of age as the time to market, heritabilities were above moderate and similar among the three growth traits. Figure 4 presents a three-dimensional plot of estimates of genetic correlations between pairwise traits measured at different ages for each of the three growth traits from 60 to 140 days of age. Patterns of the genetic correlations over time were similar for BWE, BL, and BD. Genetic correlation estimates between measurements at adjacent days of age were close to 1 but decreased monotonously as the lag between the days of age increased. In addition, genetic correlations between measurements at early days of age were lower than those at later days of age. The lowest genetic correlations were 0.566 for BWE, 0.633 for BL, and 0.497 for BD between 60 and 140 days of age. Table 4 includes family and phenotypic correlations between the pairwise selected growth days for BWE. Similar to the genetic correlations, phenotypic correlations approached 0.8 for the period between 75 and 120 days of age. In the earlier growth period, genetic correlations were higher than the phenotypic correlations between the pairwise days of age. The minimum family correlations were found between growth traits at 60 and 140 days of age, which were 0.573, 0.702, and 0.664 for BWE, BL, and BD, respectively. Family and phenotypic correlations for BL and BD were similar to those for BWE (see Additional file 1: Tables S1 and S2).

Correlation estimates
Patterns of the genetic correlations between the pairwise growth traits were similar to those between the pairwise days of age for the same trait, but they showed lower ridges, because the genetic correlation of the same trait on the same day is 1 (Fig. 5). Genetic correlations between traits at the same age ranged from 0.89 at 140 days of age to 0.95 at 60 days of age between BWE and BL, from 0.89 to 0.96 between BWE and BD, and from 0.85 to 0.95 between BL and BD. In contrast, genetic correlations between traits at different days of age were lowest between BWE at 140 days of age and BL at 60 days of age (0.50), between BWE at 140 days of age and BD at 60 days of age (0.47), and between BL at 60 days of age and BD at 140 days of age (0.44). At harvest (120 days of age), genetic correlations were equal to 0.92 between BWE and BL, 0.94 between BWE and BD, and 0.89 between BL and BD.
Family and phenotypic correlations between the pairwise growth traits exhibited similar patterns as the genetic correlations but different magnitudes at the same compared days of age (Fig. 6). Family correlations were consistently higher than genetic correlations between BWE and BD and between BL and BD. For the same day of age, most phenotypic correlations of pairwise growth traits were lower than the corresponding genetic and family correlations (see Additional file 1: Tables S3, S4 and S5).

Discussion
In this study, we used an MRRM to estimate heritabilities of BWE, BL, and BD based on repeated measurements during growth period, as well as genetic correlations between pairwise growth traits at specific days of age. Legendre polynomials were chosen to characterize the influence of fixed and random effects, i.e. family, additive genetic, and permanent environmental effects, on growth curves. The best MRRM was established by separately optimizing univariate RRM for each growth trait, according to the PAL. Based on the final MRRM, we found that heritabilities increased with age for all traits and that they were slightly lower than those obtained by using univariate RRM. In addition, for each trait, genetic correlations between measurements decreased monotonously with increased lag in days of age. For measurements at the same days of age, estimates of heritabilities and genetic correlations were close to those previously reported for tilapia [4-6, 8, 10-16, 18, 20-22, 33, 34].
For the selected RRM, changes in the genetic effects with day of age were optimally modeled for all analyzed traits by using Legendre polynomials of three orders. In contrast, Rutten et al. [9] and Turra et al. [27] used polynomials of two orders to fit the fixed and random effects, without selecting the orders of the polynomials for the different fixed and random effects. The non-significant variance for the permanent environmental effects showed that they have no impact on growth traits. Ratios of family variance [9] and of permanent environmental variance [27] to phenotypic variance are largely underestimated in a population with multiple families, which may be caused by the collinearity between the family and permanent environmental effects in the RRM used.
Repeated measurements during growth are required to estimate changes in both the fixed and the random effects with age when applying the RRM to genetic analysis of growth traits. More longitudinal measurements Estimates of genetic correlations between traits measured at different days of age (ra) for the three growth traits per individual would help to model growth curves with a higher goodness-of-fit and to more robustly estimate genetic parameters for random regression. However, more longitudinal measurements not only increase the experimental costs, but also affect fish growth, especially when all experimental individuals are measured together. In our experiment, the simultaneous measurement of 1600 tilapia fish incurred high labor costs and a measuring frequency of once every 15 days influenced fish growth to some extent. In future trials, the experimental population could be divided into several separately reared subpopulations, and each subpopulation could be observed in turn over more growth time-points. The extension of such measurements to more families from multiple generations would also help to identify maternal effects on early growth in tilapia within the MRRM framework.

Conclusions
Using repeated records of growth duration from multiple families in one generation, first we introduced MRRM to genetically analyze growth curves of body weight and main morphological traits in genetically improved farmed tilapia. The optimal MRRM was chosen by comparing the univariate RRM for the analyzed traits separately via PAL. By using the covariance functions of the MRRM, changes in heritabilities with days of age and genetic correlations between lags in days of age were estimated, which could be used to carry out early selection for each trait analyzed. More importantly, for phenotypes that are measured repeatedly in aquaculture, an MRRM can enhance the efficiency of the comprehensive selection for BWE and the main morphological traits at a specific age. Age  60  65  70  75  80  85  90  95  100  105  110  115  120  125  130  135