Genetic Parameter Estimates for Growth of Hāpuku (Groper, Polyprion oxygeneios) in Land-Based Aquaculture Using Random Regression Models

Abstract

Hāpuku (Polyprion oxygeneios) is a promising candidate for aquaculture production in New Zealand. Methods for spawning, juvenile production, and growout to harvest entirely on land, where water quality, pathogens, environmental impacts, and genetic “pollution” can be tightly controlled, have been developed, and genetic improvement to optimise land-based production is the obvious next step. However, estimates of genetic parameters are required to design a rigorous, disciplined, and effective selective breeding program. By using existing data consisting of irregularly spaced repeated measurements of fork length and live body weight collected on wild-collected founders and two generations of captively reared progeny, we evaluated the species’ genetic potential for improvement in growth. We first tested a range of univariate random regression models to identify the best-fitting models for these data. Subsequently, using a bivariate model, we estimated variance components for growth trajectories of fork length and whole body weight. With one to six records available per fish, the best-fitting univariate models included only a fixed effect for contemporary groups and fixed and random genetic third-order Legendre polynomials. More complex models that included full-sib family and/or permanent environmental effects produced unacceptable constrained and/or non-positive-definite solutions. Both traits are moderately heritable at all stages of the growout phase (~0.4–0.5), and the genetic correlation patterns between daily breeding values estimated via the covariance function are different for length and weight. Genetic correlations for length between all pairs of age-specific breeding values are positive and strong (>0.7) and change gradually and smoothly with increasing temporal separation. For weight, these correlations deteriorate more rapidly with increasing time lags between measurements and become negative for some age pairings. We conclude that random regression analyses are a valuable tool for extracting genetic information from irregularly spaced repeated measurements of fish size, speculate that emerging technologies for high-throughput genotyping and phenotyping will add to the value of this approach in the near future, and reason that a breeding strategy that rigorously takes into account the potentially unfavourable genetic correlations between breeding values for weight at some ages will further adapt hāpuku to land-based systems and enhance the profitability commercial-scale production.

1. Introduction

Polyprion oxygeneios (Polyprionidae; Nelson, 1994 [1]), one of only two species in the genus, is generally known in New Zealand by its Māori name, hāpuku, or its common name, groper. The species is widely distributed in the southern hemisphere, occurring in southern Australia, southern Brazil, South Africa, Chile, and New Zealand at depths between 50 and 900 m [2,3,4], and has also been reported in the south-west Atlantic [5]. Hāpuku is a large, slow-growing species that can live for 60 years, reach up to 160 cm in total length, and weigh over 100 kg. Hāpuku reaches sexual maturity at around seven years of age [4,6], and growth is rapid during the first three years of life but slows considerably thereafter [7,8]. The species is highly valued by recreational and commercial fishers for its good flavour, flesh quality, and texture. However, catches in New Zealand are relatively low, with a total allowable commercial catch of ~1000–2000 metric tons annually under the New Zealand Fisheries quota management system [9,10].

There is global interest in Polyprion aquaculture [11,12,13], and hāpuku’s rapid early growth, high market value, and limited supply have attracted interest from the New Zealand aquaculture sector, which is striving to reach a government-set sales target of NZD 3 billion by 2035. Achieving this ambitious objective will require diversifying both the species farmed and the production systems used to grow them, with a focus on high-value species for domestic and export markets. While offshore sea cage farming is receiving a great deal of attention in New Zealand [14], the risks of escapees breeding with wild fish, the potential for detrimental environmental impacts from waste and damaged gear, and concerns over animal welfare, some countries are limiting or banning the practice in favour of closed pens and land-based production systems where biosecurity, pathogens, water quality, and waste products can be tightly controlled. To enable and de-risk the potential for land-based hāpuku production, New Zealand’s National Institute of Water and Atmospheric Research (NIWA) has, since 2002, pioneered reliable methods to condition and spawn captively reared broodstock, as well as land-based hatchery, juvenile, and growout methods (see [8] for a review).

To date, NIWA’s captive breeding and domestication program for hāpuku has generated over 100 G1 broodstock from wild-collected G0 parental stocks [8]. This, in conjunction with breakthroughs in spawning and larval-rearing techniques, provides a solid platform for genetically rigorous and disciplined selective breeding. From an economic perspective, growth rate and feed conversion efficiency are generally the highest priority traits for improvement in farmed fish species, at least initially [15]. Faster growth shortens the time to market and, in many fish species, leads to more efficient feed conversion [16,17,18,19,20,21], although there are notable exceptions [22].

Growth is a dynamically expressed quantitative trait [23] that varies temporally and responds to a range of environmental conditions such as temperature, water quality parameters, and feeding on a continuous time scale. Such traits are also often referred to as infinite-dimension traits that can be efficiently evaluated using mixed model random regression (RR) approaches to characterise their dynamics at the genetic level [24]. However, most fish breeding programs select for growth based on genetic variation in endpoint assessments, typically whole fish weight, filet yield and other production traits measured at harvest.

Furthermore, RR approaches can estimate the heritabilities and breeding values of growth traits at any age estimated via covariance functions and improve the accuracy of breeding values by exploiting more information [25,26]. Importantly, the random coefficients estimated using these models do not correspond directly to the genetic value for a specific phenotype, but rather the temporal pattern in predicted additive genetic value as a function of a continuous independent variable (typically age), and thus predict individual-level genetic trajectories of the trait(s) over time. In cultured fish, random regression analyses have been applied mostly to tilapia [27,28,29,30,31,32,33] but also to rainbow trout [34,35], Japanese flounder [36], and turbot [37]. These studies all conclude that this approach can be more efficient than measuring growth-related traits at specific ages. Furthermore, emerging technologies, in particular high-throughput phenomics, will soon be generating high-resolution time-series data for cultured fish amenable to random regression analysis (see [38] for a review).

The objective of this study is to explore the potential for implementing random regression models to estimate variance components and breeding values for growth traits in hāpuku raised in a land-based aquaculture system using pre-existing weight and length data collected at irregularly spaced points over an extended period during the initial phases of domestication of the species [39].

2. Materials and Methods

2.1. Experimental Population and Rearing Conditions

All fish rearing and data collection were conducted at the NIWA Northland Aquaculture Centre, Ruakākā, New Zealand (previously known as the Northland Marine Research Centre and Bream Bay Aquaculture Park when the data were collected).

We produced the experimental fish in three spawning seasons. Parents for the first two spawning seasons (2008 and 2010) were collected from the wild using hook and line from the east coast of northern New Zealand and from Cook Strait and spawned naturally in groups during their natural antipodal spring/summer spawning season (September–December). Parents for the 2013 spawning season were the captively reared progeny of the 2008 wild fish that had been raised in 20 to 70 m³ circular tanks with a nearly natural photoperiod and thermal regime (either two weeks or eight weeks advanced compared to ambient/natural). Sea water temperature ranged from a minimum of 10.0–12.5 °C prior to spawning to a maximum of 17–19 °C in summer. During spawning, broodstock was maintained at between 10 and 13.5 °C. Fertilised eggs were collected from the surface of the tanks by directing the outflow from the group spawning tank through a fine mesh screen suspended in an overflow box.

We incubated eggs and yolk-sac larvae in 400 L conical-bottomed tanks supplied with filtered, UV-disinfected, and oxygenated seawater at 14–19 °C. After the onset of exogenous feeding, we transferred larvae to circular nursery tanks (5 or 10 m³) and initially fed them enriched rotifers using a semi-static green water rearing protocol [40,41,42,43] followed by enriched Artemia. Following the transition to Artemia, water flows increased, and the addition of algae ceased. Rotifer and Artemia enrichments included Red Pepper (Bern Aqua, Olen Belgium) and Selco^® S.presso (Inve, Dendermonde, Belgium).

At ~50 days post-hatch (dph), the feeding larvae were weaned onto commercial dry feed (O-range [Inve, Dendermonde, Belgium]). Weaned juveniles were then transferred to 10 m³ flow-through tanks under ambient photothermal conditions and fed continuously to satiation on O-range (Inve, Dendermonde, Belgium) and/or NovaME (Skretting, Boxmeer, The Netherlands) commercial diets. We implanted the fish with passive integrated transponder (PIT) tags at ~60 g and periodically split/transferred the groups to larger tanks (20 to 70 m³) as they grew larger to avoid over-crowding, but the spawning-event specific progeny groups were never graded or mixed. Consequently, all fish measured from the same tank on the same date had experienced the same environmental and feeding conditions for their entire lives and were treated as contemporary groups in the genetic analyses (see Section 3.5 below).

2.2. Data Collection

Individual fork length and whole-body weight (hereafter referred to simply as length and weight) were collected from each fish starting from tagging. Before we measured them, we anesthetised all fish using either 2-phenoxyethanol at 200–400 ppm or Aqui-S at 20–40 ppm. Weight measurements were to the nearest g for all fish, and length to the nearest mm for fish under 300 g, and to the nearest 0.5 mm thereafter. We measured fish several times throughout their lives, but due to wider operational requirements at the site and other contemporaneous projects, these measurements were irregularly spaced, and the number of measurements varied between groups.

2.3. Genotyping and Parentage Assignment

Genotyping and parentage assignment protocols have been published previously [39], but briefly, we extracted DNA from fin clips collected from parents and progeny and preserved in ≥70% ethanol. Fin clips were then sent to AgResearch (Mosgiel, New Zealand) for processing, genotyping, and parentage assignment.

The marker panel consisted of nine multiplexed microsatellite markers. The PCR annealing temperature was 56 °C with a MgCl₂ concentration of 2.0 mM. The primer concentrations ranged from 0.1 µM to 0.6 µM (see Supplementary Table S1 for details). AgResearch analysed the amplification products using fluorescently labelled primers and an ABI3730 genetic analyser (Applied Biosystems, Foster City, CA, USA) They scored allele sizes using the GeneScan™ –500 LIZ^® Size Standard and GeneMapper Software v.3.7 (both Applied Biosystems, Foster City, CA, USA).

AgResearch also reconstructed pedigrees using a proprietary pedigree analysis program (K.G. Dodds, pers. comm.) that compares the DNA profiles of the progeny against all combinations of parents within the spawning tank. The probabilities and limits of detection were calculated to assist in the parentage analysis. Over 95% of progeny could be assigned to a unique parental pair.

2.4. Data Quality Control and Editing

We used the R statistical computing software version 4.3 [44] for all data processing, statistical analyses, and plots.

The raw dataset consisted of 25,103 records on 5343 individual fish. After they were older than 600 days post-hatching (dph), many of these fish were used in experiments aimed at inducing sexual maturity and spawning that involved hormonal and photoperiod manipulation [45,46]. Furthermore, hāpuku reached the target harvest size in 12–18 months. Consequently, we dropped all records with dph > 600 from the dataset. In addition, exploratory analyses revealed that the weight and length records for many of the fish weighing less than 50 g and 150 mm in length showed no within-assessment variation and were not individual measurements, but contemporary group means. We, therefore, also removed records for these small fish. Finally, we removed all records for contemporary groups with 10 or fewer records because these records were from intermittent (and in some cases lethal) sub-sampling rather than full assessments and contemporary group means based on such small sample sizes are unreliable. After these edits, the final dataset consisted of 16,733 records on 5252 fish derived from 16 spawning events involving a total of 27 sires and 20 dams. These spawns generated 108 unique full-sib (FS) families (

Table 1. Summary of spawning/hatching events that produced the evaluated fish. The number of sires and dams represented in the group, number of full-sib (FS) families, and number of individual progeny produced. Because some sires and dams contributed to multiple spawning events, potentially even mating with the same partners, the sires, dams, and family totals are not column sums.

Hatch Date	Sires	Dams	FS Families	Individuals
2008-09-04	7	2	14	1020
2008-09-21	7	3	14	676
2008-10-02	4	1	4	95
2008-10-22	5	3	8	125
2008-11-02	7	7	27	3908
2010-09-25	4	1	4	628
2010-09-26	5	3	11	1516
2010-10-06	6	7	20	1680
2010-10-14	5	6	18	598
2010-11-11	8	6	31	2346
2010-12-01	8	6	30	1639
2013-09-18	7	6	19	1672
2013-09-21	5	5	16	323
2013-09-26	1	1	1	6
2013-11-19	2	1	2	237
2013-11-28	5	1	5	264
Totals	27	20	108	5252

). We retained data from 79 of the original 129 contemporary groups/sampling events between 2009 and 2014 (Supplementary Table S2).

We used the R-package pedantics (https://github.com/cran/pedantics, accessed on 23 June 2024) to visualise the pedigree and summarise the relationships for the 5252 fish with phenotypic records and their wild (G0) parents without records. Figure 1 presents the individual-level pedigree. The vast majority of the measured fish were G1 progeny produced in the 2008 and 2010 spawns of G0 parents. Consequently, most of the information available is in the form of full and half-sib relationships within the G1 generation, although the data included some trans-generational relationships involving the G2 progeny produced in 2013. Supplementary Table S3 summarises the number of individuals in all relationship classes present in the pedigree.

2.5. Model Development and Genetic Parameter Estimation

We tested a total of 48 univariate random regression models (24 for each of the two growth traits) using ASReml-R v 4.2 [47] and the asremlPlus R- package (https://CRAN.R-project.org/package=asremlPlus, accessed on 20 June 2024). All 48 models included fixed Legendre polynomials representing the population mean growth curve [35] and a categorical fixed effect of the contemporary group (assessment × tank) to account for systematic differences in group means caused by shared differences in rearing conditions. All models also included random direct additive genetic Legendre polynomials with structure from the pedigree-derived numerator relationship matrix. Beyond these effects, we systematically varied the inclusion of (1) random categorical full-sib family effects and (2) random permanent environmental effects, which we attempted to estimate as either diagonal covariance structures that assume the correlations among the polynomials are zero or as unstructured matrices to estimate these correlations. The full-sib family effect combines maternal (genetic), non-additive genetic, and common environmental effects. Permanent environmental effects account for the residual covariance between repeated measurements. Finally, we varied the order of both the fixed and random Legendre polynomials from zero to three, with the orders of all polynomials included in a single model being equal. The general form of these models is as follows:

$$y=Xb+Z_{1}a+W_{1}fs+Q_{1}pe+\epsilon$$

(1)

in which $y$ is the vector of observations; $b$ is the vector of fixed effects; a is the vector of random coefficients for direct additive genetic effects; $fs$ is the vector of random coefficients for full-sib families; $pe$ is the vector of permanent environmental effects; ɛ is the vector of residual effects; and $X$, $Z_1$, $W_1$ and $Q_1$ are incidence matrices that associate individuals with specific levels of the fixed and random effects. This model is based on the following assumptions:

$$E=\left[\begin{array}{c}\begin{array}{c}y\\ a\\ fs\end{array}\\ pe\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{X}_b\\ 0\\ 0\end{array}\\ 0\end{array}\right];V=\left[\begin{array}{c}\begin{array}{c}a\\ fs\\ pe\end{array}\\ c\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{K}_a\otimes A\\ 0\\ 0\end{array}\\ 0\end{array}\begin{array}{c}\begin{array}{c}0\\ K_{fs}\otimes I_{N_A}\\ 0\end{array}\\ 0\end{array}\begin{array}{c}\begin{array}{c}0\\ 0\\ K_{pe}\otimes I_{N_A} \end{array}\\ 0\end{array}\begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ R\end{array}\right]$$

(2)

in which $E$ is the matrix of expectations for fixed and random effects; $V$ is the matrix of (co)variances; $K_a$ $K_{fs}$ and $K_{pe}$ are (co)variance matrices between random regression coefficients (of order equal to the order of the polynomial used) for direct additive genetic, full-sib family, and permanent environmental effects, respectively; $A$ is the additive genetic numerator relationship matrix; $I$ is an identity matrix; $N_a$is the number of animals for which records are available; ⊗ is the Kronecker product between matrices; and $R$ is a matrix of residual variances. Not all models included all the random terms, and as mentioned above, we used two different structures for $K_{pe}$ (diagonal and unstructured). The complete range of models evaluated is presented in

Table 2. Summary of the models evaluated. Columns represent the full range of possible fixed and random Legendre (Leg) polymonials. Rows represent specific models with the highest degree for each of the possible model terms. PE = permanent environmental effects; US = unstructured; DIAG = diagonal structure.

Model #	Degree of Fixed Leg	Degree of Genetic Leg	Degree of Family Leg	Degree of US PE Leg	Degree of DIAG PE Leg
1	0	0	-	-	-
2	1	1	-	-	-
3	2	2	-	-	-
4	3	3	-	-	-
5	0	0	0	-	-
6	1	1	1	-	-
7	2	2	2	-	-
8	3	3	3	-	-
9	0	0	0	-	-
10	1	1	1	-	-
11	2	2	2	-	-
12	3	3	3	-	-
13	0	0	0	-	-
14	1	1	1	-	-
15	2	2	2	-	-
16	3	3	3	-	-
17	0	0	-	0	-
18	1	1	-	1	-
19	2	2	-	2	-
20	3	3	-	3	-
21	0	0	-	-	0
22	1	1	-	-	1
23	2	2	-	-	2
24	3	3	-	-	3

Supplementary File S1 provides example R scripts for fitting the various models.

To compare the goodness of fit between models with different fixed effects (i.e., different orders of the fixed Legendre polynomials), we estimated the Aikeke information criterion (AIC) and Bayesian information criterion (BIC) for each model based on its full likelihood using the approach developed by Verbyla [48] and implemented as the infoCriteria() function in the asremlPlus package (https://CRAN.R-project.org/package=asremlPlus, accessed on 23 June 2024). Based on these criteria and on the constraints imposed on variance components by the ASReml-R solver caused by limitations in the dataset (see Section 3), we identified the preferred univariate model for each trait (lowest AIC and BIC; highest log-likelihood; no constrained or non-positive-definite components).

To more accurately estimate variance components and breeding values, we next fit a bivariate model using the same order Legendre polynomials as the preferred univariate models, again using ASReml-R v 4.2 [47]. After thus estimating the covariance matrices for the random Legendre polynomials, we estimated the genetic and phenotypic variances and heritabilities for weight and length at $dph$ranging from 50 to 400 at 10-day intervals based on the covariance function [49,50,51].

We used variations on R functions available at http://morotalab.org/UFV2019/day3/day3.html#9 (accessed on 10 June 2024) (Supplementary File S2) to transform the $dph$data to a −1 to 1 scale and calculate the normalised orthogonal Legendre polynomials according to the following expressions:

Standardised $dph$ at i sampling points:

$${dphs}_i=-1+2\left({\displaystyle \frac{{dph}_i-min\left(dph\right)}{\mathit{max}\left(dph\right)-min\left(dph\right)}}\right)$$

(3)

where max(dph) and min(dph) are the maximum and minimum values of $dph$ in the dataset.

Legendre polynomials (Φ) of order 0–3 calculated for all values of ${dphs}_i$ as follows:

$${\Phi }_k\left({dphs}_i\right)={\displaystyle \frac{1}{2^k}}\sqrt{{\displaystyle \frac{2k+1}{2}}}\sum _{m=0}^{\left[{\displaystyle \raisebox{1ex}{$k$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]}{(-1)}^m\left({\displaystyle \frac{k}{m}}\right)\left(\begin{array}{c}2k-2m\\ k\end{array}\right){\left({dphs}_i\right)}^{k-2m}$$

(4)

where ${dphs}_i$is any of the (standardised) $dph$ values at which the trait was measured, and $k$ is the order of the mth polynomial of order $k$.

We estimated the additive genetic (co)variance between any pair of transformed dph values as follows:

$$\begin{gathered} {\left(Co\right)V}_g={\sigma }_{a(i,j)}^2=Z_{1i} G Z_{i1} \\ {Z}_{1i}=\left[{\Phi }_0\left(i\right) {\Phi }_1\left(i\right) {\Phi }_2\left(i\right) {\Phi }_3\left(i\right)\right] \\ {Z\mathrm{\prime }}_{1j}=\left[\begin{array}{c}\begin{array}{c}{\Phi }_0\left(j\right)\\ {\Phi }_1\left(j\right)\\ {\Phi }_2\left(j\right)\end{array}\\ {\Phi }_3\left(j\right)\end{array}\right] \end{gathered}$$

(5)

where ${\Phi }_m\left(k\right)$ is the mth degree Legendre polynomial evaluated at dph = k. If $i=j,$ the result is the genetic variance for dph_i, and genetic covariances between dph_i and dph_j if $i\ne j$.

Because the bivariate model we fit estimated a single residual variance component for each trait, we separately calculated the heritabilities (h²) for length and weight at each dph as follows:

$$h^2\left(i\right)={\displaystyle \frac{{\sigma }_{a(i,i)}^2}{{\sigma }_{a(i,i)}^2+R}}$$

(6)

We further estimated the genetic correlations between all pairs of $dph$as follows:

$$r_g(i,j)={\displaystyle \frac{{\sigma }_{a(i,j)}^2}{\sqrt{{\sigma }_{a(i,i)}^2}\sqrt{{\sigma }_{a(j,j)}^2}}}$$

(7)

where i and j are indices for $dph$at all sampling points.

3. Results and Discussion

3.1. Mating Success and Parental Contributions

Figure 2 shows the numbers of progeny and the number of mates for each of the sires and dams that produced surviving progeny at tagging. The contributions of individual parents are highly heterogeneous (Figure 2a,b), with small numbers of sires and dams highly represented and the majority making lower genetic contributions to the progeny population. The number of mates per dam (Figure 2d) is more evenly distributed than mates per sire (Figure 2c), although a few dams mated with small numbers of sires. See Figure 1 for the pedigree.

While the observed level of variation in parental contributions and numbers of matings is less than ideal for a selective breeding program, it is manageable. The historical data we used for our analyses were generated using now outdated and expensive genotyping technology (nine multiplexed microsatellite markers and Sanger sequencing). Since these data were collected, higher throughput sequencing technologies have made it possible to characterise many more single nucleotide polymorphisms much more cost-effectively (see [52] for a recent review). These technologies include shotgun and targeted variations on genotyping-by-sequencing, low-coverage whole genome sequencing, and low-density SNP chips.

While developing alternative spawning strategies to balance parental contributions would be preferable, one potential approach for managing variable parental contributions in group mating species like hāpuku could be genotyping a large number of prospective selection candidates at an early stage using a low-cost genotyping approach to either reconstruct their pedigree or estimate all pairwise relationships using markers and using those data to assemble a smaller pool of selection candidates with more balanced contributions at tagging using, for example, optimal contributions approaches [53,54,55].

3.2. Phenotypic Data

A majority of the 5252 individual fish in the experimental population were measured between three and six times, but a substantial fraction was only measured once or twice (Figure 3). One of the advantages of random regression models is that they can make efficient use of repeated records data collected at irregular intervals with variable numbers of records for individuals. The historical data we analysed consisted of a limited and variable number of measurements on individuals, but we were nonetheless able to estimate genetic parameters (see below). Looking forward, however, as with genotyping, emerging high-throughput phenotyping methods, particularly those based on image analyses, are poised to revolutionise fish breeding programs. Several fish breeding programs have tested or implemented image-based phenotyping using out-of-water systems that drastically reduce the time required and costs of regular assessments, and in-water systems that collect dense time series data are under development (J. Bastiaansen, pers. comm.). Random regression approaches will undoubtedly be key to efficiently exploiting high-resolution time series data for genetic improvement.

Figure 4 presents the phenotypic growth trajectories for length (a) and weight (b) of the fish in the final dataset. Visually, the growth trajectory for weight is potentially slightly sigmoid/cubic, whereas the growth trajectory for length appears more quadratic. These overall patterns of growth suggest that relatively low-order polynomial regressions should be sufficient to model the genetic components of variation.

3.3. Genetic Model Selection

Table 3. Estimated additive genetic (Vg) and residual (Vr) variances, genetic correlations (Rg), and percent of additive genetic variance explained by the component (AV%) for Legendre coefficients (standard error) from a univariate (Model #4) and bivariate model with third-order fixed and genetic Legendre polynomials for length (A) and weight (B).

A. Length
Component	Univariate Estimate (SE)	AV%	Bivariate Estimate (SE)	AV%
Vg(L0)	589.102 (14.436)	83.1%	570.796 (13.938)	80.6%
Vg(L1)	99.761 (4.702)	14.1%	89.843 (4.302)	12.7%
Vg(L2)	40.571 (2.381)	5.7%	36.161 (2.18)	5.1%
Vg(L3)	15.598 (1.322)	2.2%	11.7 (1.154)	1.7%
Vr	17.316 (0.325)	-	19.683 (0.361)	-
Rg(L1,0)	0.57 (0.019)	-	0.603 (0.018)	-
Rg(L2,0)	−0.217 (0.03)	-	−0.249 (0.031)	-
Rg(L2,1)	0.425 (0.027)	-	0.45 (0.027)	-
Rg(L3,0)	0.243 (0.033)	-	0.233 (0.037)	-
Rg(L3,1)	0.05 (0.044)	-	0.147 (0.047)	-
Rg(L3,2)	0.506 (0.036)	-	0.503 (0.04)	-
B. Weight
Component	Univariate Estimate (SE)	AV%	Bivariate Estimate (SE)	AV%
Vg(L0)	50,451.476 (1339.018)	0.565	48,818.61 (1284.097)	0.547
Vg(L1)	33,364.535 (1091.894)	0.374	31,747.182 (1033.043)	0.356
Vg(L2)	8010.154 (333.67)	0.09	7615.405 (314.228)	0.085
Vg(L3)	1176.908 (65.241)	0.013	1088.831 (61.122)	0.012
Vr	233.312 (4.767)	-	268.081 (5.335)	-
Rg(L1,0)	0.926 (0.003)	-	0.928 (0.003)	-
Rg(L2,0)	0.378 (0.02)	-	0.366 (0.02)	-
Rg(L2,1)	0.691 (0.012)	-	0.681 (0.012)	-
Rg(L3,0)	0.115 (0.027)	-	0.094 (0.028)	-
Rg(L3,1)	0.368 (0.025)	-	0.37 (0.025)	-
Rg(L3,2)	0.814 (0.011)	-	0.834 (0.011)	-

summarises the various goodness of fit statistics for the 24 models we evaluated for each of the two traits. Model numbers and descriptions follow

Table 4. Summary of model fitting parameters for the 24 models evaluated for length (A) and weight (B). DF = degrees of freedom for fixed effects and variance components; US = unstructured. DIAG = diagonal; PE = permanent environmental. The “Constrained Boundary Effects” column is the number of variance components that ASREML-R coded as “B”, and the “Non Positive Definite Effects” column is the number of variance components coded as “?” (see text for further explanation). The rows for the 6 best-fitting models (smallest AIC and BIC, largest log-likelihood) are in bold italic font.

A. Length
Model #	Degree of Fixed Leg	Degree of Genetic Leg	Degree of Family Leg	Degree of US PE Leg	Degree of DIAG PE Leg	DF for Fixed Effects	DF for Variance Components	Constrained Boundary Effects	Non Positive Definite Effects	Full Model AIC	Full Model BIC	loglik
1	0	0	-	-	-	79	2	-	-	93,902.74	94,528.47	−46,870.37
2	1	1	-	-	-	80	4	-	-	89,498.55	90,147.46	−44,665.27
3	2	2	-	-	-	81	7	-	-	87,589.86	88,269.67	−43,706.93
4	3	3	-	-	-	82	11	-	-	87,044.04	87,762.48	−43,429.02
5	0	0	0	-	-	79	3	-	-	93,854.33	94,487.79	−46,845.17
6	1	1	1	-	-	80	7	-	-	89,458.45	90,130.54	−44,642.22
7	2	2	2	-	-	81	13	-	-	87,541.21	88,267.38	−43,676.61
8	3	3	3	-	-	82	21	-	20	86,986.52	87,782.21	−43,390.26
9	0	0	0	-	-	79	3	1	-	93,854.33	94,487.80	−46,845.17
10	1	1	1	-	-	80	8	1	-	89,460.03	90,139.84	−44,642.02
11	2	2	2	-	-	81	14	3	6	87,542.32	88,276.21	−43,676.16
12	3	3	3	-	-	82	22	2	20	86,986.22	87,789.63	−43,389.11
13	0	0	0	-	-	79	3	1	-	93,854.33	94,487.80	−46,845.17
14	1	1	1	-	-	80	10	-	3	89,451.59	90,146.85	−44,635.79
15	2	2	2	-	-	81	19	-	5	87,528.49	88,301.01	−43,664.25
16	3	3	3	-	-	82	31	-	30	87,033.82	87,906.76	−43,403.91
17	0	0	-	0	-	79	3	-	-	93,892.96	94,526.43	−46,864.48
18	1	1	-	1	-	80	7	-	3	89,484.60	90,156.68	−44,655.30
19	2	2	-	2	-	81	13	-	6	87,574.32	88,300.49	−43,693.16
20	3	3	-	3	-	82	21	-	10	86,995.70	87,791.39	−43,394.85
21	0	0	-	-	0	79	3	-	-	93,892.96	94,526.43	−46,864.48
22	1	1	-	-	1	80	5	1	-	89,491.82	90,148.46	−44,660.91
23	2	2	-	-	2	81	8	2	-	87,589.51	88,277.04	−43,705.75
24	3	3	-	-	3	82	11	4	-	87,044.03	87,762.47	−43,429.02
B. Weight
Model #	Degree of Fixed Leg	Degree of Genetic Leg	Degree of Family Leg	Degree of US PE Leg	Degree of DIAG PE Leg	DF for Fixed Effects	DF for Variance Components	Constrained Boundary Effects	Non Positive Definite Effects	Full Model AIC	Full Model BIC	loglik
1	0	0	-	-	-	79	2	-	-	167,669.07	168,294.80	−83,753.53
2	1	1	-	-	-	80	4	-	-	151,932.96	152,581.88	−75,882.48
3	2	2	-	-	-	81	7	-	-	142,313.28	142,993.09	−71,068.64
4	3	3	-	-	-	82	11	-	-	140,589.00	141,307.44	−70,201.50
5	0	0	0	-	-	79	3	-	-	167,576.64	168,210.10	−83,706.32
6	1	1	1	-	-	80	7	-	-	151,867.45	152,539.54	−75,846.73
7	2	2	2	-	-	81	13	-	6	142,276.59	143,002.75	−71,044.29
8	3	3	3	-	-	82	11	10	-	141,635.83	142,354.27	−70,724.92
9	0	0	0	-	-	79	3	1	-	167,576.64	168,210.10	−83,706.32
10	1	1	1	-	-	80	7	2	-	151,867.46	152,539.54	−75,846.73
11	2	2	2	-	-	81	13	2	-	142,276.69	143,002.85	−71,044.34
12	3	3	3	-	-	69	11	14	-	141,366.51	141,984.46	−70,603.26
13	0	0	0	-	-	79	3	1	-	167,576.64	168,210.10	−83,706.32
14	1	1	1	-	-	80	10	-	3	151,869.24	152,564.50	−75,844.62
15	2	2	2	-	-	81	19	-	12	142,277.00	143,049.52	−71,038.50
16	3	3	3	-	-	82	21	20	10	140,610.18	141,405.87	−70,202.09
17	0	0	-	0	-	79	3	6	-	167,670.91	168,304.37	−83,753.46
18	1	1	-	1	-	80	7	-	3	151,922.12	152,594.20	−75,874.06
19	2	2	-	2	-	81	13	-	6	142,308.60	143,034.76	−71,060.30
20	3	3	-	3	-	82	21	-	10	140,557.07	141,352.76	−70,175.54
21	0	0	-	-	0	79	3	-	-	167,670.91	168,304.37	−83,753.46
22	1	1	-	-	1	80	4	2	-	151,932.96	152,581.88	−75,882.48
23	2	2	-	-	2	81	7	3	-	142,313.28	142,993.09	−71,068.64
24	3	3	-	-	3	82	11	4	-	140,589.01	141,307.45	−70,201.50

. “Constrained Boundary Effects” are variance components that the ASReml-R solver codes as “B” indicating that it prevented a solution falling outside the “boundary” of allowable values (e.g., negative variance components) by assigning them a small acceptable value and “Non Positive Definite Effects” are variance components the solver coded as “?”, indicating that the solution matrix is not positive-definite and, therefore, must be rejected. This is most likely due to inherent limitations of the data, specifically insufficient information to estimate some components and/or multicollinearity between them.

The AIC and BIC presented are based on the full REML likelihood method [48], which is appropriate for comparing models with different fixed effects, such as different orders of the fixed polynomials. For both traits, the rows for the six best-fitting models (smallest AIC and BIC, largest log-likelihood) are shown in bold italic font (Models # 4, 8, 12, 16, 20, and 24). All of these models use third-order Legendre polynomials for the fixed and random regression effects (although the specifics of these effects vary between models). However, all of these models, except for Model # 4, must be rejected because at least some of the estimated (co)variances are constrained at the boundary and/or non-positive-definite. Based on this, we chose Model # 4, which uses third-order fixed and genetic Legendre polynomials as the best model for both traits and used a bivariate version of this model for the rest of our analyses.

While the RR approach was not able to fit all of the potential models without encountering data limitations that resulted in boundary constraints and/or non-positive-definite solutions (particularly models that included full-sib family and/or permanent environmental effects with polynomial orders > 0), it was able to fit simpler models with first- to third-order Legendre polynomials for fixed and direct additive genetic effects, producing estimates that aligned well with our expectations based on the overall shapes of the phenotypic growth curves (Figure 4). Excluding the potential non-additive genetic (full-sib families) and permanent environmental effects may overestimate the additive genetic (co)variances and heritabilities.

3.4. Genetic Parameter Estimates

Table 3 summarises the estimated genetic parameters for the Legendre coefficients for both traits based on both the third-order univariate and bivariate analyses. The first pattern to note is that the univariate and bivariate models produce similar genetic variance component estimates (Vg) for the polynomial coefficients for both traits, but the genetic variance estimates from the bivariate models are slightly smaller than those from the univariate models and have somewhat smaller standard errors (SE). This result is expected as the bivariate model accounts for the genetic correlations between the two traits, and the univariate models do not.

Most of the genetic variation (AV% for Vg) in the growth curves for both length and weight resides in the intercept term of the polynomials (L0), but the proportion of genetic variance explained by the linear term (L1) is larger for length than for weight (80.6 vs. 54.7%). The proportion of genetic variance accounted for by the quadratic term for length is nearly one-third that for weight (12.7 vs. 35.6%), and for both traits, the cubic terms account for even less variance (<10%). Furthermore, the intercept and linear terms show positive genetic correlations (Rg) that are moderate for length and very large for weight.

These results align well with our expectations based on the phenotypic growth curves (Figure 4), which are more linear for length than for weight. This implies that the growth curve for weight is more amenable to selection to change its shape because the higher-order polynomials that contribute to non-linearity are more genetically variable for weight than length.

Figure 5 plots the estimated genetic variances derived from the covariance functions for weight and length over the growout period analysed. For both traits, there is a consistent pattern of increasing genetic variance over time. Similarly, Figure 6 plots the estimated heritabilities for weight and length over the growout period. Both traits are moderately heritable, with the heritability for weight being slightly higher and more stable across the growout period than the heritability for length. These temporal patterns of genetic variation imply that selection on either length or weight at any age would be effective but do not address how selection at any specific age would affect other ages (see below).

Figure 7 shows the daily breeding values for 200 randomly chosen fish. These plots illustrate two important patterns: (1) individual fish re-rank extensively over time, particularly after approximately 300 dph, and (2) the individual-level temporal patterns in breeding values show substantial variation in both their elevations and shapes, with weight showing more variation in shape than length. This is a visualisation of the differences in the AV% for the degree 0, 1, and 2 Legendre polynomials for the two traits (Table 3). Together, these two patterns suggest that selection at different ages could result in genetic improvement not only in the two growth traits at specific ages but also in the shape of the mean growth curve in the post-selection population. This has important on-farm implications, such as how to optimise feeding rates over time and how often groups of fish must be split to avoid excessive densities.

Figure 8 presents the genetic correlations between all pairs of measurement dates from 50 to 400 days at 50-day intervals and reinforces these visual patterns quantitatively. For length (Figure 8a), the genetic correlations between all pairs of dph are positive and relatively strong (>0.7) and change gradually and smoothly with distance from the diagonal (i.e., increasing time separation between the two dph values). For weight (Figure 8b), the pattern is less regular. The correlations deteriorate more rapidly with distance from the diagonal and even become negative between some pairings of early (50, 100) and later (250–350) dph. This result differs from those of Zhao et al. [36] for turbot and of He et al. [27] for tilapia, all of whom found similar temporal patterns in time-lagged genetic correlations for weight and length and no evidence of negative correlations for either trait.

These temporal patterns in the genetic correlations indicate that while direct selection on length at any age will produce a reasonably large correlated response at all other ages, this is not the case for weight. Based on these estimates, the selection of weight early in the growth cycle will result in a negative response at 250–350 dph, the expected harvest time/size.

3.5. Synthesis and Prospectus

Taken as a whole, our analyses demonstrate that random regression models are a highly useful tool for estimating genetic parameters efficiently from even relatively sparse repeated records collected on multiple spawning cohorts at irregular intervals and that genetic improvement of growth in hāpuku farmed in land-based systems is likely to be highly effective. We could only fit simple models that ignored potential full-sib family (non-additive and maternal genetic) and permanent environmental variance components. Despite these limitations, our analyses support the hypothesis that the genetic variances and heritabilities (0.46–0.5) for length and weight are sufficient to support a response to selection at all stages of the growth cycle. It should be noted, however, that excluding the potential non-additive genetic and permanent environmental effects has likely led to some overestimation of the additive genetic (co)variances and heritabilities.

We also found substantial genetic variation in the shapes of temporal patterns in breeding values (i.e., the second- and third-degree polynomial coefficients; Table 3, Figure 7), particularly for weight. This indicates that it should be possible to select not only size at specific time points but also on the path by which fish achieve that size by directly selecting the breeding values of the higher-order Legendre regression coefficients.

In addition, random regression models allowed us to estimate the genetic correlations between measurements taken at different points in the growth cycle (Figure 8), a temporal analogue of the norm of reaction in genotype-by-environment interaction studies, and an important consideration for designing a disciplined and rigorous genetic improvement strategy. Interestingly, the patterns of these correlations differ markedly for the two traits. For both traits, time points separated by 50–100 days showed moderate to high genetic correlations, mostly > 0.9 for length and > 0.55 for weight. The high genetic correlations for length suggest that selecting length at any stage will produce favourable correlated responses at all stages. However, more caution will be required for selection on weight due to the more rapid decay of the genetic correlation with temporal distance and even negative genetic correlations between early and late measurements.

Finally, looking forward, although collecting repeated measurement data for cultured fish, such as the data we analysed, is currently costly, time-consuming, and stressful for the animals, this is changing rapidly. High-throughput, automated morphometric, health status, and deformity phenotyping using image analysis and other technologies are developing rapidly e.g., [56,57,58,59,60,61,62,63,64,65]). This emerging technology will undoubtedly become standard practice in the very near future, and in-cage and in-tank systems will likely produce irregularly spaced time-series data, making random regression methods an invaluable technique in fish breeders’ toolkits.

4. Conclusions

Random regression methods are seldom used but are powerful tools for estimating genetic parameters and breeding values using irregularly spaced repeated measurements of fish growth. These methods efficiently exploit all available information and make it possible to estimate key parameters for selective breeding at any age via covariance functions. In hāpuku (Polyprion oxygeneios), both length and weight are highly heritable at all phases of growout, but the genetic correlation structures between phases are different. All time-lagged genetic correlations are high and positive for length, but for weight, early and late growth is weakly or negatively correlated. Furthermore, genetic variation is second- and third-degree random regression coefficients that reveal genetic variation in the shape of growth trajectories. These methods make it possible to select the pattern of growth in addition to harvest weight in a disciplined and rigorous genetic improvement program based on repeated measurements. While acquiring such data is currently laborious and expensive, recent and near-future developments in high-throughput phenotyping, in particular automated imaging technologies, are potential game-changers that will make random regression analyses an increasingly important tool for fish breeders.