Lactation modeling and the effects of rotational crossbreeding on milk production traits and milk-spectra-predicted enteric methane emissions

Rotational crossbreeding has not been widely studied in relation to the enteric methane emissions of dairy cows, nor has the variation in emissions during lactation been modeled. Milk infrared spectra could be used to predict proxies of methane emissions in dairy cows. Therefore, the objective of this work was to study the effects of crossbreeding on the predicted infrared proxies of methane emissions and the variation in the latter during lactation. Milk samples were taken once from 1,059 cows reared in 2 herds, and infrared spectra of the milk were used to predict milk fat (mean ± SD; 3.79 ± 0.81%) and protein (3.68 ± 0.36%) concentrations, yield (21.4 ± 1.5 g/kg dry matter intake), methane intensity (14.2 ± 2.0 g/kg corrected milk), and daily methane production (358 ± 108 g/d). Of these cows, 620 were obtained from a 3-breed (Holstein, Montbé-liarde, and Viking Red) rotational mating system, and the rest were purebred Holsteins. Milk production data and methane traits were analyzed using a nonlinear model that included the fixed effects of herd, genetic group, and parity, and the 4 parameters ( a , b , c , and k ) of a lactation curve modeled using the Wilmink function. Milk infrared spectra were found to be useful for direct prediction of qualitative proxies, such as methane yield and intensity, but not quantitative traits, such as daily methane production, which appears to be better estimated (450 ± 125 g/d) by multiplying a measured daily milk yield by infrared-predicted methane intensity. Lactation modeling of methane traits showed daily methane production to have a zenith curve, similar to that of milk yield but with a delayed peak (53 vs. 37 d in milk), whereas methane intensity is characterized by an upward curve that increases rapidly during the first third of lactation and then slowly till the end of lactation (10.5 g/kg at 1 d in milk to 15.2 g/kg at 300 d in milk). However, lactation modeling was not useful in explaining methane yield, which is almost constant during lactation. Lastly, the methane yield and intensity of cows from 3-breed rotational crossbreeding are not greater, and their methane production is lower than that of purebred Holsteins (452 vs. 477 g/d). Given the greater longevity of crossbred cows, and their lower replacement rate, rotational crossbreeding could be a way of mitigating the environmental impact of milk production.


INTRODUCTION
In recent years, the study of traits associated with global climate change has become increasingly important (Moumen et al., 2016).Enteric methane emissions (EME) from ruminants are credited with being the most impactful source of greenhouse gases from food production at the world level (Roos et al., 2017;FAO, 2020), even though they have a short life in the atmosphere and their impact is much lower than previously thought (Place et al., 2022).Strategies have been developed for predicting traits related to EME of dairy cows by adjusting and calibrating measurements obtained from collected milk samples using Fourier-transform infrared (FTIR) spectroscopy.This approach is becoming increasingly common, and FTIR is currently one of the most routinely used technologies worldwide (de Haas et al., 2017).
The main metabolic link between the fermentation activity of the rumen microbiota and the infrared spectrum of milk is represented by the relationship between the proportion of VFA (Williams et al., 2019) produced in the rumen (acetate production is directly related to hydrogen and methane production) and the proportion of some fatty acids among the triglycerides of milk (van Gastelen and Dijkstra, 2016;Bougouin et al., 2019;Pitta et al., 2022) and other dairy products (Bittante and Bergamaschi, 2020).It is common knowledge that several milk fatty acids can be predicted by FTIR spectrometry (Engelke et al., 2018), and their proportions in milk vary during lactation in relation to the energy balance of cows and the mobilization or storage of their body reserves.The phenotypic and genetic relationships between some milk fatty acids and EME traits also vary during lactation (Vanrobays et al., 2016;Bittante and Cecchinato, 2020).
Prediction of EME proxies through FTIR is a promising technique, as it is much simpler and less expensive than methods based on the sampling and analysis of gases in respiration chambers, automatic feeders, automatic milking, and so on, and is not surrounded by ethical controversy over animal welfare (Garnsworthy et al., 2019).The predictions based on FTIR allow the EME proxies of a large number of cows to be estimated, and the effects of different factors, such as herd, breed, parity, and lactation stage, to be tested.
The EME are represented by traits with different meanings and properties (de Haas et al., 2017).The first trait to consider is methane production, which is the total daily production of methane per cow (dCH 4 , g/d).This is a quantitative trait affected by cow's feed intake, milk yield, and body size; therefore, dCH 4 is not a direct measure of the ecological efficiency but rather of milk production and feed intake (de Haas et al., 2017).Other traits that appear to be more closely related to the concept of efficiency (van Lingen et al., 2014) are methane yield, expressed as the methane emitted per kilogram of DM ingested by the cow (CH 4 /DMI, g/ kg), and methane intensity, expressed per kilogram of fat-and protein-corrected milk (CH 4 /CM, g/kg) or per kilogram of cheese obtained from the milk (Bittante et al., 2018).Finally, there are other traits measured by sniffer devices, which are based on the ratio of methane to carbon dioxide or the concentration of methane in the air (CH 4 /air; Sorg et al., 2018;Lee et al., 2022;Wang and Bovenhuis, 2019).All these traits are expected to be differently affected by the major sources of variation in milk yield, including genetics and lactation stage (Brito et al., 2018).
Variations in EME traits throughout lactation have not been investigated in depth, as trials using measured gas emissions are often relatively short.Methane emissions during lactation depend on feed intake and diet composition, which change during lactation in relationship to the variation of daily milk yield and nutrient requirements of lactating cows (de Ondarza and Tricarico, 2017;Beauchemin et al., 2022).Further knowledge could be obtained by comparing the pattern of variations during lactation of EME traits with those of milk yield.An interesting approach would be to test the use of parametric models developed to study milk yield throughout lactation to analyze EME traits (Bouallegue and M'Hamdi, 2019;Rekaya et al., 2000).The Wilmink model (Wilmink, 1987) has been widely used for modeling lactation curves, not because of its better statistical performance compared with other models, but because of its flexibility, simplicity (it is based on only 4 parameters), and ease with which function parameters can be interpreted, attributing them a physiological meaning (Macciotta et al., 2011).This last aspect could be very interesting if the physiological meanings used for explaining the development of milk production during lactation could be applied also to EME patterns.
Crossbreeding is increasingly used in the dairy sector (Guinan et al., 2019) due to its favorable effects on cow fertility (Malchiodi et al., 2014a;Hazel et al., 2020), longevity (Buckley et al., 2014;Clasen et al., 2017), and milk quality (Malchiodi et al., 2014b), which results in increased profitability despite reduced milk yield (Dezetter et al., 2017;Clasen et al., 2020;Hazel et al., 2021).The prevalent mating strategy in pasturebased dairy systems is a 2-breed alternating scheme using Holstein (HO) and Jersey bulls (Vance et al., 2012;Ferris et al., 2018).More common in indoor intensive dairy systems is the 3-breed rotational mating scheme, especially that using HO, Montbéliarde (MB), and Viking Red (VR) breeds.Following early research carried out by the University of Minnesota (Heins and Hansen, 2012;Hazel et al., 2017), this mating scheme is now also used in European countries (Malchiodi et al., 2014b;Dezetter et al., 2017;Saha et al., 2017), although its effects on the ecological footprint, and specifically on EMEs, is not well known.Therefore, the objectives of this study are (1) to describe and compare the lactation patterns of EME proxies predicted by FTIR spectrometry (CH 4 /DMI, CH 4 /CM, and dCH 4 ) using the models and physiological interpretation of the parameters tested for milk production traits; and (2) to assess the influence of either a 3-breed rotational crossbreeding system or HO purebred group on EME traits.

Experimental Design
The present study is part of the project titled, "Improving the ecological footprint of dairy farms through a three-way crossbreeding program" (ProCROSS/ Genesi project).The ProCROSS rotational crossbreeding scheme (ProCROSS, 2023)  The experimental design involved the selection of 2 large commercial herds belonging to the 2 main dairy farming systems of Northern Italy.The first farm (GP, 331 cows) is in the Lombardy region and follows the production regulations for Grana Padano cheese, the most important European Protected Designation of Origin (PDO) cheese in terms of quantity produced and revenue.The milk aimed at Grana Padano cheesemaking comes mainly from specialized intensive dairy herds, such as the GP herd, located in the lowlands of the northern part of the Po Valley, and from dairy cows generally fed TMR based on corn silage, concentrates, soybean meal, and roughages.The second farm (PR, 785 cows) follows the production regulations for Parmigiano Reggiano PDO cheese, the second most important Italian PDO cheese.The milk aimed at Parmigiano Reggiano production comes mainly from the Emilia-Romagna region, and is produced in specialized farms located on the plains and in the hills of the southern part of the Po Valley.Production regulations forbid the use of silages in the diet, so the feed consists of dry forages (especially alfalfa and meadow hay) and concentrates administered separately, or as TMR given dry or moistened with water.
Both herds have both been using the ProCROSS rotational crossbreeding mating scheme for at least 12 yr.The 2 farms have kept both genetics (HO and CR) throughout this period and continue to mate several purebred HO cows with MB or VR bulls every year, so they both have 5 crossbreeding generations, including first-and second-generation cows also among the younger animals.The GP farm followed the rotational sequence VR-MO-HO-VR-…, breed of sires, beginning initially from purebred HO cows; the PR farm followed the same sequence, but also the reverse sequence MO-VR-HO-MO-….Mating is not seasonal, but takes place all year round, and different parities were represented in both genetic types and the different crossbreeding generations.
On both farms, the HO and CR cows were kept together, managed as one herd, and fed the same TMR.As a consequence, the effects of the farm, genetic type, parity, and lactation stage could be considered substantially independent.

Animals, Milk Samples, and Acquisition of FTIR Spectra
Data for this study come from previous research aimed to investigate production traits and the cheese-making properties of milk of 1,116 purebred HO and crossbred cows (Saha et al., 2020), where details about crossbreeding generation and rotational sequence can be found.In total, 475 cows were HO and 641 were CR (185 belonging to the first generation of crossbreeding, 164 to the second, 222 to the third, and 70 to the fourth).The distribution of parity was 423 cows in first lactation, 343 in second lactation, and 350 in third and later lactations.The DIM (mean ± SD) was 155 ± 87 d.
Milk samples (100 mL) were collected once in all cows during the evening milking in 2 sampling dates (one per herd) according to the sampling protocol described by Saha et al. (2020).The FTIR spectra were acquired from milk samples and analyzed with a MilkoScan FT+ 6000 (Foss A/S).The instrument and its operations are validated and certified according to ISO 9622:2013/ IDF 141:2013(ISO, 2013), recently updated to ISO 21543:2020/IDF 201:2020(ISO, 2013, 2020;Niermöller and Holroyd, 2019).The complete spectrum of every milk sample analyzed was stored in the experimental database.A total of 1,060 absorbance values were recorded for each milk sample covering the infrared wavenumbers ranging from 5,000/cm (corresponding to a wavelength of 2.0 μm) in the near-infrared subdivision of the infrared area, through the mid-infrared to wave number 930/cm (corresponding to a wavelength of 10.8 μm) in the far-infrared subdivision (Bittante and Cecchinato, 2013).
For this study, data on milk yield and composition of all lactating cows with 5 to 365 DIM have been matched with FTIR spectra, resulting in an overall dataset of 1,059 cows: 439 HO (160, 132, and 147 for first, second, and third and more parity, respectively) and 620 CR (239,200,181 for first, second, and third and more parity, respectively).

Data Editing
A principal component analysis was performed on the FTIR spectra with Mahalanobis distances; samples with a probability level of <0.01 were considered outliers and removed from the analysis.A detailed description of this methodology is reported in Toledo-Alvarado et al. (2018).After editing, the dataset contained the complete spectra of 1,059 milk samples collected during routine test-day milk recording paired with the measured daily milk yield obtained from milk-recording data (dMY MR , kg/d), and the fat (%) and protein (%) concentrations predicted from the same spectra according to the official International Committee for Animal Recording procedure (Orlandini, 2020).From these traits, a daily fat-and protein-corrected milk yield (dCMY, kg/d) was calculated (CVB, 2008;van Lingen et al., 2014).Only records within x ± 3 standard deviations for each trait were subsequently used.The parity and DIM of each cow on the sampling date were also reported.

Prediction of EME
Direct predictions of EME traits from the milk FTIR spectra were made using the calibration equations described in a previous study (Bittante and Cipolat-Gotet, 2018).The EME prediction equations were obtained in the previous study by using, as reference, the values estimated from milk fatty acids according to a meta-analysis of trials carried out in respiration chambers (van Lingen et al., 2014).The predicted EME traits from infrared spectra were CH 4 /DMI, CH 4 /CM, and methane production per cow and per day (dCH 4-IR , g/d).
The predicting equations used were characterized by a calibration accuracy of 80%, 84%, and 69%, a cross-validation accuracy of 70%, 75%, and 60%, and a residual error of cross-validation of 1.18 g/kg, 1.17 g/ kg, and 86 g/d, respectively, for CH 4 /DMI, CH 4 /CM, and dCH 4-IR (Bittante and Cipolat-Gotet, 2018).Direct FTIR prediction of daily milk yield (dMY IR ) was also obtained for comparison with the measured dMY MR and with the direct FTIR-predicted dCH 4-IR .Lastly, an indirect prediction of dCH 4 (dCH 4-CMY , g/d) was made by multiplying the dCMY (kg/d) of each cow by its corresponding CH 4 /CM (g/kg).

Lactation Curve Modeling
The Wilmink's function initially used for modeling dMY MR (a zenith curve with negative values for the b and c parameters) was expanded to include positive and negative values (Amalfitano et al., 2021): where y t is the response of the trait (dMY MR , fat, protein, dCMY, CH 4 /DMI, CH 4 /CM, dCH 4-IR , dCH 4-CMY, and dMY IR ) at time t (corresponding to DIM); a is a parameter representing the potential value of the trait at the beginning of lactation, in the absence of adaptation; b is the parameter that quantifies the shortterm adaptation fraction soon after parturition; c is the parameter that describes the pattern of the curve after the peak, interpreted as the long-term variation (persistency of lactation); k is a parameter associated with the time of the peak of lactation, and interpreted as the speed of adaptation after parturition.
We used the PROC GLM of SAS (release 9.4, SAS Institute Inc., Cary, NC) to analyze each trait using DIM to define the initial parameters (a, b, and c), and with the value of the parameter k set to 0.05, as sug-gested by Wilmink (1987).These values were used as references for the preliminary statistical analysis to obtain suitable parameters to apply to the population.

Pattern and Derived Variables of the Curves
The patterns of the curves for each trait were classified according to a combination of the signs of parameters b and c (a and k are always positive).They were as follows (Amalfitano et al., 2021): • Zenith curve, when b < 0 and c < 0 (i.e., a curve increasing from calving to the maximum [peak] of lactation, then decreasing until the end); • Nadir curve, when b > 0 and c > 0 (i.e., a curve decreasing from calving to the minimum [negative peak] of lactation, then increasing until the end); • Upward curve, when b < 0 and c > 0 (i.e., a curve continually increasing from the beginning to the end of lactation [no peak]); and • Downward curve, when b > 0 and c < 0, i.e., a curve continually decreasing from the beginning to the end of lactation (no peak).
The detailed behavior of all the patterns can be consulted in Macciotta et al. (2005).
Other variables were also calculated from the 4 parameters of the lactation model to describe the lactation curves (Amalfitano et al., 2021) as follows: • Initial value, which is the sum of the maximum production potential (parameter a, long-term compartment) and the adaptation production (parameter b, short-term compartment) when DIM is equal to 0, [a+b]; • Peak value, which is the maximum (zenith curve) or minimum (nadir curve) value of the trait dur-

Statistical Analysis
All production and EME traits were analyzed with a nonlinear regression model using PROC TRANSREG and PROC NLIN of SAS (release 9.4; SAS Institute Inc.) with the Gauss-Newton algorithm.For each trait, the designed model was the following: where y hij is the response of the trait (dMY MR , fat, protein, dCMY, CH 4 /DMI, CH 4 /CM, dCH 4-IR , dCH 4-CMY , and dMY IR ); a, b, c, and k are the parameters of the lactation model; t is the DIM from 1 to 300; herd h is the fixed effect of the herd (2 levels: GP and PR herds); breed i is the fixed effect of cow genetic type (2 levels; HO, and CR); parity j is the fixed effect of the cows' parity (3 levels; first, second, and third and more parity); and e hij is the effect of the residual.

Milk Production and Quality Traits
Table 1 summarizes the descriptive statistics of the measured and FTIR-predicted traits analyzed in this study.The production traits from milk recording (dMY MR , and the fat and protein concentrations of the milk samples) were the subject, together with milk coagulation and cheesemaking properties, of a previous study on the effects of rotational crossbreeding (Saha et al., 2020), but in that previous study, lactation modeling was not tested.The results of the analysis of milk production and quality traits obtained are important to characterize the production environment of this study and to test the lactation modeling adopted here on traits that have a well-known lactation pattern (control traits) before moving on to the main objective of the study (i.e., EME proxies [experimental traits]).The distribution of dCMY was, as expected, very similar to that of the uncorrected daily milk yield (dMY MR ).
Effects of Farm, Parity, and Rotational Crossbreeding on Milk Traits.The 2 selected farms belong to 2 different farming systems: the GP herd belongs to the the Grana Padano system, with production mainly on intensive indoor farms in the lowlands using TMR based on corn silage, concentrates, and some roughages; and the PR herd belongs to the Parmigiano Reggiano system, a less intensive system with production more frequently in the hills, and using dry rations without silages.As expected (Mucchetti et al., 2017), the cows kept in the GP herd showed higher milk yield (Table 2), and the difference in milk yield between the 2 herds was close to 12%.In contrast, the milk of cows kept in PR herd had a 5.3% higher protein concentration than the milk of cows kept in GP herd.The effects of the parity of the cows in this study also confirm the expected increase in milk yield with advancing number of lactations (Table 2), and the small variation in quality traits.Lastly, the results obtained here confirm that the crossbred cows reared in intensive indoor farming systems tend to have a lower fluid milk yield (see dMY MR and dCMY in Table 2) and a higher fat and protein concentrations than the purebred HO (Hazel et al., 2020).

Effects of Lactation Modeling on Measured Milk
Traits.The lactation model is significant (P < 0.05) for all production traits, with at least 3 of the 4 parameters significantly different from the null value (Table 2).Consequently, Table 2 also reports the variables that describe the lactation curves derived from the model parameters.
The resulting lactation curves are depicted in Figure 1A, 1C, 1E, and 1G.The lactation modeling adopted here fully reflects the expectations regarding these traits.It is worth noting that the daily yield traits (dMY MR and dCMY) have a zenith curve with a positive peak in the first part of the lactation curve (negative short-term adaptation compartment b, and negative persistency coefficient c).On the contrary, the 2 qualitative traits (fat and protein concentration) of milk are characterized by a nadir curve with a negative peak (positive parameters b and c).The zenith curve shape of recorded daily production (dMY MR and dCMY) was well established by Wilmink's (1987) original research and subsequently by many other authors (e.g., Macciotta et al., 2005).
As is well known, the model is based on the estimation of 2 compartments, one long term and the other short term.The first represents lactation as an initial potential value (parameter a, theoretical intercept at DIM 0) and a linear slope (parameter c, persistency), which is negative for milk yield and depicts the decrease in milk yield toward the end of lactation.The second is a transitory compartment that represents the cow-udder effect of the adaptation after calving.Adaptation is represented as the difference between the theoretical and actual initial production (b parameter).This parameter is negative for milk yield because after parturition the cow needs to improve its feed intake capacity and the volume and productivity of the udder before reaching its potential milk production.The speed   of this adaptation process is represented by the fourth parameter, k, in the form of an instant rate constant causing the disappearance of the second compartment, and then the increase in milk yield, with advancing lactation.After parturition, the daily milk yield increases because the favorable effect of adaptation is greater than the unfavorable effect of the decrease in persistency.The 2 opposite effects are equal at the peak of lactation, whereas the decrease in persistency prevails after the peak till the end of lactation (zenith curve).
The first, long-term, compartment parameters, a and c, are greater in this study, and the average daily milk yield (dMY) is almost twice that reported by Wilmink (1987) and reflect the improvement of cow genetics, feeding, and management.Regarding the second, shortterm compartment, it should be noted that due to the difficulty in simultaneously calculating all 4 parameters of the model, Wilmink (1987) assumed a priori that k = 0.05, as did many other authors in subsequent studies.The a priori assumption of a k value consequently determines the estimated value of the b parameter.As our modeling was nonlinear, we directly estimated the k value, that resulted much larger (k = 0.103, Table 2) than Wilmink's one, and consequently also b parameter was greater.However, the result is that we estimated the peak would be reached at 37 DIM, which is in the range of estimates by Wilmink (1987).In a previous study on Brown Swiss cows using the same modeling, the peak was estimated at 32 DIM (Amalfitano et al., 2021).
Milk quality traits are well known to be characterized by a curve that is the inverse of the production traits curve, and this was confirmed in this study.Therefore, the same lactation model could be adapted to describe the variation in the fat and protein concentrations of milk, provided that the sign of b and c parameters can change.In the case of fat and protein concentration, as expected, both the persistency (c) and adaptation (b) parameters are positive.Then the shape of the curve is inverse of that of milk yield, with a high initial actual value decreasing till a minimum value and then increasing till the end of lactation (nadir curve).Few studies have used this lactation modeling for quality traits, and among these, the previous study of Amalfitano et al. (2021) on Brown Swiss cows yielded similar model parameters.

Effects of the Fixed Effects and Lactation Modeling on dMY IR .
As a further control trait to evaluate the EME traits, and particularly dCH 4 , we also predicted dMY IR .The results were very different from those for the corresponding trait dMY MR .First of all, average production was much lower for the FTIRpredicted trait than for the measured trait (24.6 kg/d for dMY IR vs. 32.7 for dMY MR ; Table 1).It is worth noting that 24.6 kg/d is very similar to the average dMY of the dataset used to calibrate and validate the prediction equations (Bittante and Cipolat-Gotet, 2018).Moreover, with regard to the effects of the 3 fixed factors considered, dMY IR yielded more variable results than dMY MR .With respect to the measured production (dMY MR ), Table 2 shows that the prediction of milk yield based only on milk FTIR spectra (dMY IR ) confirms the higher level of the GP herd over the PR herd, but the higher level of the HO cows over the CR cows was not, and the effect of parity reversed, with an unrealistic prediction of decreasing production from the first lactation to maturity.The lactation modeling for dMY IR differed little from that for dMY MR , but the values were much lower and the pattern flatter (Figure 2).This means that FTIR predictions of a quantitative trait, such as dMY, are unreliable, even though the spectra capture some of the variations in production, especially when they are (negatively) correlated with quality traits (such as in early and late lactation).On the other side, there is no theoretical justification linking a qualitative picture of milk, such as its FTIR spectrum, to the quantity of milk produced.

EME Traits
According to the main aims of this study, an important topic of discussion is the use of EME proxies predicted from FTIR spectra acquired from milk samples (Negussie et al., 2017;Shetty et al., 2017;Bittante et al., 2020).The 4 EME traits had very different sources of variation and patterns during lactation according to their nature and the estimation method.
Methane Yield.The average CH 4 /DMI of HO and CR cows in this study (21.4 g/kg; Table 1) is within the interval (20.1 ± 4.3) obtained by Hristov et al. (2018) in an extensive review covering different experimental methods and dairy systems.The same authors reported that of the various EME traits, this one had the lowest coefficient of variation across studies (21.6%).In a previous study on Brown Swiss cows, Bittante et al. (2018) obtained an average CH 4 /DMI of 21.3 g/kg, whereas van Lingen et al. ( 2014) reported an average of 21.5 g/ kg in their meta-analysis of several trials with HO cows.It is worth noting from Table 1 that the coefficient of variation of CH 4 /DMI (7.0%) is the lowest among all the traits studied.
Table 3 summarizes the results of the nonlinear statistical analysis of the data on the EME traits directly predicted from FTIR spectra (CH 4 /DM, CH 4 /CM, and dCH 4-IR ), and from the milk-recording data and FTIR spectra (dCH 4-CMY ).Methane yield (CH 4 /DMI) was affected only by herd, whereas breed, parity, and lactation stage did not affect this trait.The PR farm could be expected to have a higher CH 4 /DMI due to the greater use of forages in the diets, while there is no apparent reason to expect differences between the HO and CR cows for this trait, as on each farm they received the same diet.Given the same diet as in the present study, the fermentation pattern in the rumen (microbiota composition, proportions of VFA, and so on) is not expected to differ between HO and CR cows.Therefore, potential differences in the rumen fermentation patterns of the 2 genetic types could arise from differences in DMI, whereas similar BW and modest differences in milk yield are not expected to have much effect on this trait.It is worth noting that no scientific article the authors are aware of reports the opportunity of different prediction equations for milk from cows of different dairy breeds reared in intensive farming systems.The same factors (the same diet and modest differences in DMI) could explain the absence of any significant effect of cow parity on CH 4 /DMI.
The pattern of this trait during lactation as depicted in the lactation modeling is almost flat (Figure 1B, 1D, 1F, and 1H).In fact, none of the lactation model parameters of this trait was significant (Table 3), and so the derived traits were not calculated.Some small variations could be a consequence of variations in the amount of DM ingested daily by the cows, which is expected to increase until about mid lactation and to decrease thereafter (Mäntysaari et al., 2012).In previous studies on Brown Swiss cows, Bittante and Cipolat-Gotet (2018) and Bittante et al. (2018) found that DIM classes affected the CH 4 /DMI trait, but the maximum difference among the least squares means values was only 0.7 g of methane produced per kilogram of DM ingested by cows.
Methane Intensity.The average CH 4 /CM found in this study (14.2 g/kg; Table 1) is within the interval (15.6 ± 8.5) obtained in the articles reviewed by Hristov et al. (2018), who also found that this trait had the highest coefficient of variation across the different studies (54.4%).This average value is very similar to those previously found by Bittante et al. (2018) for Brown Swiss cows and by van Lingen et al. ( 2014) for HO cows.It is probable that this trait has an inverse relationship with milk production (denominator of the ratio), so the lower values found in the GP herd and in mature cows (Table 3) are expected.Similarly, a marginally lower value for HO cows is expected compared with CR cows (no significant difference).
Regarding the lactation curves, in a previous study on milk protein fractions (Amalfitano et al., 2021), where the model used for dMY was adapted to several other milk traits, 2 other curve patterns without peaks were obtained: an upward curve with negative b and positive c parameters and a downward curve with positive b and negative c parameters.Here (Table 3), CH 4 / CM exhibited an upward pattern, increasing rapidly at the beginning of lactation and more slowly toward the end (Figure 1).This curve is very similar to the polynomial trend found in the previous work on Brown Swiss cows of Bittante et al. (2018).The curve intercept in that work (10.9 g/kg) is almost identical to that in the present study in HO and CR cows (10.5 g/kg), and the values calculated after 300 DIM are also very similar (15.5 and 15.2 g/kg, respectively).Obviously, the linear, quadratic, and cubic slopes of the polynomial curve cannot be interpreted from a physiological point of view.Moreover, after 300 DIM the value obtained from the polynomial increased very rapidly toward infinity, whereas in the lactation modeling presented here it was almost stable.
The CH 4 /CM is expected to be low after parturition with a tendency to increase thereafter, rapidly at beginning and slowly toward the end of lactation (Kandel et al., 2017).In fact, we expect that dCH 4 will reflect the amount of feed ingested by the cow.However, at the beginning of lactation, milk production is much higher than allowed by the cow's DMI (negative energy balance), and part of the milk is produced using energy and nutrients mobilized from the cow's body reserves (Roche et al., 2009).The use of these reserves does not affect the actual production of methane in the rumen, but reflects part of the methane emitted by the cow before parturition when the body reserves were stored.The short-term compartment of the CH 4 /CM lactation curve can then be interpreted as the quantity of milk produced at the beginning of lactation from body reserves without producing methane.At the same time, the positive persistency of CH 4 /CM in the long-term compartment could be interpreted as the cow's need to eat more than required for maintenance and lactation to restore the body reserves and support pregnancy.The values of the variables adaptation (−25%) and persistency (+8%) shown in Table 3 are therefore a numerical indication of these opposing phenomena.
Methane Production Directly or Indirectly Predicted from FTIR Spectra.The average dCH 4-IR found in this study (358 g/d) is identical to the average value (358 ± 105 g/d) reviewed by Hristov et al. (2018), and almost identical to the average value found in an earlier study by Bittante and Cipolat-Gotet (2018) on Brown Swiss cows (357 g/d), but lower than the value for HO cows reported by van Lingen et al. (2014).Hristov et al. (2018) also showed this trait to have an intermediate coefficient of variation across studies (29.2%).
It is worth noting, however, that by dividing this value of dCH 4-IR by dCMY (Table 1), we obtain an indirect estimate of CH 4 /CM of only 10.2 g/kg, much lower, as previously seen, than the values obtained in the present study through FTIR prediction for HO and CR cows (14.2 g/kg), in the Hristov et al. (2018) review of different breeds and methods (15.6 g/kg), in van Lingen et al. (2014) review of HO cows in respiratory chambers, and for Brown Swiss cows (14.2 g/kg) using the milk fatty acid profile of milk (Bittante et al., 2018).
As previously seen, dCH 4 is expected to be closely related to DMI, which in turn is expected to reflect the   cow's energy requirement (except at the beginning of lactation).Milk energy is by far the largest and most variable component of the total energy requirement of a dairy cow.We expect a close relationship between dMY (or rather dCMY) and dCH 4 (Martínez-Marín et al., 2023).Summarizing the trials carried out in respiration chambers, Hristov et al. (2018) 2 and 3).This means that the inability of FTIR spectrometry in capturing information on the quantity of milk produced daily, is paralleled by its inability in capturing the quantity of methane produced daily (both quantities are unreliably low, but their ratio seems correct).
In the present study, dCH 4-IR was found to be higher for cows on the PR farm, cows obtained from rotational crossbreeding (not significantly) and those in their first and second lactations (Table 3), which is similar to dMY IR , but contrary to expectations based on dMY MR (Hristov et al., 2018).
Regarding the lactation modeling, we obtained an upward pattern (Figure 1) and this is also contrary to the expectation of a zenith curve simulating the classical lactation curve and DMI curve.In the previous study on Brown Swiss cows, Bittante et al. (2018) obtained a cubic polynomial curve characterized by a growing phase during the first 3 mo of lactation followed by a decrease till end of lactation, not very different than the expected zenith curve.
It is evident that direct prediction of a quantitative trait, such as dCH 4-IR , using FTIR spectra alone is as unreliable as prediction of dMY IR .A similar conclusion was also reached recently by Denninger et al. (2020).This is also confirmed by the much lower R 2 of the equation predicting this trait compared with the equations predicting qualitative EME traits (Bittante and Cipolat-Gotet, 2018;Coppa et al., 2022).Vanlierde et al. (2015) also found that dCH4 predictions using FTIR spectra alone produced a lactation pattern (nadir-type curve), contrary to expectations.The inclusion of DIM information (linear and quadratic Legendre polynomials) in the prediction equations did not increase the R 2 of cross-validation, but yielded a lactation pattern according to expectations (zenith-type curve) and increased the robustness of the predictions when adopted at the population level.The need to combine spectra information with a variable expressing quantitative information was further demonstrated by Vanlierde et al. (2021), who concluded that the explanatory models based on FTIR spectra and DIM were improved by also including dMY, parity and breed information, and also reported that when these 3 variables were included in the predictive equations the predictions of dCH 4 were in closer agreement with the literature.
In contrast, it is evident that after adding DIM, dMY, parity, and breed to the prediction equation, the contribution of the FTIR spectra to the accuracy of the results becomes less certain.
Methane Production Indirectly Predicted from FTIR Spectra and Measured dCMY.In the previous study on the development of FTIR-based EME prediction equations (Bittante and Cipolat-Gotet, 2018), it was concluded that direct FTIR predictions of CH 4 / DMI and CH 4 /CM yielded informative results, but not dCH 4 .We suggested that a more informative result could be achieved by estimating methane production by simply multiplying the measured dCMY from milk recording by the FTIR-predicted CH 4 /CM.This is why we calculated this EME trait (dCH 4-MR ) here.
The first difference between dCH 4-MR and dCH 4-IR (Table 1) is that the average value of the former is much larger than that of the latter (450 vs. 358 g/d), and is closer to that found for HO cows with a high dMY, as with ours.The second difference is that the fixed effects included in the dCH 4-MR statistical model (Table 3) were all in agreement with expectations (whereas those of dCH 4-IR were all the inverse): GP > PR herd, HO > CR cows, and third > second > first lactation.The third difference is that, as expected, dCH 4-MR showed a zenith curve and not an upward curve as dCH 4-IR (Figure 1).

Effects of 3-Breed Rotational Crossbreeding on EME Traits
Some information is available on the effects on some EME traits of cow breed (Lassen et al., 2012;Vanlierde et al., 2021) and of HO-Jersey crossbreds (Xue et al., 2011;Hynes et al., 2016), while the authors are not aware of any information on the effects of rotational crossbreeding on the ecological footprint of the intensive dairy production chain.The EME is only one aspect of the ecological footprint, and EME during lactation is only a fraction of this aspect.Nonetheless, EME is known to be the most important part of the ecological footprint of ruminants, and lactation the most important segment of the cow's lifetime (Berton et al., 2020).
It is generally accepted that the environmental impact is reduced, first of all, by increasing produc-tion efficiency, and that efficiency is related to animal productivity (Kandel et al., 2017).Crossbreeding is expected to reduce the daily milk yield of cows, so it can be expected to reduce feed intake, but also production efficiency.Possible crossbreeding expectations for EME traits could include the following (1) no effect on CH 4 /DMI, if the diet is the same as that given to HO cows; (2) improvement (decrease) in dCH 4 , due to an expected lower intake of feed for milk production (but a constant intake for maintenance and pregnancy); and (3) worsening (increase in) CH 4 /CM, because the numerator (dCH 4 ) is expected to decrease less than the denominator (dCMY).
In the present study, the first expectation is confirmed because the LSM of the purebred HO and the 3-breed CR for CH 4 /DM were not significantly different.The second expectation is also confirmed because HO cows have a higher dCH 4-CMY than CR cows.The third expectation was not confirmed because the HO and CR cows presented very similar LSMs for CH 4 / CM.This better-than-expected CH 4 /CM in CR cows could be due to the better feed efficiency of these cows, which can compensate for their slightly lower dCMY (Pereira et al., 2022).Shonka-Martin et al. (2019) found that VR-MB-HO rotational crossbred cows are both nutritionally and economically more efficient than purebred HO cows.
Therefore, we can conclude from this study that rotational crossbreeding does not appear to increase the production of methane per kilogram of milk produced during lactation, although it should be remembered that the ecological footprint of dairy production also includes the EME of dry cows and replacement heifers.Rotational crossbreeding has proven to be particularly effective in improving the fertility (Malchiodi et al., 2014a) and longevity of cows, reducing the risk of culling (Piazza et al., 2023a).This means that the (environmental) costs of heifer production are mitigated by a longer productive lifetime and overall milk production in the case of CR cows (Piazza et al., 2023b).Taking a life cycle assessment approach (Beauchemin et al., 2022) at the individual cow, we found that overall methane production (including not only lactating cows, but also heifers and dry cows) per kilogram of milk produced in the career is lower in CR than in purebred HO cows (Piazza et al., 2022).

CONCLUSIONS
Milk FTIR spectra have proven to be useful for the direct prediction of proxies of qualitative EME traits, such as CH 4 /DMI and CH 4 /CM, but not for quantitative traits, such as dCH 4 .Estimation of the daily production of methane seems better done by multiplying a measured dMCY by a FTIR-predicted CH 4 /CM.Lactation modeling of EME traits confirms that dCH 4-MR has a zenith curve, similar to that of milk yield but with a delayed peak, as expected for the cows' DMI.The CH 4 / CM is characterized by an upward curve that increases rapidly after calving in parallel with the reduction in the negative energy balance, then slowly until the end of lactation, in parallel with the recovery of body reserves and the progression of pregnancy.In contrast, lactation modeling was not useful for explaining CH 4 / DMI, which is almost constant during lactation.Lastly, 3-breed rotational crossbreeding involving HO, MO, and VR did not unfavorably affect CH 4 /DMI and CH 4 / CM during lactation, and reduced dCH 4 compared with purebred HO.Taking into account the greater longevity of CR cows and their lower replacement rate, rotational crossbreeding could be seen as a way to mitigate the environmental impact of milk production.
which is the number of days from parturition to reach the peak of lactation, [-(1/k) × ln(c/(k × b))]; • Final value, which is the production value when lactation ends, calculated as the sum of the maximum potential (parameter a) and the persistency (parameter c) multiplied by the defined final DIM (DIM f ; for this study the limit was 300 DIM), [a-c × DIM f ]; • Adaptation fraction, which is the proportion of the adaptation in relation to the maximum potential expressed as a percentage, [b/a × 100]; and • Persistency fraction, which is the percentage of the production variation up to the end of lactation in relation to the maximum production potential, [(c × DIM f )/a × 100].Martínez-Marín et al.: EFFECTS OF ROTATIONAL CROSSBREEDING

2GP=
Grana Padano herd (using corn silage), PR = Parmigiano Reggiano herd (no silages).HO = Holsteins; CR = 3-breed rotational crossbreds.3 Trait obtained the measured dMY MR corrected for FTIR-predicted fat and protein concentrations.4Parameters: a = potential value of trait at the beginning of lactation, in the absence of adaptation; b = short-term adaptation fraction soon after parturition; c = pattern of curve after the peak, interpreted as longterm variation (persistency of lactation); k = time of lactation peak, interpreted as speed of adaptation after parturition.*P < 0.05; **P < 0.01; ***P < 0.001.

Figure 1 .
Figure 1.Lactation patterns of milk production traits (panels A, C, E, and G) and enteric methane emission traits (solid lines, panels B, D, F, and H).Dashed lines refer to the long-term compartment (y = a + c, where y is the independent variable, a is the potential value at the beginning of lactation, and c is the persistency during lactation).The difference between the solid line and the dashed line represents the effect of the short-term adaptation compartment depending on b (adaptation fraction at the beginning of lactation) and k (speed of adaptation) parameters.dMY MR = measured daily milk yield (kg/d); dCMY = daily corrected milk yield (kg/d); CH 4 /DMI = methane yield per unit of DMI (g/kg); CH 4 /CM = methane intensity per kilogram of corrected milk (g/kg); dCH 4-IR = daily methane production predicted by infrared spectroscopy (g/d); dCH 4-CMY = daily methane production obtained indirectly by multiplying dCMY by FTIR-predicted CH 4 /CM (g/d).FTIR = Fourier-transform infrared spectroscopy.

Figure 2 .
Figure 2. Comparison between the lactation patterns of milk yield measured during milk recording (dMY MR , in blue) and milk yield predicted by Fourier-transform infrared spectroscopy milk spectra (dMY IR , in green).Dashed lines refer to the long-term compartment (y = a + c, where y is the independent variable, a is the potential value at the beginning of lactation, and c is the persistency during lactation).The difference between the solid line and the dashed one represents the effect of the short-term adaptation compartment depending on b (adaptation fraction at the beginning of lactation) and k (speed of adaptation) parameters.
is managed by the leading suppliers of MB genetics (Coopex Montbéliarde, Roulans, France) and by the Scandinavian organization leading in VR genetics (Viking Genetics, Randers, Den-Martínez-Marín et al.: EFFECTS OF ROTATIONAL CROSSBREEDING 1487 mark).Because no human or animal subjects were used for the current analysis, approval by an Institutional Animal Care and Use Committee or Institutional Review Board was not required.

Table 1 .
Martínez-Marín et al.: EFFECTS OF ROTATIONAL CROSSBREEDING Descriptive statistics of milk production traits and enteric methane emission (EME) traits 1 tion per cow per day predicted directly from FTIR (g/d); dCH 4-CMY = daily methane production obtained indirectly by multiplying dCMY by FTIR-predicted CH 4 /CM (g/d); FTIR = Fourier-transform infrared spectroscopy.2 Trait obtained from the measured dMY MR corrected for FTIR-predicted fat and protein concentrations.3 Trait obtained by multiplying the mixed dCMY by FTIR-predicted CH 4 /CM.

Table 2 .
Martínez-Marín et al.: EFFECTS OF ROTATIONAL CROSSBREEDING Significance levels and least squares means of the fixed effects (herd, breed, and parity), and lactation modeling of the milk production traits 1 1 dMY MR = measured daily milk yield (kg/d); dMY IR = Fourier-transform infrared spectroscopy (FTIR) prediction of daily milk yield; dCMY = daily corrected milk yield (kg/d).

Table 3 .
Significance levels and least squares means of the fixed effects (herd, breed, and parity), and lactation modeling of the enteric methane emission (EME) proxies 1 IR by dMY -CMY is unreliable (10.2 g/kg), the CH 4 /CM became fully congruent (14.5 g/kg) when we divide the dCH 4-IR by dMY IR (Tables