Introduction

Attempts to reduce greenhouse gases (GHG) are aimed at mitigating the causes of climate change. On an international level, GHG inventories measure progress made by different countries in reducing emissions in production sectors. In 2010, Mexico’s GHG emissions totaled 748,252 Gg carbon dioxide equivalents - 1 Gg is equal to 10⁹ g¹. The overall mass of gas emitted into the atmosphere and the radiative forcing effect of each GHG in this mass determines the importance of greenhouse gases such as CO₂, CH₄, N₂O, and O₃, among others. Methane (CH₄) deserves particular attention among the GHG because its radiative forcing effect is 25 to 34 times greater than that of CO₂, that is, 1 g CH₄ equals at least 25 g Ceq².

The role of livestock production in GHG emissions is actively debated, but clearly significant^3,4,5. In Mexico, for example, total CH₄ emissions from enteric fermentation accounted for 20.1 % methane emissions while fugitive emissions from petroleum production represented most of the remainder⁶. Enteric and manure-related CH₄ is responsible for 75 to 95 % of total GHG emissions at dairy farms in 38 countries⁷. Worldwide, Mexico was in eleventh place as far as enteric CH₄ emissions (22,244 Gg) in 1998⁸. By 2010, however, emission from enteric fermentation was 37,961 Gg Ceq and that from feces management was 7,553 Gg Ceq, which together accounted for almost half the 92,184 Gg Ceq total CH₄ emissions from the agricultural sector¹. However, other studies have generated disparate figures for the total amount of CH₄ emitted. One study reported emissions of 144,000 Gg CH₄ in 1970 and 270,000 Gg CH₄ in 2010 for the dairy sector in Mexico⁹, whereas another study reported emissions in 1995 to have been 285 Gg CH₄¹⁰. This wide range of figures highlights the need to reduce uncertainty in emissions estimations by studying emission factors and intensities.

Signatory countries of the United Nations Framework Convention on Climate Change (UNFCCC) agreed to annually report GHG emission inventories following Intergovernmental Panel on Climate Change (IPCC) guidelines, using either of three Tier or precision levels¹¹. For CH₄ emissions originating in enteric fermentation, Tier I considers activity levels and emission factors (EF) per head of livestock; the former are taken from livestock inventory data and the latter from default values specified for continental regions and productive systems. The CH₄ emissions level is then calculated by multiplying the activity level by the EF. In Tier II, the EF are calculated from net energy intake for maintenance, activity, growth, gestation and lactation. Tier III determines the EF using direct GHG emissions measurement methodologies.

The Tier I methodology, based on calculation of EF per head, discourages changes in productivity as a mitigation strategy while GHG emissions per product unit decrease as productivity increases¹². Emission intensity varies among producers in response to agroecological conditions, grazing and rearing practices, feed regimes and feed type. The concept of emission intensity is useful in evaluating mitigation strategies at the productive system level and for guiding feed intake patterns¹³.

Some data are not normally used in GHG inventory generation. Integrating them into inventory calculation can improve the cost:benefit ratio of inventory generation compared to other options such as running experiments, designing new measurement procedures or modifying existing methodologies or equations. For dairy herds in Mexico, an excellent source of this type of data is maintained by the Holstein Association of Mexico (Asociación Holstein de México - AHM), which collects, processes, certifies and publishes data on registered dairy herds. The AHM generates production indicators aimed at improving cattle genetics and archives genealogical data for production in Holstein-Friesian cows. These records are sufficient to allow calculation of enteric fermentation CH₄ emissions using the Tier II methodology.

The objective of the present study was to calculate the enteric CH₄ emissions for lactation and dry periods of dairy cows using data from 2007 in the state of Querétaro, Mexico. In this same year, the National Institute of Statistics, Geography and Data Processing (Instituto Nacional de Estadística, Geografía e Informática - INEGI) conducted the Agricultural, Livestock and Forestry Census. The Tier II methodology was applied using AHM data representing confined, intensive system dairy herds in the state of Querétaro, Mexico.

Material and methods

The lactation data used in the calculations included herd and cow identification; date of birth; calving date; lactation number; lactation end date; milk yield (kg) and weighing date; and milk fat and protein percentage. The AHM laboratory is certified by Quality Certification Services Inc. and runs fat and protein measurements using a mid infrared analyzer (Bentley 2000, Bentley Instruments Inc., USA). Each lactation record contained up to ten weighing dates at an average interval of 33 d. A total of 11,092 records were used in the analysis, including those with at least one lactation day in 2007, lactations longer than 100 d and shorter than 500 d. Limiting the analysis to lactations of less than 500 d results in a reasonable compromise between cows with lactations of 305 d and those with longer lactations, as only 5 % of cows lactate longer than 500 d¹⁴, but 55 % lactate for more than 305 d¹⁵.

Milk yield determines a significant proportion of a cow’s energy requirements; daily milk yield can vary during lactation and by lactation number. The incomplete gamma function adequately represents the lactation curve characterized by an ascending phase, a peak and a gradual decline¹⁶; indeed, it has an error variance equal to or less than other formulas (e.g. Gaines, Rook, Dijkstra or Pollott) used to represent lactation curves in herds in central Mexico¹⁷. Given the relative predictive capacity of complex lactation curve models, those like the incomplete gamma function are recommended, unless mechanistic details are needed for extended lactation curves (> 800 d)¹⁵. The formula used here was:

Y _t = a · t^b · e^{- ct}

where Y_t is milk yield on the t^th day of lactation; e is the base of the natural logarithms; and a, b and c are the parameters defining the curves to fit. Yield immediately after calving is represented by a; the rate of change from initial production to maximum yield is b; and c is the rate of change from maximum yield to end of lactation. Persistence (s) is defined as the point at which peak milk yield is maintained; it is calculated as “(b+1) · ln(c)”. Milk yield per 305 d lactation (Y) was calculated as the sum of Y_t to d 305. Maximum milk yield (Y_max) was determined with the formula a · (b/c) ^b · e^-b , and mean duration of lactation (DL) was the average number of days in each lactation.

The lactation curve model was generated using the package Table Curve 2D (Systat Inc.). Significance was confirmed with the F test for goodness of fit error, normality and homogeneous variance of residuals. Only significant models that met the assumptions of normality and homoscedasticity were considered in the analysis. Residuals independence was confirmed with the autocorrelation function to the 10^th order, and a t-test using the Minitab v. 16 pacakege (Minitab Inc.). Modeled variables were milk yield (Y_t) and fat percentage (G_t). Cows within each herd were subsamples, while each herd was a replicate. This structure allowed us to infer about the herds, and confirm goodness of fit considering the herd as a characteristic management condition of each group of animals. Repetitions were n = 16, 17, 16, 16, 13, 14 for the first to the sixth lactation (L1 to L6). The sixth lactation included records for cows with six or more lactations. Parameter estimator means for the incomplete gamma function, DL, s, Y, Y_max and peak production were compared with a Tukey-Kramer test for imbalanced models¹⁸, with lactation number as the explanatory variable. The contrast between primiparous and multiparous cows was also tested. The confidence level in all analyses was 95% and analyses were run with the JMP v7 package (SAS Institute Inc.).

Methane (CH₄) emissions were calculated for each herd and lactation number using the Tier II methodology¹¹. Emissions from enteric fermentation were calculated in terms of gross energy (GE) content and the conversion rate of energy to methane (Y_m). The only feed property considered was energy digestibility:

Emissions= EF · Population

Where Emissions is the release of CH₄ from enteric fermentation (kg CH₄ yr^-1); EF is the emission factors corresponding to a specific population (kg CH₄ cow-1 yr^-1); and Population is the number of cows. The EF is calculated as:

EF= (GE · Y_m · 365 d) / (55.65 MJ / kg CH₄)

Where Y_m is the CH₄ conversion rate from GE with a value of 0.06 equivalent to 60 kJ MJ^-1(11. However, Y_m can vary in response to milk yield and feed properties. For example, Y_m is 61 kJMJ^-1 in high yielding cows (7,200 kg L cow^-1 yr^-1) in Germany¹⁹, and in high-yield cows in the United States Y_m was 5.63 % of GE²⁰. Absorption of GE in MJ cow^-1 d^-1 is derived from the sum of net energy needs and feed digestibility characteristics (DE= digestible energy, etc.):

GE= {[(E_m + E_a + E_L + E_p) / (E_ma / DE)]} / (DE/100)

Net maintenance energy (E_m) in MJ d^-1 is: E_m= Cf_i · (weight)^0.75

Where Cf_i= 0.335, based on the coefficients for calculating E_m¹¹. For the purpose of calculating E_m, the body weight corresponding to L1, L2 and L3 were calculated from a sample of 165 cows from eleven corrals. The resulting averages were 526 kg for L1, 590 kg for L2 and 640 kg for L3; the weight 680 kg was used for lactations 4, 5 and 6 since it represents large cows²¹.

Net energy for milk synthesis (E_M), in MJ d^-1, is:

E_m = Y_t · (1.47 + 0.40 · G_t)

Net energy for pregnancy (E_M), in MJ d^-1 is:

E_p= C_p · E_m

where C_p= 0.10, which is the pregnancy coefficient.

The ratio between maintenance energy in a diet and digestible energy (DE) intake (E_ma/ DE), in MJ d^-1 is:

E_ma / DE= 1.123 - (4.092 · 10^-3 · DE) + [1.126 · 10^-5 · (DE)²] - (25.4/DE)

Feed energy concentration was expressed as digestible energy (DE), calculated in Mcal kg^-1 and converted to MJ kg^-1(21:

DE= 0.1847371 · TDN (%)

Where TDN is total digestible nutrients. This is equal to 68 % of TDN in the diet of high-yield, large (average weight = 680 kg) cows, and includes a 3 % reduction in digestion efficiency in the TDN value when feeding at three times the maintenance level²¹. The 12.56 MJ kg^-1 value used for DE is similar to the 18.45 MJ kg^-1 default value for GE concentration²². Only one energy concentration value was used for all lactations, and changes in live weight during lactation were not considered. Daily CH₄ emissions were calculated up to d 305 to generate an EF comparable to the 305 d adjusted yield (Y), and also until the total DL to obtain the EF during lactation. Emission intensity per unit produced (EI, CH₄ kg^-1 milk) was calculated per lactation by adding the CH₄ produced in all the herds during all the lactation days and dividing the result by the number of lactating cows on said days. This figure was divided by average milk yield in the lactation period.

Results and discussion

Means of the lactation curve parameters were different between lactation numbers (Table 1). In terms of milk yield, differences were found between primiparous and multiparous cows in the a and b parameters of the incomplete gamma function (P<0.0001), s (P=0.05) and Ymax (P=0.017); no differences were present for c, DL or Y. Lactation persistence (s) was high due to a low rate of decrease in production after Ymax which generated relatively flat curves. In intensive system dairy herds, the a and b parameters differ between the first, second and third lactations¹⁷. These variations are probably better explained through analysis of seasonal variability, herd management conditions or genotype.

The number of modeled lactations was highest for L1 and then decreased thereafter (Table 2); this is probably an artifact of the voluntary early culling typical of intensive production dairy herds^15,23. The incomplete gamma function represented most (89.8 %) of the studied lactations, but all the L1 and L2 lactations had good fit (P<0.05). Lack of fit appeared when records were incomplete, or when the number of records increased as lactation number progressed in some herds. In contrast, the models were not significant when the curve was flat and the b and c values were near zero (Table 1). The capacity of the incomplete gamma function to represent L1 or L2 agreed with reports for dairy herds in central Mexico in which a good fit was attained using this formula for the first, second and third lactations¹⁷.

The long average DL (365 d) indicates the presence of extended lactations in the analyzed database (Table 2). Use of extended lactations is a management strategy that capitalizing on high lactation persistence and the financially profitable yield towards the end of a lactation; however, it can also be a consequence of reproductive failure¹⁵. Extended and induced lactations are common in Mexican dairy herds, and are maintained by using bovine somatotropin²³. This is important to consider when calculating the annual GHG emissions inventory for a state or country since methane emissions will be greater in the absence of a dry period in any given year.

Values for 305 d milk yield (Y) were similar between the analyzed lactations despite differences identified between the a and b estimates. Average Y was 9,985 kg in the database modeled with the incomplete gamma function. In terms of the variation coefficient, variation in Y in the database (28 %) was acceptable. However, as part of the uncertainty in the calculations, the standard error values must be included in the input parameters when calculating methane emissions inventories²⁴. The 95% confidence interval for Y in the present data ranged from 8,550 to 11,419 kg, which differs from the 7,078 L (7,312 kg corrected to 3.5 % fat) reported for Y in intensive diary production systems in Mexico for operations located in the states of Aguascalientes, Hidalgo, Guanajuato, Jalisco, Michoacán, Puebla and Querétaro²⁵.

Using the average DL figure of 365 d, the 95% confidence level was 11,061 to 12,793 kg. These yields are higher than the 6,700 kg milk yield described in the highest emission factors for dairy cattle using the Tier I methodology: in North America, high yield cows fed on grains and high quality forage have EF of 118 kg CH₄ cow^-1 yr^-1(22.

Enteric fermentation EF values differed between primiparous and multiparous cows when the adjustment was made at 305 d (P=0.0025) to 365 d (P=0.0011), although emission intensity remained between 17.6 and 18.6 g CH₄ kg milk (Table 3). For primiparous and multiparous cows, methane EF differed depending on the milk and fat production curve modeled for each herd and lactation number (Figure 1). Milk yield was the most valuable factor for describing methane emissions in the analyzed database, although its representation was simplified when expressed as emission intensity per production unit.

Table 3 Enteric methane emission factors (EF) and intensity (EI) per lactation based on duration adjustment from 365 or 305 days 

ab Different letters in the same column indicate difference (P<0.0025). SE= Standard error.

Figure 1 Emission of enteric methane in primiparous (A; n = 16) and multiparous cows (B; n = 76). Lines represent the mean and standard error 

The present study contributed to characterize the intensive system milk production in Mexico and reducing uncertainty in application of the IPCC Tier II methodology. However, methane conversion factor (Ym), feed digestibility and feed composition values were taken from the literature, following IPCC best practices for inventory generation^11,22. Future research will need to include the wide diversity of management practices and environmental conditions present in any given country. Validation of the enteric CH₄ emissions estimates used here could only be by using scaled data for CH₄ measurements in individual animals or by disaggregating atmospheric measurements to the farm level by applying gas measurement and turbulent flow systems. Indeed, the Tier II methodology contains greater estimation and systemic errors than do mechanistic models of enteric CH₄ emission^20,26. This methodology is intended for reporting regional and country inventories, and does not consider emissions associated with electric energy use, processing, fuels, etc. Grain import and transport of liquid milk are important aspects of the dairy market structure in Mexico. When applied in the future, the IPCC methodologies will need to be paired with direct measurements, comparisons with mechanistic models and calculation of lifecycle emissions for production systems.

The emission factors and intensities reported in the literature need to be compared considering environmental conditions, productive systems, livestock management practices, the methane production estimation model used, available databases and study depth, among other factors. The present emission intensity results (17.6 to 18.6 g CH₄ kg milk) were within previously reported ranges for high yield cows: 25 g CH₄ L^-1 milk (enteric fermentation and manure) for the dairy cow population in Mexico⁹; 15.4 to 17.6 g CH₄ kg^-1 milk (enteric fermentation) in response to changes in the proportions of grazing to corn silage in the Netherlands²⁶; and 18.0 to 19.1 CH₄ g kg^-1 in experiments using different proportions of alfalfa and corn silage²⁷.

The present EF per head results (201 to 223 kg CH₄ cow^-1 in 365 d) were higher than the average reported in some national inventories: 120.76 kg in the USA²⁸; 113 kg in the Netherlands²⁹; and 128.8 to 139 kg in Germany¹⁹. This discrepancy is probably due to the fact that the AHM database does not include all the cows in a herd, since it focuses on high-yield dairy cows and longer lactations (365 vs 305 d). Emission of GHG increases as animal yield rises, but emissions per kilogram of milk corrected for fat content decreases notably as production increases¹². In addition, Tier II uses a fixed Y_m, which can result in a 12.5 % overestimation in dairy cattle, negatively affecting its accuracy^19,20. Country-level solutions are needed in enteric methane inventories to derive feed composition and properties as a function of animal performance and thus generate more accurate calculations¹⁹.

Enteric methane emission estimates for cattle using emission factors specific to Mexico have been done. In these studies, enteric emissions are represented as a function of energy intake in the feed and typical body mass of cattle classified by age and climatic region, but milk production levels are not considered¹⁰. One national inventory applied a modified Tier I methodology using a 305 d lactation, a 60 d dry period, a 4,400 kg average yield per lactation in a typical 400 kg cow; and a lactation curve model³⁰. However, this lactation model was obtained from data fitted for yield from 100 Friesian cows in the third lactation from Reading, UK^31,32. Another national inventory used milk yield from FAO reports (2,800 to 4,541 kg milk cow^-1 yr^-1), milk composition and feed rations taken from previous studies, and assumed that all dairy cows were confined⁹.

The historical contribution of methane emissions is significant, but simplifications need to be made in the model so that available databases can be used; for example, averaging the dairy and non-dairy cattle emission factors to then multiply it by the cattle population of each country⁸. Published studies expose the limitations in the availability and quality of cattle activity and population data, highlighting the need for further research in an effort to adequately quantify how much of GHG emissions originate in the livestock sector^3,4,5 . Another challenge is difficulty in data analysis caused by inconsistency among the units used to report GHG emission figures, particularly when SI (Sistème international d’unités) units are not employed. For instance, in one study the units Gg (10⁹ g) and Mt (an abbreviation interpreted as megaton - equivalent to 10¹² g, or as metric ton equivalent to 10⁶ g) are used interchangeably and produce the same methane emission result for dairy cattle in Mexico⁹. In another example, worldwide methane concentration from enteric fermentation from 1890 -1998 was calculated as 1,183,095 Gg CH₄⁸, but the same authors later published the figure of 8 million Gg Ceq³³ for the same concept. These two figures are not comparable using acceptable values for conversion of methane into carbon dioxide equivalents². Finally, before any attempt at scaling is made to calculate GHG inventories, comparisons will be needed between the emission factors or intensities of different animal classes and productive models.

Conclusions and implications

Analysis with the incomplete gamma formula and lactation curves showed high-yield dairy cows in Querétaro state to differ between primiparous and multiparous individuals for the studied variables. Average lactation was one-year (365 d) long. Enteric methane emission during lactation differed between primiparous and multiparous cows, whether at 365 d or adjusted to 305 d. However, lactation number or duration had no effect on average enteric methane emission intensity per kilogram of milk produced. The uncertainty present in enteric methane emissions inventories can be reduced by incorporating data for milk yield, lactation duration, and herd composition in terms of primiparous and multiparous cows. Milk yield can be represented using the incomplete gamma formula. However, an acceptable simplification would be to multiply average daily milk yield by the number of lactation days. This was shown to be viable in the present study in which relatively flat, highly precise curves were identified for high-yield cows. If enteric methane emission on a dairy farm is assumed to be 75 % of its GHG emissions, and the emission factors are accounted for, then a dairy producer managing more than 2,500 cows would surpass the 25,000 Gg Ceq threshold established in the Rules of the General Climate Change Law for the National Register of Emissions for the Republic of Mexico (Reglamento de la Ley General de Cambio Climático en Materia del Registro Nacional de Emisiones para la República Mexicana - RENE), and the farm would be forced to report its emissions beginning in 2018. In response to applicable laws (e.g. NMX-SAA-14064/1-IMNC- 2007, NMX-SAA-14064/2-IMNC-2007, among others), the Mexican Certification Organization (Entidad Mexicana de Acreditación) maintains current lists of the organizations that implement GHG emissions verification and validation to meet RENE stipulations.