Genetic parameter estimates for methane emission during lactation from breath and potential inaccuracies in reliabilities assuming a repeatability versus random regression model

Methane (CH 4 ) emissions will be added to many national ruminant breeding programs in the coming years. Little is known about the covariance structure of CH 4 traits over a lactation, which is important for optimizing recording strategies and to establish optimal genetic evaluation models. Our aim was to study CH 4 over a lactation using random regression (RR) models, and to compare the accuracy to a fixed regression repeatability model under different phenotyping strategies. Data were available from repeated measurements of CH 4 concentrations (ppm), recorded in the feed bins of milking robots, on 52 commercial dairy farms in the Netherlands. In total, 36,370 averaged weekly records were available from 4,664 cows. Genetic parameters were estimated using a fixed regression model, and a RR model with 1st to 5th order Legendre polynomials for the additive genetic and within lactation permanent environmental effect. The mean heritability was 0.17 ± 0.04, and the mean within lactation repeatability was 0.56 ± 0.03. The genetic correlations between days in milk were high and ranged from 0.34 ± 0.36 to 1.00 ± < 0.01. Permanent environmental correlations showed large deviations and ranged from −0.73 ± 0.08 to 1.00 ± < 0.01. With a large number of full lactation daughter CH 4 records per bull, the reliability was not sensitive to using the fixed versus RR model. However, when shorter periods were recorded at the start and end of the lactation, the fixed regression model resulted in a loss of reliability up to 28% for bulls. Assuming the fixed model when the true (co)variance structure is reflected by the RR model, more than twice as long recording from the start of lactation was required to achieve maximum reliability for a bull. Thus, a too simplistic model could result in implementing too little recording, and lower genetic gains than predicted from the reliability.


INTRODUCTION
Enteric fermentation by ruminants is the source of approximately 6% of all anthropogenic greenhouse gas emissions, and ruminant emissions are expected to increase because of the increasing global demand for meat and milk and therefore more ruminants will be kept worldwide (Beauchemin et al., 2020).Ruminant greenhouse gas emissions, with methane (CH 4 ) as the main contributor, have a significant effect on climate change, as CH 4 has a global warming potential approximately 27 times greater than CO 2 over a 100 year lifespan (IPCC, 2021).There is a worldwide interest in applying selective breeding to lower the environmental impact of enteric CH 4 emissions from cattle, because it offers a cost-effective solution and the effect is permanent and cumulative (Lassen and Difford, 2020).
To apply selective breeding, models that accurately estimate breeding values for individual animals and correlations with other important breeding goal traits will be needed.The goal of selective breeding is to reduce CH 4 emissions along full lactations.However, recording CH 4 is challenging and requires high levels of labor, knowledge, and expensive specialized equipment, and measurements of the total amount of CH 4 that is emitted during full lactations are not easily available for large number of cows (Garnsworthy et al., 2019).Equipment is often routinely exchanged between farms for short periods of time to maximize the number of cows that can be recorded.As a result, it is important to determine the optimal moment and length of recording, to ensure highest accuracy with CH 4 emissions for the full lactation.Also, the best model to be used for genetic evaluations is critically dependant on the Genetic parameter estimates for methane emission during lactation from breath and potential inaccuracies in reliabilities assuming a repeatability versus random regression model genetic covariance structure of methane records collected during part of the lactation.Most initial studies that estimated genetic parameters for CH 4 emissions of dairy cattle applied repeatability models (Lassen and Lovendahl, 2016;van Engelen et al., 2018), which generally assume that the genetic correlation between records taken during different stages of a lactation are equal to unity.However, previous studies that used random regression (RR) models, which allow for different genetic variances and covariances between time points, have shown that the repeatability and heritability of CH 4 concentrations measured from breath change over a lactation (Breider et al., 2019;Manzanilla-Pech et al., 2022;Pszczola et al., 2017;Sypniewski et al., 2021).Differences in variances and correlations over the lactation may have important implications for the genetic gains that can be predicted for, and expected from, future breeding programs that aim to reduce CH 4 emissions.For example, the genetic correlation between lactation stages is a good indicator for if CH 4 concentrations should be modeled with random regressions, as multiple traits, or if using a single trait in a repeatability model is sufficient.
For the repeatability of daily CH 4 concentrations, Pszczola et al. (2017) reported a steep decrease in repeatability in mid lactation, ranging from 0.40 to 0.18, whereas Manzanilla-Pech et al. ( 2022) reported the highest repeatabilities in mid lactation for weekly CH 4 concentrations, ranging from 0.63 to 0.86 (median SE = 0.03).For the heritability of CH 4 concentrations, Manzanilla-Pech et al. (2022) and Pszczola et al. (2017) reported relatively stable heritabilities over a lactation which reached the maximum in mid lactation of 0.10 to 0.28 (median SE = 0.05) and 0.27 ± 0.12 to 0.30 ± 0.08, respectively.In contradiction, Breider et al. (2019) reported that heritabilities for weekly CH 4 concentrations were lowest in mid lactation and reached the maximum in late lactation (h 2 = 0.12 ± 0.16 to 0.45 ± 0.11), and Sypniewski et al. (2021) reported a steady increase of heritability of daily CH 4 concentrations over the lactation (h 2 = 0 to 0.14), with the heritability being zero at the start of the lactation.Thus, a consensus has not been reached to date, on the shape of the curve of the repeatability or the heritability estimated over a lactation.Differences in results may have arisen from including records on a relatively small number of cows, ranging from 184 to 575, coming from at most 2 farms.As a result, when reported, standard errors were high.Furthermore, for 3 of the 4 studies the data recording period did not cover a full lactation.As explained above, being able to model the covariance structure over the full lactation correctly helps to draw inferences from longitudinal data.Therefore, our objective was to study the heritability and repeatability of measured CH 4 concentrations over a lactation using a RR model, and to compare the accuracy of using a fixed regression repeatability model or the RR model under different phenotyping strategies.For the analyses we used a novel data set from emissions measured on the largest number of cows to date from across the Netherlands, which has previously only been analyzed using repeatability models (van Breukelen et al., 2022).

Methane recording
Data were collected on 7,097 Holstein cows, from 54 farms (Table 1), using non-invasive sniffers (WD-WUR v1.0 and v2.0, manufactured by Carltech BV), that sampled CH 4 and carbon dioxide (CO 2 ) concentrations (ppm) from the feed bin of automated milking systems (AMS).Ethical approval was not needed for this study because no animal procedures were performed.The number of cows recorded per farm ranged from 43 to 299.Various types of AMS systems were present in the study, manufactured by: DeLaval (DeLaval BV), Fullwood (Fullwood Packo BV), GEA (GEA Group), Lely (Lely Industries NV), and SAC (SAC BV).On each farm, at most one AMS unit was equipped with a sniffer.This strategy was decided on, to maximize the numbers of farms with phenotyping, and thereby maximizing the number of cows that were phenotyped.The sniffers were calibrated to measure CH 4 concentrations (CH 4 c) in a range of 0 to 2,000 ppm.Measured concentrations can be used as a proxy for CH 4 emissions (g/day), similar to what has been reported in previous studies (Difford et al., 2020;van Breukelen et al., 2023).Data were recorded between March 2019and March 2023. Between March 2019and February 2021, v1.0 sniffers were used of which the data collection process is described in more detail in van Breukelen et al. (2022).Between December 2021 and March 2023, v2.0 sniffers were used.The 2 versions of sniffers functioned similarly, however, the v2.0 sniffers measured concentrations every 5 s, opposed to the longer recording intervals of the v1.0 sniffers (ranging from 10 to 35 s).Furthermore, the v2.0 sniffers had improvements to the housing and data sharing, which benefitted the ongoing data collection in the barn environment but did not change how CH 4 was measured.

Data editing
The sniffer records were filtered to exclude data which were biologically improbable, for example due to blocking of the sampling tube in the AMS feed bin by dust from pellets.To do so, for each hour within a farm, data recorded within that hour would be discarded if: 1) the mean was below 30 ppm CH 4 , 2) the inter quartile range was below 200 ppm CH 4 , 3) the maximum was above 3,500 ppm CH 4 , or 4) if at least 30% of the data would fall within the range of the first and second mode, plus and minus 10 ppm CH 4 .Furthermore, individual outliers were discarded, which were defined as being lower than −200 ppm CH 4 , or values outside of the upper and lower 0.001th quantile.A threshold below zero was used, because sensors could measure below zero if the calibration drifted.Nonetheless, calibration lines would remain linear and could therefore still provide important information about variation in emissions on that farm.After filtering, background concentrations were estimated from the 0.001 lowest quantile of each day, and subtracted from each measurement.
The timestamps (on the level date-time, with time as hours, minutes, and seconds) from sniffer data and AMS data (also provided by CRV) were used to connect cow ID to each individual sniffer record.The AMS and v2.0 sniffers automatically synchronised to real time through the Odido (Odido Netherlands BV) mobile network, to prevent drift of the clocks; for the data recorded by v1.0 sniffers the alignment process is described in van Breukelen et al. (2022).In addition, the alignments were confirmed visually for each farm, by observing if CH 4 concentration would increase after the start of a milking and remain low when no cow was in the AMS.After the alignment, the mean CH 4 concentration per AMS visit was calculated for each individual visit.In calculating the mean, only measurements recorded from the first and up to the fifth minute of milking were used.The other measurements were discarded to account for the delay in the air sample reaching the sniffer, and to exclude records when cows are likely to have finished eating the pellets provided in the AMS feed bin, which makes it more likely that cows move the head away from the sampling inlet.Thereafter, means that were derived from less than 2 and a half minutes of milking were discarded, to ensure that multiple records were collected during the milking visit, and that these records included belching events (van Soest, 1994).The aligned data were aver-aged per AMS visit, and contained 661,917 records of mean CH 4 c (Table 1).
A linear model was used to correct each record for diurnal variation in measured CH 4 concentrations within farm, with a Fourier series approach (Lassen and Lovendahl, 2016;Lovendahl and Bjerring, 2006), using the following model: where y ij is the phenotype for the mean CH 4 c per AMS visit; Farm is the fixed effect for the ith farm; θ is a decimal fraction of the time of measurement, following a 24-h diurnal cycle (i.e., θ = hour at measurement / 24), where j is the order of regression; and e is the residual error, following e i ~ N 0 2 , , Iσ e ( ) where I represents the identity matrix and σ e 2 the residual variance.To derive the corrected estimates for each visit, the estimated fixed effects were subtracted from the corresponding measurement.
The mean CH 4 c per AMS visit were then summarized as weekly mean concentrations (Table 1).If the weekly mean of a cow consisted of less than 7 AMS visits, then the weekly record was discarded.Two farms had a sniffer installed for only a short period of time and had issues with data collection, and because of the limited number of records per cow, weekly means were not available and the 2 farms were discarded (n = 52).On average, a weekly mean record consisted of 11 recorded visits (min-max: 7-33).For each weekly measurement, the associated days in milk were reported as the first day that was included in the weekly measurement.Large differences were observed in the mean and standard deviation of records taken with the different versions of sniffers.Therefore, CH 4 c was scaled to a standard deviation of one and centered within version of sniffer, by subtracting the mean of all records, measured with the corresponding sniffer version, from each measurement and dividing the result by the standard deviation of all measurements with the corresponding sniffer version (e.g., a CH 4 c record taken by a v1.0 sniffer was scaled with the overall v1.0 sniffer mean and standard deviation).The data were not log-transformed, as this had little impact on normality of the data in previous analyses (van Breukelen et al., 2022), and the residuals from the analyses appeared normally distributed.For the genetic analyses, the data set was filtered to include only records up to 405 d in milk (DIM) and the animals that were at least 75% Holstein (Table 1).The number of recorded and cows after each filtering step, and the final numbers used for the analyses are reported in Ta-  1.
The reduced number of records in late lactation appear to be related to fewer visits to the milking robot, which aligns with the natural decline in milk yield for cows later in their lactation, possibly in combination with other factors.On average, the final data set included 8 weekly mean CH 4 c records per cow, with minimum of 1 record and a maximum of 37 records.Per farm, on average 19 weeks were recorded, with a minimum of 1 week and a maximum of 63 weeks.

Pedigree data
Pedigree and other cow information were provided by CRV (Arnhem, the Netherlands).The pedigree was pruned to include all cows with a CH 4 record and their ancestors using the R-package "pedigreeTools" in R v4.2.0.The pruned pedigree was 25 generations deep and contained 48,926 animals.

Genetic parameter estimation
To investigate if measurements from the 2 sniffer versions are interchangeable, fixed regression bivariate analyses were performed that applied CH 4 c from the v1.0 or v2.0 sniffers as different traits, to estimate genetic correlations, using the model defined in Equation (2) as a bivariate model.The bivariate model excluded the across parity permanent environmental effect due to convergence issues.Nonetheless, a univariate model with the across parity permanent environmental random effect was also used to estimate the heritability for each sniffer version.Since this effect was included in the further analysis.
Thereafter, a single-trait RR model was used to estimate variance components, using a restricted maximum likelihood method in ASReml 4.2ng (Gilmour et al., 2015).Different orders of Legendre polynomials for the random genetic and permanent environmental effect were fitted, ranging from the 0th to the 5th order, and compared using the Loglikelihood, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC).The model using Legendre polynomials of the 0th order is equal to using a fixed regression repeatability model.The following model was used to estimate genetic parameters: where y ijlk is the scaled and centered phenotype for the mean CH 4 c per week; µ is the mean; HYW is the fixed effect for the interaction of herd, year, and week of measurement i; Par is the fixed effect for parity (j = 1 to 4, where 4 includes parity 4 or higher), and is fitted as an interaction with β, which is a fixed regression coefficient with third-order Legendre polynomials, measured at t DIM, for the kth regression coefficient of animal l; Breed is the fixed effect for the second breed, and is fitted as an interaction with β, which is a fixed regression coefficient with second-order Legendre polynomials measured as a breed fraction u that is other than Holstein derived from the pedigree; a lk and pepar lk are RR coefficients for the genetic effect and the permanent environmental effect within parity, respectively, and ø is the term for the nth order Legendre polynomial (ranging from zero to 5) at t DIM; pe is the random permanent environmental effect across lactations; and e ijl is the random residual.The residual error was assumed to have heterogeneous variances and was divided into 5 classes for DIM (0-59, 60-119, 120-239, 240-359, and 360-405).Although Legendre polynomials give a rank reduction compared with a full rank matrix between all DIM, higher order Legendre polynomials might still give convergence issues in ASReml.This is due to the average information algorithm struggling with close to non-positive definite matrices, and the expectationmaximization algorithm being notoriously slow.Therefore, factor analytical modeling, using XFA (extended factor analysis) inflation factors, was used in ASReml, to decompose and reduce the rank of the variance-covariance matrix of the Legendre polynomials (Thompson et al., 2003).The results reported in the paper applied several XFA equal to the order of Legendre polynomials, except when fitting polynomials of the fourth and fifth order.Because of convergence issues, here the third and fourth XFA were fitted for the models with fourth, and fifth order Legendre polynomials, respectively.Estimated (co)variance components were used together with the Legendre polynomial coefficients at each DIM to estimate the full (co)variances matrices and the genetic parameters (heritability, repeatability, and correlations between DIM) and their approximate standard errors (SE), as described by Fischer et al. (2004).In short, the genetic, permanent environmental, and phenotypic (co)variances were estimated by: G = ΦKΦ′, PE = ΦKPEΦ′, and where G is the genetic (co)variance matrix; PE is the permanent environmental (co)variance matrix, to which the across parity permanent environmental variance was added to all (co)variance elements in the matrix; and P is the phenotypic (co)variance matrix per DIM (n*n, where n is the level of DIM, consisting of 58 classes of 7 DIM, up to 400 DIM); Φ is a matrix of order t*n, where t is equal to the number of orthogonal polynomial coefficients; K and KPE are matrices of order t*t, which contain the estimated covariance functions that describe the genetic (co)variance components and permanent environmental (co)variance components, respectively, for each RR coefficient; and σ 2 e is the residual variance that is added to the diagonal of the phenotypic (co)variance matrix, with residual (co)variances assumed to be zero.
The heritability (h 2 ) was defined as: the within lactation repeatability (r) was defined as: where σ a 2 is the additive genetic variance; σ wpe 2 is the permanent environmental variance within lactations; σ ape 2 is the permanent environmental variance across lactations; and σ e 2 is the residual variance.

Selection index calculations
To predict the reliability using genetic evaluation models assuming different orders of Legendre polynomials for the random effects, the associated selection index coefficients were estimated, following: where b is the selection index coefficient, P −1 is the inverse of the phenotypic (co)variance matrix among observations in the selection index (n*n, where n consists of 58 potential classes of weeks in milk where recording takes place), G is the n*m matrix of genetic (co)variance matrix among the n weeks in P and all m traits in the aggregate genotype (where m consists of 58 classes of weeks in milk), and v is a vector with weights, which were kept equal to for al 58 weeks in milk.Because for all orders the smallest eigenvalues of the G matrix were slightly negative, the matrix was first made positive semidefinite.This was done by changing negative eigenvalues to 0.001.Thereafter, the matrix was again constructed by multiplying the eigenvectors with the transformed diagonal of the matrix and the transpose of the eigenvectors.Nonetheless, the non-positive definite matrices yielded similar results (not shown).
Thereafter, the accuracy of the indexes (r HI ) were estimated using the following formula: Which applied the parameters from the previous formula, and the matrix C which is a m*m matrix with genetic (co)variances in the aggregate genotype.From this, reliabilities were predicted for the different models as the squared accuracies (r 2 ).These reliabilities reflect the model reliabilities that would be published for breeding values given that the assumed model was used (regardless of the true genetic parameters).
The reliability was also calculated for an additional scenario, i.e., where the repeatability model was used, but in fact the estimated parameters from the RR reflect the true parameters more closely.In that scenario the b-values reflect a "sub-optimal" index using the fixed regression repeatability model, whereas the RR model gives the "optimal" index.
Optimal indices and models also depend on the recording strategy, and the number of records, since P will be affected.Therefore several scenarios were evaluated.For the first set of scenarios, it is assumed that for each recorded week in the lactation one record is available per cow (scenarios "Cows").The P matrices simply reflect the phenotypic (co)variances in this case.In the second set of scenarios, the P matrices were replaced with the (co)variances of the G matrix, simulating sires with a large number of daughters recorded, (scenario "Bulls").The large number of records provide close to the true genotype for a bull at the recorded moments.
Taking the 2 models, and the Cow and Bull situation, scenarios were investigated where observations were only available for parts of the lactation.This is relevant for example when cows are phenotyped for CH 4 emissions with expensive equipment, which is routinely exchanged between farms for short periods of time.The scenarios were: 1) with records on different parts of the lactation, where always 8 weeks (56 d) were recorded sequentially; 2) with records on different parts of the lactation, where always 3 weeks (21 d) were recorded sequentially; and (3) starting with one observation at the start of the lactation, which cumulatively increased by one week of recording until 400 DIM was reached, again changing the P and G matrices based on the DIM with records.

Exploratory analyses
Before standardisation, the mean CH 4 c of v1.0 sniffer sensors was 512 ppm, with a standard deviation of 172 ppm, and from v2.0 sniffers the mean was 719 ppm, with a standard deviation of 313 ppm (Table 2).The coefficients of variation were 34% and 44%, for mean CH 4 c measured by the v1.0 and v2.0 sniffer sensors, respectively.The overall mean CH 4 c per week was 668 ppm, with a standard deviation of 299 ppm.For the overall mean, the coefficient of variation was 45%.After standardisation and filtering of the data, the mean scaled CH4c of v1.0 sniffer sensors was 3.74, with a standard deviation of 1.01, and from v2.0 sniffers the mean scaled CH4c was 2.03, with a standard deviation of 1.00.The overall mean scaled CH4c per week was 3.7, with a standard deviation of 1.24.The total number of recorded DIM for an individual cow ranged from 1 to 271, and the total number records per DIM are plotted in Figure 1.
From the data including only measurements from v1.0 sniffers, the heritability and repeatability were 0.27 ± 0.03 and 0.62 ± 0.01, and for v2.0 they were 0.37 ± 0.03 and 0.76 ± 0.01.The genetic corelation between the 2 sensors was 0.99 ± 0.09.Nonetheless, there was a considerable difference in the estimated variance components.For v1.0 and v2.0 data, respectively, the phenotypic variances were 15,759 ± 468 and 52,662 ± 1,208, and the genetic variances were 4,180 ± 620 and 19,368 ± 1,936.Thus, in general the measurements taken with v2.0 sniffers were associated with higher variances compared with v1.0 sniffers, justifying the need to standardize the sniffer data to a common phenotypic variance.But because of the high genetic correlation and not majorly different heritabilities, it was decided to further analyze the data using univariate models, where we assume that after phenotypic standardization the records taken by either version of sniffer is genetically the same trait.

Order of the random effects
The random effects of the RR model were modeled with different orders or Legendre polynomials (0 to 5th).The fit of the models was investigated based on the Log-likelihood, AIC and BIC.The models that applied higher orders of Legendre polynomials had consecutively better goodness of fit, indicated by the lower AIC and BIC, and the higher Log-likelihood (Table 3).

Genetic parameter estimates
The estimated heritability, using a fixed regression repeatability model was 0.34 ± 0.02, and the repeatability was 0.73 ± 0.01.However, when including an across lactation permanent environmental effect (as is applied in the further analyses), the heritability and repeatability were lower.Here, the heritability was 0.18 ± 0.03, and ranged between 0.17 and 0.18 for the different groups of residual variances.The average within lactation repeatability was 0.48 ± 0.03, and ranged between 0.46 and 0.49 for different groups of residual variances.The average across lactation repeatability was 0.41 ± 0.02, and ranged between 0.39 and 0.42 for different groups of residual variances.
The lactation pattern for the heritability of the random effects modeled with different orders of Legendre polynomials was similar between orders (Figure 2, top).The heritability of the fixed regression repeatability model (Leg 0) was higher in the beginning and end of the lactation, whereas the models with first or higher order Legendre polynomials had a higher heritability throughout mid lactation.For the fourth order RR model, the heritability of weekly mean CH 4 c was on average 0.17 ± 0.04 and increased steeply in the first 130 DIM.The maximum heritability was 0.21 ± 0.03 at 134 DIM.Thereafter, the heritability decreased, first slowly and then steeply, to a minimum of 0.10 ± 0.04 at 358 DIM.The heritability of the model with Legendre polynomials of the fifth order deviated from the other estimates at the end of the lactation, where they peaked at a level of 0.33 at 393 DIM.
Similarly, the lactation pattern for the within lactation repeatability of the random effects modeled with different orders of Legendre polynomials was similar and were within a range of 0.45 to 0.83 (Figure 2, bottom).The repeatability from the fixed regression repeatability model (Leg 0) was lowest throughout the full lactation.Throughout mid lactation, the repeatability of the second, third, fourth and fifth order models were close to identical, whereas the first order model had a lower repeatability between 75 and 210 DIM, and a higher repeatability after 250 to 380 DIM.For the fourth order RR model, the repeatability remained stable during mid lactation, at approximately 0.56 ± 0.03.The minimum repeatability was 0.49 ± 0.04 at 351 DIM, and the end of the lactation the repeatability peaked to 0.70 ± 0.04 at 400 DIM.The across lactation repeatability of the fourth order model ranged from 0.28 ± 0.09 to 0.45 ± 0.03, with an average repeatability of 0.41 ± 0.03.
The underlying variance estimates, using the second to the fifth order models, are similar (Figure 3).The additive genetic variance estimates of the third to fifth order models spiked at the end of the lactation after approximately 320 to 375 DIM.This spike was the largest for the model with fifth order regressions.The within parity permanent environmental variance was similar for the first to the fifth order, and also spiked at the end of the lactation.For the model of the fifth order, the permanent environmental variance did not increase greatly at the end of the lactation and had a maximum of 0.22.
The residual variance, which was modeled in 5 classes of DIM, consistently decreased for higher orders of Legendre polynomials that were fitted (Figure 4).Especially the fixed regression repeatability model (Leg 0) and the model using a first order RR had somewhat higher residual variances.The estimates for second to fifth order RR were similar, especially during mid lactation.

Genetic correlations between DIM
The genetic correlations for the mean CH 4 c per week between DIM were high and many approximated one, except for measurements at the start and end of the lactation (Figure 5, left).The genetic correlation was on average 0.91 ± 0.08.The lowest correlation was 0.34 between the 302nd and 400th DIM, and had a large SE of 0.36.
The within lactation permanent environmental correlations between mean CH 4 c per week between DIM were high between measurements taken close in time  (up to 1.00 ± < 0.01).Whereas the permanent environmental correlations were low, and in some cases negative, for measurements taken further apart in time (up to −0.73 ± 0.08, between 1 and 393 DIM) (Figure 5, right).The mean permanent environmental correlation was 0.39 ± 0.06.

Selection index calculations
Selection index coefficients and the model reliability of the index (i.e., published breeding values) were predicted to evaluate the effect of using a fixed or a RR model in the genetic evaluations.First assuming that recording happened during the whole lactation period.The reliability of the selection index (r 2 ) in the scenarios for cows (using the phenotypic (co)variance matrix as P) was highest at 0.43 for the model using fifth order RR (Table 4).This was followed by the model with second, third, and fourth order regressions at 0.29.The model with the lowest reliability was the fixed regression repeatability model (Leg 0) at 0.24.Thus the published reliability will be lower with the fixed regression model.However, the reliability, when using the b values for the fixed regression model (sub-optimal), Leg 0, when assuming the (co)variance estimates from the fifth order RR model was still 0.31.Thus when using the fixed regression repeatability model in the genetic evaluation model, only 71% of the reliability would be realized when we assume that the fifth order RR model yields true (co)variances.However, for the first to fourth order RR model, the difference with the .Heritability (top) and within lactation repeatability (bottom) over a lactation (DIM; days in milk) for the mean CH 4 concentration per week, with standard errors, using a random regression model with random effects fitted with Legendre polynomials of the 0th (Leg 0; blue), 1st (Leg 1; green), 2nd (Leg 2; gray), 3rd (Leg 3; orange), 4th (Leg 4; black), and 5th (Leg 5; yellow) order fixed regression model were smaller, and ranged from 92% to 98%.
The reliabilities when we assumed a large number of daughter records are available for each individual bull (the bulls scenario) were all one.Similarly, when using the b values of the fixed regression repeatability model with the genetic (co)variance matrix instead of the phenotypic (co)variance matrix (sub-optimal) the repeatabilities in all scenarios were also one, thus there was little loss using the simple repeatability model.
The results reported above assume that full lactations are recorded for individual cows.In future phenotyping strategies, it is possible that data on weekly CH 4 c, or other CH 4 traits, are not available throughout a full lactation, but still the interest is in a breeding value predicting the full lactation.Therefore, we predicted the reliability of the breeding values for the full lactation, with limited data available during various parts of the lactation (using a fourth order RR model, and a fixed regression repeatability model).Overall, reliabilities for the breeding values from measurements taken during mid lactation (50 to 302 DIM) were predicted to be higher than the reliabilities of measurements taken at the start and end of the lactation (Table 5).When using the fixed regression repeatability model, and assuming the (co)variances estimated from the RR model are reality (sub-optimal), when measuring during the first 57 DIM only 93% of the reliability would be realized.
Using shorter recording periods, of 3 weeks, also showed the higher value of mid lactation records (Table 6, see Appendix Table A1 for the full table).However,   4. The reliability of the selection index (r 2 ) for models using different orders of Legendre polynomials (Leg) for random effects (optimal (opt)), when using the phenotypic (co)variance matrix (Cows), or genetic (co)variance matrix as phenotypes (Bulls), and the two scenario's while using the selection index coefficients coming from the fixed regression repeatability model (sub optimal (sub opt)), with the difference between the two as a percentage the estimated reliabilities for bulls decreased notably at the start and end of the lactation.When using the fixed regression repeatability model (Leg 0), while assuming the estimated (co)variances from the RR are reality, only 73% of the reliability was realized with measurements taken in the first 3 weeks of the lactations.Nonetheless, between 29 and 372 DIM the differences were small, and high reliabilities for bulls could be achieved with either model (0.97 -0.98).
When looking at the cumulative trend of collecting weekly measurements over a lactation on the reliability, large differences can be observed in the increase in reliability at the start of the lactation (Figure 6).When using the Leg 0 model, the reliabilities will be inflated at the start of the lactation, both in the Cows and especially the Bulls scenarios (Figure 6).Predicting a reliability of one after one week of recording many daughters, whereas the true reliability is much lower.Also, it takes more weeks to achieve a true reliability of selection, when using the suboptimal fixed regression model versus the RR model (approximately 90 versus 30 weeks).

DISCUSSION
The aim of this research was to study the heritability and repeatability of measured CH 4 concentrations over a lactation using a RR model, and to compare using the fixed regression repeatability model and a RR model for genetic evaluations under different phenotyping strategies.The data used was a novel data set from emissions measured on 4,664 cows from 52 farms located across the Netherlands, with up to 37 weeks of the lactation recorded, which has previously only been analyzed using fixed regression repeatability models.

Order of polynomials for RR models
In this study we estimated genetic parameters, using different orders of Legendre polynomials in a RR model, ranging from the 1st to the 5th order.In our study, the models with higher order polynomials had better goodness of fit (Table 3).This is similar to what has been shown in studies on random regressions for milk yield using test-day models (Li et al., 2020;Pool van Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION Table 5.The reliability of the selection index (r 2 ) with phenotypes only being available within sequential ranges of 56 d in milk (DIM) (eight weeks), when using the phenotypic (co)variance matrix (Cows), and genetic (co)variance matrix as phenotypes (Bulls) for the random regression model of the fourth order (optimal (opt)), with using the selection index coefficients coming from the fixed regression repeatability model (sub optimal (sub opt)) and the (co)variance matrices from the RR model, and the difference between the two as a percentage  6.The reliability of the selection index (r 2 ) with phenotypes only being available the first and last three weeks of lactation as sequential ranges of 21 d in milk (DIM), when using the phenotypic (co)variance matrix (Cows), and genetic (co)variance matrix as phenotypes (Bulls) for the random regression model of the fourth order (optimal (opt)), with using the selection index coefficients coming from the fixed regression repeatability model (sub optimal (sub opt)) and the (co)variance matrices from the RR model, and the difference between the two as a percentage (For the full  et al., 2000;Van Der Werf et al., 1998).However, as these studies also highlight, fitting higher order polynomials comes at a higher computational cost.Since for the higher orders fitted the increase in goodness of fit becomes less, the common consensus between the previous studies is that fitting at least third order polynomials are considered sufficient when using heterogeneous residual variances.Furthermore, fitting higher order polynomials can induce oscillatory patterns along the lactation, which are unlikely to be biological (Pool et al., 2000).This highlights that more parsimonious models should be preferred when the improvement in model fit is limited.In the final results of this paper, the model fitted with fourth order Legendre polynomials was applied for the random additive genetic and within lactation permanent environmental effects, as the fifth order model resulted in large deviations in the heritability at the end of the lactation (Figure 2) and in large deviations in genetic correlations in the beginning and end of the lactation, which are difficult to explain (result not shown).

Heritability and repeatability within lactations
The heritability for weekly mean CH 4 c increased steeply in early lactation, peaked at 134 DIM, whereafter it steadily decreased again (Figure 2, top).The average heritability was 0.17, with a low standard error of 0.04.The standard errors of the estimates are lower than of previous studies, giving more confidence in the estimates of this study (Breider et al., 2019;Manzanilla-Pech et al., 2022;Pszczola et al., 2017).The pattern of the heritability within a lactation showed similarities to the study by Pszczola et al. (2017), who reported that the heritability for CH 4 p (CH 4 production, estimated from sniffer CH 4 c measurements) was highest in mid lactation.Although, the reported heritability estimates were less variable and ranged between 0.23 ± 0.12 and 0.3 ± 0.08.A similar pattern was also observed in the study by Manzanilla-Pech et al. (2022), where the heritabilities for CH 4 c ranged between 0.10 and 0.28.In a study by Breider et al. (2019), a different pattern was observed for the heritability of CH 4 c over a lactation, compared with the afore mentioned studies.Breider et al. (2019) showed that the heritability was the lowest in mid lactation, and increased toward the end of the lactation.The study by Manzanilla-Pech et al. (2022) discussed that this was possibly an outcome of the inclusion of multiple lactations in the analyses by Breider et al. (2019), whereas the former study analyzed first and second lactation separately.However, in this study multiple lactations were also combined.Another factor that could have influenced the heritability estimates was that Breider et al. (2019) applied 2 bivariate models, with CH 4 c modeled with milk yield and body weight, respectively.Also, the study by Sypniewski et al. (2021) reported a different pattern of heritability over the lactation.There, the heritability for CH 4 c increased throughout the full lactation, whereas the pattern for the heritability of CH 4 p was similar to this study.In the study by Sypniewski et al. (2021) the residual was modeled with homogeneous variances, which may have resulted in a different partitioning of variance between the genetic and permanent environmental variance within the lactation compared with this study.
For all order models fitted, the within lactation repeatability was relatively stable throughout the lactation (at approximately 0.58, Figure 2, bottom).This indicates that the within animal environmental variance is similar for all measures taken throughout the lactation, and that there are no parts of the lactation that would benefit form a larger number of measurements than other parts.The within lactation repeatability reported in this study was of similar pattern, although slightly lower, than what was reported in the study by Manzanilla-Pech et al. (2022) for CH 4 c, where they ranged between 0.63 and 0.86.However, in the study by Manzanilla-Pech et al. (2022) a separate across lactation permanent environment effect was not fitted.Including the across lactation permanent environmental effect, as the total repeatability, the repeatability would be 0.80 and thereby similar to what was reported previously.Other estimates for the repeatability of CH 4 c within lactations have not been previously reported.
The increase in heritability at the start of the lactation was a result of lower additive genetic variance at the beginning of the lactation, which can be observed in Figure 3.However, caution should be taken in interpreting the first and last days of the lactation, as models fitted with polynomials are known to deviate at the far extremes (Pool et al., 2000).The additive genetic and the permanent environmental variance started to deviate before approximately 50 DIM and after 300 DIM (for the models fitted with second to fifth order Legendre polynomials).A similar pattern has been reported in studies on milk yield using RR models (Pool et al., 2000;Strabel et al., 2005;Van Der Werf et al., 1998).Higher permanent environmental variances at the beginning and end of the lactation may be attributed to a lower number of records that often are available at the start and end of the lactation.It is therefore expected that the large increase in permanent environmental variance at the end of the lactation in this study was in part an artifact of the data.From Figure 1 (left) it can be observed that the number of records at the beginning and at the end of the lactation were low, with under 50 weekly mean CH 4 c records per DIM after 300 DIM, and decreasing.The deviations influenced the heritability before 10 DIM and after 300 DIM, which is underlined by the higher associated standard errors at the beginning and toward the end of the lactation.The extent of the effect of a lower number of records at the beginning and end of the lactation should be further investigated, and could be overcome, for example, by fitting models using a smoothing factor such as splines (White et al., 1999).

Genetic and permanent environmental correlations within lactations
The genetic correlations between DIM were high, except between some DIM at the far ends of the lactation (Figure 5, left).High and positive genetic correlations indicate that, even for most records taken far apart in time, the direction of selecting for lower CH 4 c would be the same.In addition, the magnitude of the genetic correlations have an influence on the reliability of the selection index of measurements taken along the lactation, which will be discussed later.The genetic correlations of weekly mean CH 4 c between DIM were similar to what has been reported in previous studies (Manzanilla-Pech et al., 2022;Pszczola et al., 2017).Nonetheless, there were some differences.The study by Manzanilla-Pech et al. (2022) reported that the genetic correlations were higher and close to one along the full lactation in the second lactation, in contrast to the first lactation where the genetic correlations reduced for weeks that were the furthest apart in time.This should be investigated further, as it may implicate that a different phenotyping strategy should be maintained for first or later lactation cows.Nonetheless, the model in this study included a between parity random effect and a fixed effect for the lactation curve as an interaction with parity, to correct for differences between parities.
The permanent environmental correlations between weekly mean CH 4 c measurements closely together in time were high, but the correlations reduced and became negative for measurements further apart in time (Figure 5,right).This confirms that modeling CH 4 including a permanent environmental effect is the most appropriate, ideally with random regressions, as different lactation stages are associated with differences in the permanent environmental variance.Permanent environmental correlations have not been previously reported in the literature.
Similar to the heritability and repeatability, the genetic and permanent environmental correlations may have been influenced by the lower number of records at the start and toward the end of the lactation.Therefore, caution should be taken in interpreting the results before 10 DIM and after 300 DIM, and the results should be confirmed using a data set that includes a

Selection index calculations
The selection index coefficient is directly linked to the maximum expected response in the aggregate genotype, and the reliability that would be published next to the breeding values.Multiple scenarios were simulated, and each scenario was simulated in 2-fold: 1) using the phenotypic (co)variance matrices (Cows), and 2) replacing the phenotypic relationship matrix by the genetic (co)variance matrices (Bulls).The first represented scenarios where breeding values would be estimated for a cow with single measurements at each class of DIM, whereas the second represented scenarios where breeding values would be estimated for a bull with a large or infinite number of daughters with many records.In the second case it is expected that the phenotypic (co)variance matrix approximates the genetic (co)variance matrix, and thus breeding values would become 100% reliable.This was indeed the case when full lactations would be measured, as can be observed in Table 4.
The reliabilities for a cow with measurements of a full lactation were higher for the RR model with high orders.Especially, the fifth order RR model resulted in a large increase in the estimated reliability and was 0.43.Nonetheless, this high reliability should be interpreted with care, and should be confirmed with using a model that applies a smoothing factor such as splines, as the results for the fifth order RR model showed deviations from the other models.Therefore, in this study we focused on the RR model using fourth order Legendre polynomials.
When using the sub-optimal selection index coefficients of the fixed regression repeatability model, with (co)variance estimates of a fourth order RR repeatability model, an almost identical reliability was reached as when using the optimal selection index coefficients of RR model (98% for cows and 100% for bulls, Table 4).That would suggest that there is not a large gain in using a RR model over a fixed regression repeatability model.Nonetheless, for the bulls scenario the sub-optimal fixed regression repeatability selection index coefficients always realized 100% of the reliability compared with the RR model, suggesting that for bulls with a large number of daughters recorded over the full lactation there is never an advantage in using a RR model over a fixed regression repeatability model.Although this might not be a practical scenario, it reflects the loss in selection reliability, based on the genetic architecture of the trait.That is obviously small when the true breeding value is known at any moment in time, because the genetic correlation estimates between the time periods play no role in the genetic evaluation.
When measurements are available over only a short period of time, the RR repeatability model will most likely perform better compared with a fixed regression repeatability model.To investigate the effect of shorter recording period, 3 simulations were performed, with: 1) sequential periods of 8 weeks of recording, 2) sequential periods of 3 weeks of recording, and 3) starting at one DIM and cumulatively adding one week of recording up to 274 DIM.This is relevant for phenotyping strategies, with only limited logistic and/ or financial resources available, by which they are not able to record full lactations, and is especially relevant in farming systems with seasonal calving.When recording short periods of the lactation, the realized reliability was lower than when a full lactation would be recorded, with increasingly lower reliabilities for shorter recording periods.The difference in reliabilities was smaller for records taken during mid lactation (Table 5 and 6).For short recording periods, from the data reported in this study it can be concluded that records taken between approximately 100 DIM and 200 DIM will most likely yield the highest reliability and thereby the highest genetic progress, regardless of the model chosen.Thus, a fixed regression repeatability model in this case can perform sufficient, and might be the better choice when there is limited data available or limited computational power.
For short recording periods at the start of the lactation (1 to 100 DIM), the reliabilities of the fixed regression repeatability model were higher than the reliabilities of the RR model for the cows and bulls scenarios (Figure 6).Here, the fixed regression model likely overestimates the reliability, as this model assumes equal genetic correlations between DIM, whereas the RR model showed that the first DIM had lower genetic correlations with the rest of the lactation.This has important implications for practice, and the results suggest that when breeding values are estimated from a fixed regression model, and are only recorded during the first ± 50 DIM, the realized genetic gain may become lower than expected.Therefore, in these situations using a RR model should be preferred.

CONCLUSIONS
Knowledge about the covariance structure of CH 4 over a lactation is needed to model CH 4 for genetic evaluations under many different phenotyping strategies.We showed that the heritability was highest mid lactation (on average 0.17 ± 0.04), and genetic correlations between lactation stages were high (0.34 ± 0.36 to 0.91 ± 0.08).Permanent environmental correlations deviated greatly over a lactation and ranged between van Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION −0.73 ± 0.08 and 1.00 ± < 0.01, which highlights that it is most appropriate to model CH 4 c with a RR model including a random permanent environmental effect.With a large number of full lactation CH 4 recorded daughters, the reliability of bulls using the fixed versus RR model was similar.However, when shorter periods were recorded at the start and end of the lactation, the fixed regression model could achieve reliabilities for bulls of just 72% of the RR model reliability.Additionally, assuming the fixed model whereas the true (co)variance structure is reflected by the RR model, more than twice as long recording from the start of lactation was required to achieve maximum reliability for a bull.Thus, using a simplistic model could result in implementing too little recording, and lower genetic gains than predicted from the reliability.

van
Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION

vanFigure 1 .
Figure 1.The total number of records (left) and phenotypic pattern of weekly mean methane (CH 4 ) concentration (right) over a lactation as days in milk (DIM)

van
Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION

van
Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION

van
Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION

Figure 5 .
Figure 5.The additive genetic correlation (left) and within lactation permanent environmental correlation (right) between different stages of the lactation as days in milk (DIM) for weekly mean CH 4 c, modeled with a fourth order random regression model

vanFigure 6 .
Figure 6.The reliability of the selection index (r 2 ) for the random regression model using Legendre polynomials of the fourth order (Leg 4), the fixed regression repeatability model (Leg 0), and with using the selection index coefficients coming from the fixed regression repeatability model (b Leg 0) and the (co)variance matrices from Leg 4, on Bulls or Cows, with phenotypes being available for a cumulative number of days in milk (DIM) ranging from only 1 DIM to 274 DIM

van
Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION similar number of records throughout the lactation, or should be modeled using splines.

Table 1 .
van Breukelen et al.: GENETIC PARAMETER ESTIMATES FOR METHANE EMISSION The number of farms, number of recorded cows, and total number of records for visit and weekly mean methane CH 4 concentrations, after each step of data editing

Table 2 .
The mean, standard deviation (SD), minimum, maximum, and coefficient of variation (CV) of the weekly mean CH 4 c phenotype in parts per million (ppm), and after standardisation, recorded by two different versions of sniffers (v1.0 and v2.0) and as all measurements combined

Table A1 .
THE RELIABILITY OF THE SELECTION INDEX (R2) WITH PHENOTYPES ONLY BEING AVAILABLE AS SEQUENTIAL RANGES OF 21 DAYS IN MILK (DIM) (THREE WEEKS), WHEN USING THE PHENOTYPIC (CO)VARIANCE MATRIX (COWS), AND GENETIC (CO)VARIANCE MATRIX AS PHENOTYPES (BULLS) FOR THE RANDOM REGRESSION MODEL OF THE FOURTH ORDER (RR), WITH USING THE SELECTION INDEX COEFFICIENTS COMING FROM THE REPEATABILITY MODEL (B LEG 0) AND THE (CO)VARIANCE MATRICES FROM THE RR MODEL, AND THE DIFFERENCE BETWEEN THE TWO AS A PERCENTAGE