Genetic variance for uniformity of body weight in lumpfish (Cyclopterus lumpus) used a double hierarchical generalized linear model

Lumpfish (Cyclopterus lumpus) is the most widely use biological pest control agent to delouse sea lice on Atlantic salmon (Salmo salar). From 2012 to 2017, the production of lumpfish has increased from less than two million individuals to more than 30 million individuals in Norway alone. To meet the increasing market demand for lumpfish and to ensure sustainable production in the future, we need to domesticate the wild stocks and establish a breeding program for lumpfish. One of the most important traits in lumpfish, the genetic property for uniformity of body weight was investigated. In 2015, wild broodstock of lumpfish were collected from Southern and Northern coastal regions of Norway. A total of 68 full-sib families were produced and 7588 lumpfish juvenile were measured for body weight at tagging. A sire-dam double hierarchical generalized linear model was used to estimate the genetic parameters for body weight and its uniformity. Body weight was transformed by standardized and log scales. The heritabilities for uniformity of body weight in both scales, were low at 0.014 and 0.021 for the standardized and log scales, respectively. However, the genetic coefficient of variation for uniformity were relatively high, 45.8% and 63.5% (for standardized and log transformation), indicating that uniformity of body weight for lumpfish are substantially under genetic control. The genetic correlation between body weight and its uniformity was estimated. For the standardized transformation, the genetic correlation was −0.055 ± 0.252 (close to zero), and for the log transformation, the estimate was −0.813 ± 0.114. Both genetic correlations suggest that it is possible to simultaneously select for body weight and uniformity of body weight without a trade-off in genetic gain. In conclusion, there is substantial genetic variation for uniformity of body weight in lumpfish and the potential to select on this trait is promising.


Introduction
Cleaner fish as biological pest control agent for sea lice on Atlantic salmon (Salmo salar) has become one of the best alternatives to medicinal solution (Imsland et al., 2016a). Currently, lumpfish (Cyclopterus lumpus), the most suitable cleaner fish, is the third most valuable species in Norway, after Atlantic salmon (Salmo salar) and rainbow trout (Oncorhynchus mykiss) (Holmyard, 2018). In Norway, in terms of sale to Atlantic salmon and rainbow trout sea farms, lumpfish production has increased from 0.4 million individuals in 2012 to 30.3 million individuals in 2017 (Norway Directorate of Fisheries, 2018). With the increasing demand, domestication of wild stocks and establishing a breeding program for lumpfish are top priorities. These will allow the wild stocks to adapt to rearing environments, make it possible to select for economically important traits, and to minimize rate of inbreeding. However, before establishing a breeding program, information about important production traits and their genetic variations need to be understood.
The variation in growth rate for commercially produced lumpfish is usually large. Generally, the deployment size of lumpfish to commercial sea-cages is around 20 g (Imsland et al., 2016a). Imsland et al. (2016aImsland et al. ( , 2016b showed that lice eating ability in lumpfish decreases with increasing size of lumpfish. In accordance with this, it is a common view among salmon farmers that lumpfish above 400 g are inefficient lice eaters. Thus, it is important to produce more uniform stocks of lumpfish before deploying to Atlantic salmon cages. Large variation in body size T among individuals may also increase competitive interactions, contributing to dominant behavior and stress in subordinate fish, consequently reducing animal welfare (Ahvenharju and Ruohonen, 2007;Dou et al., 2004;Nielsen et al., 2014). From a selective breeding point of view, optimal growth (not too fast growth nor too slow growth) may improve sea-lice eating activity of lumpfish in Atlantic salmon cages.
Phenotypic variability is explained by environmental variance which indicates presence of micro-environmental sensitivity. Microenvironmental sensitivity can be defined as a change in environmental variance due to sensitivity to unknown environmental factors (i.e. unidentifiable factors that cannot be quantified or categorized) and subtle development (i.e. internal abnormal development that cannot be observed externally and will not cause phenotypic variation), which contributes to increased phenotypic variation, leading to increased size variation within a population (Sae-Lim et al., 2016). If genetic heterogeneity of environmental variance exists, selection for lower environmental variance will enable us to increase uniformity in growth of the fish. Heterogeneity of environmental variance has not previously been estimated for growth in lumpfish. Hence, the aim of this study was to quantify genetic parameters for body weight and uniformity of body weight in lumpfish.

Data
The data used in this study was generated from Center for Marine Aquaculture (CMAQ), Tromsø, Norway. In 2015, fertilized eggs from wild broodstocks were obtained from commercial lumpfish producers, Norsk Oppdrettservice AS and Akvaplan Niva. The broodfish used were caught by local fishermen from Southern (Flekkefjord, 58.2°N, 6.6°E) and Northern coastal regions (Sommarøy 69.3°N, 18.0°E) of Norway. A total of 23 sires and 60 dams were mated using a nested mating design to produce 68 full-sib families ( Table 1). Out of these 68 full-sib families, 11 were paternal half-sib and seven maternal half-sib (no overlapping of paternal and maternal half-sib families). The number of progeny produced by each full-sib family were ranged from 90 to 140 fish. After hatching, lumpfish fingerlings from the same full-sib family were separately reared in 125 l circular tanks using standard rearing protocols implemented at CMAQ until they reached the taggable size, which is at a mean body weight of approximately 20 g. Each lumpfish juvenile was tagged using passive integrated transponder (PIT) and was individually weighed using electronic balance Sartorius QC7DCE-S. The phenotypic data comprised 7588 records of individual lumpfish juvenile body weight at tagging.

Statistical analysis
For many morphological traits, a positive correlation between mean and variance, also known as scale effect, is expected. Previous studies indicated that data transformation, for example log-transformation, results in a trait variation that does not depend on the scale effect (Lande, 1977;Sonesson et al., 2013). Thus, in this study, body weight at tagging was transformed in two different ways, which were standardization and log transformation. In standardization (STDWT), individual body weight at tagging were standardized to a mean of 0 and variance of 1 to rescale the original data, which facilitates convergence . In the second transformation, individual body weight at tagging were log-transformed (LOGWT) to reduce dependency of variance on mean or the scale effect.
To estimate (co) variance components for body weight and its uniformity, a double hierarchical generalized linear model (DHGLM) was used (Rönnegård et al., 2010). Both STDWT and LOGWT data were modelled using the sire-dam DHGLM in ASReml (Felleki et al., 2012;Gilmour et al., 2009): where, y is the vector of either STDWT or LOGWT for the i th individual; Ψ is the vector of response variables for the residual variance, wherê= , which was linearized using a Taylor series approximation in ASReml (Felleki et al., 2012),ê i 2 is the squared residual of the i th STDWT, h i is the diagonal element in the hat-matrix of y (predicted value matrix) (Hoaglin and Welsch, 1978), andσ e 2 i is the predicted residual variance of the i th observation in the previous iteration of ASReml; X and X v are incidence matrices of the fixed effect (covariate of age at tagging) for the trait mean and its uniformity, respectively; b (b v ) is the solution vector for the corresponding fixed effect; Z s and Z d are incidence matrices for the random sire (s) and dam (d) effects; u (u v ) is the vector of additive genetic effects of sire-dam on either STDWT or LOGWT (uniformity), which was assumed to follow a normal distribution: where the 1/4 accounts for the fact that the sire and dam each explain only a quarter of the additive genetic variance for growth (σ2 a ) and its uniformity on the exponential scale (σ2 a v , exp ). The σ a, a v , exp is the additive genetic covariance between growth and its uniformity; Q (Q v ) is the incidence matrix for the random common effects to full-sibs; c (c v ) is the vector of common effects to full-sibs: where σ 2 c and σ 2 c v , exp are common environmental variances for growth and its uniformity and σ c, c v , exp is the common environmental covariance. The residuals of y(e) and Ψ(e v ) were assumed to be independently normally distributed as follows: , and σ 2 ϵ (σ 2 ϵ v ) is a scaled variance that was expected to be 1. The sire-dam DHGLM was fitted iteratively to update ψ, diag(W) and diag(W v ) until the log-likelihood converged (Felleki et al., 2012).

Calculation of genetic parameters
Estimates of σ2 u v , exp and σ2 c v , exp for uniformity of body weight were on the exponential scale (exp) and were converted to an additive scale (σ 2 u v and σ2 c v ) using the extension of the equations of Mulder et al. (2007), as derived by Sae-Lim et al. (2015). To reduce steps of calculations, the converting equations: where, σ4 E is the square residual variance (σ 2 E = σ 2 e − 2σ2 u ) from the animal model. The σ2 e is the residual variance of body weight from a sire-dam model. For the sire-dam DHGLM, the estimated genetic variance for sire-dam (σ2 u ) was equal to a quarter of additive genetic variance. Therefore, the additive genetic variance for body weight (σ2 a ) and its uniformity (σ2 a v ) were calculated as 4σ 2 u and 4σ2 u v , respectively. Phenotypic variance (σ2 P ) of body weight was equal to 2σ 2 u + σ 2 c + σ2 e , where σ2 c is the variance component for the effect common to full-sibs. Heritability estimate for body weight (h 2 ) was calculated as σ 2 a /σ2 P while heritability estimate for uniformity of body weight (h 2 v ) on the additive scale was calculated as . The GCV is the indicator of evolvability (a population's ability to response to natural or artificial selection; Houle, 1992), which expresses the biological opportunity for response to selection, relative to the mean of phenotype. Standard errors of h 2 v and GCV were approximated using the equations derived and presented by Mulder et al. (2016). The genetic correlation (r g ) between STDWT or LOGWT and its uniformity was calculated as σ a, a v , exp divided by the product of the two corresponding additive genetic standard deviations. Likewise, the common environmental correlation (r c ) was calculated as σ c, c v , exp divided by the product of the two corresponding common environmental standard deviations.

Descriptive statistics
The descriptive statistics of the traits is presented in Table 2. The average age at tagging was 163.3 days (minimum of 149 days and maximum of 185 days). The mean of body weight at tagging was 24.12 g with standard deviation of 8.96. The body weight at tagging by family is presented in boxplot as showed in Fig. 1, which has corrected for the differences in age at tagging.

Genetic variation of body weight and its uniformity
For STDWT, the σ 2 a and σ2 c were 0.514 and 0.045, respectively (Table 3). For LOGWT, the σ 2 a and σ2 c were reduced to 0.072 and 0.007, respectively. The h 2 estimate for STDWT was relatively high, 0.596 and with rather low c 2 , 0.050. The estimated h 2 and c 2 for LOGWT were almost the same as those for STDWT, which were 0.607 and 0.059, respectively.
The h 2 v estimates for uniformity either of STDWT or LOGWT were low, 0.014 and 0.021. Although, the σ2 a v for uniformity of LOGWT was lower than the estimate of STDWT, 0.001 vs 0.025, the GCV for uniformity of LOGWT was higher, 63.5%, than the GCV for uniformity of STDWT, 45.8%. The c 2 v for both uniformities of STDWT and LOGWT were low, 0.007 and 0.003, respectively.

Genetic correlation between body weight and its uniformity
The estimate of r g between STDWT and its uniformity was zero, −0.055 ± 0.252. However, the estimate of r g between LOGWT and its uniformity was strongly negative, −0.813 ± 0.114 (Table 3).

Genetic variation for uniformity of body weight
In this study, the estimated h 2 v for uniformity of STDWT and LOGWT for lumpfish was relatively low (0.014 and 0.021). However, the magnitude of these estimates are in line with previous studies reported in Atlantic salmon (h 2 v = 0.036 by Sae-Lim et al., 2017), rainbow trout (h 2 v = 0.024 by Janhunen et al., 2012; h 2 v = 0.011 by  and in some of the livestock species (h v 2 =0.028, or ranged from 0.006 to 0.047, as reviewed by Hill and Mulder, 2010). Based on our knowledge and reported studies, in livestock and aquaculture, the h 2 v estimates are at the maximum of 5% (for example, Hill and Mulder, 2010;Janhunen et al., 2012;Sae-Lim et al., 2017).
The low magnitude of h 2 v for both STDWT and LOGWT indicates that the accuracy of selection may be low, and large number of sibs is required to obtain sufficient accuracy of selection. In this study, the accuracy of sib selection with average h 2 v (0.0175) and c 2 v (0.005) calculated from estimates of STDWT and LOGWT is 0.440, using the average full-sib family size of 111 (Table 1). When full-sib family size approaches infinity, the maximum accuracy is 0.564. Due to high fecundity of lumpfish (Gregory and Daborn, 1982), it may be possible to increase accuracy by increasing family size. An alternative approach to increase accuracy is by using genomic selection, especially for the lowly heritable traits, such as uniformity of birth weight of piglet (Sell-Kubiak et al., 2015) and somatic cell score in dairy cattle (Mulder et al., 2013). Sae-Lim et al. (2017) reported that single-step genomic evaluation with an animal DHGLM can increase predictive ability of estimated genomic breeding values for uniformity (41% to 78%) compared to genetic evaluation based on numerator relationships. To our knowledge, there is no single nucleotide polymorphism (SNP) chip available for lumpfish. Thus, currently, to increase the accuracy of selection for lumpfish breeding program, we will have to increase the family size ().
To gain more understanding about the expected genetic change in uniformity of body weight, we calculated expected change using, ΔG = ir IH GCV (as a percent of the trait mean). With 10% of selected animals (selection intensity, i = −1.755), the accuracy (r IH ) of sib selection of 0.440, and a GCV of 0.458, the expected genetic change in residual variance of body weight was −35.4% per generation. The genetic gain may be lower if the number of full-sibs is reduced to, for example, 20 fish per family (−21.2%). The predicted reduction in Table 2 Number of observations (N), simple mean (μ), minimum, maximum, standard deviation (SD) and coefficient variation (CV, %) of age at tagging (days), body weight at tagging (g) in original, standardized (STDWT) and log-transformed (LOGWT) scales. residual variance in lumpfish is higher than the predicted reduction in rainbow trout , which could be mainly due to higher GCV reported in this study.

Data transformation
The σ 2 a v for uniformity of LOGWT (which has accounted for the scale effect) was more than 50% lower than the σ 2 a v for uniformity of STDWT (0.001 vs 0.025). However, the h 2 v for uniformity of LOGWT was slightly higher than the h 2 v for uniformity of STDWT. Our finding is in agreement with the rainbow trout study by . Authors reported a slightly higher h 2 v (0.024) estimated from logtransformed body weight in rainbow trout compared to the h 2 v estimate of non-transformed body weight (0.011). On the other hand, there are studies reported reduction in h v 2 after data transformation (e.g. Sae-Lim et al., 2017). The pattern of higher h 2 v after data transformation was observed at the breeding environment, while the opposite pattern was observed in the production environment . Hence, different environments may result in different direction of change of h 2 v after accounting for the scale effect.
The GCV presented in this study also correspond to the direction of change of h 2 v . It can be concluded that the potential for selective breeding is beyond the scale effects.

Genetic correlation between body weight and its uniformity
The estimated genetic correlation between body weight and its uniformity, under standardization scale, in lumpfish (−0.055) is lower than the estimates in Nile tilapia (0.60; Marjanovic et al., 2016), Atlantic salmon (0.952;Sae-Lim et al., 2017), and rainbow trout (0.79 and 0.30; . The genetic correlation between body weight and its uniformity was influenced by the scaling effect drastically. The genetic correlation between body weight and its uniformity was close to zero (−0.055). However, after log-transformation, the genetic correlation was strongly negative (−0.813). The change in magnitude of genetic correlation after data transformation is in line with the previous studies in aquaculture and livestock species (Yang et al., 2011). In Atlantic salmon, Sonesson et al. (2013) found that Pearson correlation between estimated breeding values of body weight and its uniformity changed from 0.42 to −0.17 after log-transformation. Similarly, the genetic correlation between body weight and its uniformity reduced from 0.952 to −0.093 after log-transformation in Atlantic salmon (Sae-Lim et al., 2017). In rainbow trout,  reported changes in genetic correlations between body weight and its uniformity after data transformation in the same direction as in this study, which were from 0.30 to −0.83 in breeding environment, and from 0.79 to −0.62 in production environment.
The uniformity of body weight may combine the micro-environmental sensitivity and scale effect while the uniformity of log-transformed body weight may consider only micro-environmental sensitivity. The definition of uniformity before and after data transformation has been discussed in previous studies Sae-Lim et al., 2017). For fish breeding programs, the major goal is to increase body weight and reduce micro-environmental sensitivity resulting in more uniformity fish. In other word, genetic control of observed residual variance is more relevant for selection on more uniform growth. The reason is the uniformity on the observed scale has a direct interpretation to an actual range of body sizes which are processed by aquaculture industries. Hence, the uniformity of observed body weight should be considered as a breeding goal trait. The decision whether the uniformity of log-transformed body weight may be included in the Fig. 1. Body weight at tagging (corrected for age at tagging) for each family in the lumpfish breeding program.

Table 3
Phenotypic and genetic parameters with their standard errors (in parenthesis) for body weight at tagging and its uniformity, under standardization and logarithm transformation. selection index may depend on the genetic correlation between uniformity of observed and log-transformed weights.

Conclusion
Our study shows that the genetic potential for response to selection for uniformity of body weight in lumpfish is substantial, but the accuracy of selection may be low due to low heritability for uniformity. A larger number of relatives may be required to obtain sufficient accuracy of selection if the breeding goal is to select for more uniform lumpfish in term of body weight. The genetic correlations which were close to zero and strongly negative after log-transformation suggests that selection for both body weight and its uniformity should be possible without trade-off in genetic gains when including both traits into a selection index.

Declaration of Competing Interest
None.