Análisis de curvas de lactancia en vacas Holstein de la cuenca lechera de Chiriquí, República de Panamá

R Guerra Montenegro, A Hernández Rodríguez¹ y A Menéndez Buxadera²

Universidad de Panamá, Facultad de Ciencias Agropecuarias, Carretera Panamericana, Corregimiento de Chiriquí, Chiriquí, Panamá.
rguerram09@gmail.com
¹ Instituto de Ciencia Animal. Carretera Central, km. 47 ½. San José de las Lajas. Mayabeque. Cuba. C P 32 700.
² Consultor independiente

Resumen

Con el objeto de determinar la producción de leche y estimar la curva de lactancia en vacas Holstein de la cuenca lechera de Chiriquí, en la República de Panamá, se utilizaron 90 342 registros mensuales de la producción de leche en el día de control (PDC), de dos fincas controladas durante los años 1997 a 2016. Dicha información incluyó 9 742 lactancias pertenecientes a 3 196 vacas. Se utilizó el procedimiento GLIMMIX del SAS para obtener las producciones medias por orden de pesaje ajustadas a los factores incluidos en el modelo. Las mismas se usaron para estimar los parámetros de la curva de lactancia (β ₀, β₁ y β₂)), producción al pico (PP) y tiempo de pico (TP) por medio de la regresión no lineal, con el método modificado de Gauss-Newton, del procedimiento NLIN del SAS. Se evaluaron las funciones lineales Cuadrático logarítmica (CL) y Polinomial inversa (PInv) y no lineales Gamma incompleta (GI), Wilmink (W) y Cobby-Le Du (CLD) para determinar la mejor curva de ajuste. La producción media se inició con 24.47 kg/día, el pico de producción se alcanzó entre el segundo y tercer pesaje y el valor mínimo (16.37 kg/día) en el décimo pesaje que coincidió con la terminación de la lactancia. Las desviaciones, entre la producción de leche total estimada (PTE) y la observada fueron bajas, para todas las funciones, y la producción de leche total estimada por la curva de lactancia (6852 kg) fue próxima al valor observado (6425 kg). El PP estimado para las diferentes funciones se mantuvo entre los 28.1 y los 29.6 kg/día y el TP varió entre los 16 y 41 días de lactancia. Los coeficientes de determinación ajustados fueron mayores en las funciones W, GI y CL (99.99, 99.98 y 99.87%) y el menor valor lo obtuvo la función PInv (94.39%). En general, las funciones GI, W y CL fueron las de mejor ajuste ya que presentaron un menor error estándar de la estimación (0.002, 0.001 y 0.002) , menor suma de cuadrado de la estimación (0.16, 0.09 y 0.13), además de presentar una menor desviación de los residuos y valores de Durbin Watson cercanos a 2; por lo que se puede concluir que las funciones GI, W y CL mostraron un mejor ajuste a los criterios de selección utilizados y, por tanto, predijeron mejor la producción de leche en las condiciones de la cuenca lechera de Panamá.

Palabras clave: funciones lineales, funciones no lineales, Holstein, parámetros

Analysis of lactation curves of Holstein cows in the milk basin region of Chiriquí, Republic of Panama

Abstract

To determine the milk production and estimate the lactation curve in Holstein cows of the Chiriquí dairy basin in Panama republic, 90 342 monthly records of milk production were used on the control day (PDC), from two farms during the years 1997 to 2016. Information included 9 742 lactations belonging to 3 196 cows. The GLIMMIX procedure of SAS statistical software was used to obtain the average productions of PDC adjusted to the factors included in the model. They were used to estimate the parameters of the lactation curve (β₀, β ₁ and β₂)), peak production (PP) and peak time (TP) by means of the non-linear regression, with the modified Gauss-Newton method, of the SAS NLIN procedure. The linear functions Quadratic logarithmic (CL), Inverse polynomial (PInv) and nonlinear Gamma incomplete (GI), Wilmink (W) and Cobby-Le Du (CLD) were evaluated to determine the best model. The started average production was 24.47 kg/day, the production peak was reached between the second and third milk weighing and the minimum value (16.37 kg/day) was in the tenth weighing that coincided with the end of lactation. The deviations, between the total estimated production (PTE) and that observed were low, for all functions, and the total milk production estimated by the lactation curve (6852 kg), was close to the observed value (6425 kg). The milk production peak estimated for the different functions was maintained between 28.1 and 29.6 kg/day and peak time varied between 16 and 41 days of lactation. The adjusted coefficients of determination were greater for the W, GI and CL functions (99.99, 99.98 and 99.87%) and the lowest value was for the PInv function (94.39%). Functions CLD, W and CL were the best fit, since they presented a lower standard error of the estimate (0.002, 0.001 and 0.002) and a lower sum of square estimation (0.16, 0.09 and 0.13). In addition to presenting a lower deviation of residual and values of Durbin Watson close to 2, we can conclude that the GI, W and CL functions show the best fit to the selection criteria used, so are those that best predict milk production under the conditions of the dairy basin of Panama.

Key words: Holstein, linear functions, milk yield, non-linear functions, parameters

Introducción

La lechería en Panamá se ha considerado por años como una actividad de subsistencia, donde los productores, en su mayoría pequeños, trabajan con los requerimientos básicos de producción, sin optimizar sus recursos y hacer más eficiente dichos procesos. Estudios realizados por el Ministerio de Desarrollo Agropecuario en el 2011, indican que existe una gran población productora de leche. En 1985 existían en el país alrededor de 4 503 productores de leche en todas las categorías (A, B, C), de los cuales solo unos 250 producían leche grado A (MIDA 2011). Esta cantidad ha ido en aumento a través de los años.

Actualmente el sector lechero presenta grandes deficiencias en cuanto a su producción, ya que el mismo no cubre la demanda nacional. En Panamá se producen 195 millones de litros de leche y el mercado nacional consume 300 millones de litros anualmente (Guevara 2012), por lo que la importación se hace necesaria, a la vez que ofrece un nicho de desarrollo importante para aumentar la productividad y cerrar esta brecha.

De acuerdo a Silvestre et al (2006), el estudio de las propiedades matemáticas de la curva de lactancia proporciona información resumida acerca de la producción de ganado lechero, que es útil en la toma de decisiones de gestión, por ejemplo, la vigilancia de la salud y la alimentación que permiten optimizar la productividad de los hatos lecheros.

El único estudio encontrado de curvas de lactancia en Panamá fue reportado por Villalobos et al (2016) en vacas cruzadas Pardo Suizo, en la región de Azuero. Sin embargo, en la cuenca lechera de Chiriquí, que es la región más productora del país, no se dispone de estudios que muestren el comportamiento lechero siendo, por ello, el objetivo del presente trabajo determinar la producción de leche y estimar la curva de lactancia en vacas Holstein de la cuenca lechera de Chiriquí, Panamá.

Materiales y métodos

Se utilizaron los datos de producción láctea de dos fincas ubicadas en la cuenca lechera de Chiriquí, en la república de Panamá. En dicha región existen dos estaciones claramente definidas, la de lluvias (mayo a diciembre), y la de seca (enero a abril). La precipitación anual promedio es 3000 mm y la temperatura media anual 26 °C (FAO 2017).

Los animales estuvieron en pastoreo a voluntad, a base de Cynodon nlemfuensis, Echinochloa polystachia, Brachiaria brizantha, Pennisetum clandestinum, en praderas con periodos de ocupación de 1 día por 21 días de descanso. En épocas secas o de mucha lluvia se suplementó con pasto de corte de Pennisetum purpureum y heno de Digitaria swazilandensis el cual se ofreció durante su estancia en la sala de espera para el ordeño.

A las vacas en producción se les suministró además de pasto, ensilaje o heno; y pienso concentrado en dependencia del nivel productivo de la vaca, a razón de 0.5 kg por cada litro de leche producido. En las vacas de alta producción (más de 35 kg/día), se le adicionó a la dieta grasas sobrepasantes, proteínas sobrepasantes y sustancias bufferantes.

Las vacas se ordeñaron dos veces al día, y el lote de vacas en ordeño se dividió en tres sub lotes llamados de alta, media y baja producción, los cuales se determinaron con base en los pesajes de leche y el historial de parto, siendo el factor determinante la cantidad de leche producida según el pesaje semanal (mañana y tarde). Con base en estos grupos se estableció el tiempo de ordeño de cada vaca.

Para el estudio se analizaron 90 342 registros de la producción de leche en el día de control (PDC), de 9 742 lactancias (1 a 6), agrupándose las lactancias superiores a la sexta en esta última categoría, controladas mensualmente, pertenecientes a 3 196 vacas Holstein, en los años 1997 a 2016.

Se aplicó un proceso de depuración donde se eliminaron los registros fuera del rango de la media ± 3 desviaciones estándar, lactancias con menos de 4 pesajes mensuales controlados, lactancias que iniciaron el pesaje con más de 35 días después del parto, y aquellas que truncaron las lactancias en los 305 días de duración.

Para la depuración y determinación de estadígrafos descriptivos (media, desviación estándar y coeficiente de variación), se utilizó el paquete estadístico SAS, versión 9.3 (2010).

Se aplicó un modelo lineal generalizado mixto (Wolfinger y O’Connell 1993), con ayuda del procedimiento GLIMMIX del SAS v.9.3 (2010). Se utilizaron las funciones de enlace en estos modelos según la distribución de las variables, para obtener los valores predichos de las medias, a partir de la estimación del vector β, según la función inversa a la función de enlace.

El modelo utilizado fue el siguiente:

Y_gijklm = V_g +NL_i + OP_j + APA_k + MPA_l + F_m + (NLOP)_ij + (APAMPA)_kl + (FAPA)_mk + (FMPA)_ml + (FNL)_mi + (FOP)_mj + e _gijklm

Dónde:

Y_gijklm= PDC

V_g = efecto aleatorio de la g-ésima vaca (g=1…3 196)

NL_i = efecto fijo del i-ésimo número de lactancia (i=1…6)

OP_j = Efecto del j-ésimo orden de pesaje (j=1…10)

APA_k = efecto fijo del k-ésimo año de pesaje (k=1997…2016)

MPA_l = efecto fijo del l-ésimo mes de pesaje (l=1…12)

F_m = efecto fijo de la m-ésima finca (m=1,2)

(NLOP)_ij = interacción entre i-ésimo numero de lactancia y j-ésimo orden de pesaje

(APAMPA)_kl = la interacción entre k-ésimo año de pesaje y el l-ésima mes de pesaje

(FAPA)_mk = interacción entre m-ésima finca y k-ésimo año de pesaje

(FMPA)_ml = interacción entre m-ésima finca y l-ésimo mes de pesaje

(FNL)_mi = interacción entre m-ésima finca e i-ésimo número de lactancia

(FOP)_mj = interacción entre m-ésima finca y j-ésimo orden de pesaje

e_gijklm = error aleatorio debido a cada observación NID~ (0, s2e).

Se aplicó la prueba de Tukey-Kramer para la comparación múltiple de las medias de mínimos cuadrados (Kramer 1956).

Para determinar la curva de lactancia media se creó una nueva base de datos, con ayuda del paquete Microsoft Excel del Office (2007), con las producciones medias por orden de pesaje, que se ajustaron a los factores no genéticos y genéticos, a partir de los análisis del procedimiento GLIMMIX del SAS v.9.3 (2010).

Se evaluaron las funciones lineales Cuadrático logarítmica (CL) y Polinomial inversa (PInv) y no lineales Gamma incompleta (GI), Wilmink (W) y Cobby-Le Du (CLD) para determinar la del mejor ajuste. Las ecuaciones son las siguientes:

Cuadrático Logarítmica (Bianchini-Sobrinho 1984):

 Polinomial Inversa (Nelder 1966): 

Gamma Incompleta (Wood 1967):

Wilmink (1987):

Cobby y Le Du (1978):

En las funciones: Y (t) es la producción de leche diaria, t es tiempo en días de pesaje,β₀, β₁ y β₂, son los parámetros específicos de cada función donde:β₀ = producción inicial,β₁ = producción al pico, y β₂= descenso de la producción post pico hasta el secado.

La estimación de los parámetros de la curva de lactancia (β₀, β₁ y β₂), producción al pico (PP) y tiempo de pico (TP) se obtuvo por medio de la regresión no lineal, con el método modificado de Gauss-Newton, del procedimiento NLIN del SAS v. 9.3 (2010).

Los indicadores de PP y TP se estimaron de acuerdo a cada modelo ajustado, así en la función GI se definió el PP y TP de acuerdo con Wood (1967) como:

El punto de máximo se determinó para la CL (Bianchini-Sobrinho 1984), por:

El punto de máximo se definió, para la W (Wilmink 1987), como:

El TP y PP se determinaron en la PInv (Nelder 1966), por:

El TP Y PP, para la función CLD (Cobby y Le Du 1978), como:

Los criterios que se utilizaron para la selección de la mejor función, fueron los propuestos por Guerra et al (2003):

• Coeficiente de determinación (R²), estimado para las curvas medias.
• Coeficiente de determinación ajustado a los grados de libertad (R ²A), estimado para las curvas medias.
• Error estándar de estimación (EEE).
• Suma de cuadrados de la estimación (SCE).
• La prueba de significación del modelo y de los parámetros.
• Estadístico de Durbin Watson (DW) (Durbin y Watson 1971).
• Distribución de los residuos.

Resultados y discusión

Las constantes mínimas cuadráticas de las PDC ajustadas por orden de pesaje, de la población en estudio (Figura 1) variaron de acuerdo con la fase de lactancia. La producción media inició con 24.46 kg/día, el pico de producción se alcanzó entre el segundo y tercer pesaje (correspondiente con los días 30 y 60 de realizado el pesaje) y el valor mínimo (16.37 kg/día) en el décimo pesaje que coincidió con la terminación de la lactancia. Una característica interesante de la curva es que después del PDC₃, se observó una tendencia de declinación continua en la producción hasta el PDC₁₀.

Estos resultados se corresponden con los reportados en vacas Holstein por Vargas et al (2016) en Colombia con animales de uno a tres partos quienes observaron el pico de lactancia entre los 41 y 45 días; Mesén (1999) en las condiciones de Costa Rica con vacas Holstein de primera y segunda lactancia reportaron el pico de producción entre las 56 y 63 días, y Dorneles et al (2009) en Brasil quienes informaron el pico de lactancia en el segundo pesaje correspondiente a los 55 días de lactancia.

Figura 1. Medias de PDC ajustadas por orden de pesaje.

Las desviaciones, entre la producción total estimada (PTE) y la observada fueron bajas para todas las funciones (Tabla 1). La PP estimada para las diferentes funciones se mantuvo entre los 28.07 y los 29.64 kg/día lo cual coincide con los resultados (27 y 35 kg/día) de múltiples autores con vacas Holstein en condiciones de manejo y alimentación similares (Dorneles 2011; Cañas et al 2012; Cuatrín et al 2016).

El TP para las funciones CL, GI y W estuvieron dentro del rango de valores reportados para la raza Holstein por Fernández et al (2011) en Cuba y Cuatrín et al (2016) en Argentina quienes estimaron el TP a los 35 días; así como los resultados de Cañas et al (2012) en Colombia con un TP de 44 días; en tanto, Bouallegue et al (2013) en Túnez observaron el pico de lactancia entre los 42 y 56 días. Las funciones CLD y PInv mostraron un TP de 16 y 23 días respectivamente. Este resultado no corresponde a lo reportado por Peña et al (2005) como normal para las vacas lecheras especializadas.

Al considerar los parámetros estimados para la curva de lactancia media (Tabla 2), se encontraron diferencias (p˂0.0001) para todas las funciones analizadas. En cuanto a los criterios de selección analizados, las funciones presentaron ajustes adecuados a los mismos. Los coeficientes de determinación ajustados fueron mayores en las funciones W, GI y CL (99.99, 99.98 y 99.87%) y el menor valor lo obtuvo la función PInv (94.39%). En general, las funciones GI, W y CL fueron la de mejor ajuste ya que presentaron un menor error estándar de la estimación (0.002, 0.001 y 0.002) y menor suma de cuadrado del error (0.16, 0.09 y 0.13).

La prueba de Durbin-Watson para las funciones GI, W, CLD Y CL indicó valores próximos a dos que son los deseables lo cual significó ausencia de correlación (Hoffman y Vieira, 1987) no así la función PInv que presentó un valor inferior a 1.

Finalmente se consideró el análisis de los residuos (Figura 2), en conjunto con la prueba DW=2(1-r) siendo r el valor de las auto-correlaciones, para elegir las funciones que mejor ajuste, observándose que los residuos fueron variables en todas las funciones y se distribuyeron por encima y debajo del eje de la ordenada. Es importante considerar que los residuos consideran los efectos genéticos, pues son desviaciones de una misma vaca y los efectos del ambiente permanente.

Tomando en consideración los criterios expuestos anteriormente, las funciones GI, W y CL fueron las que mostraron menor desviación de los residuos y valores de DW cercanos a 2, lo cual coincidió con los altos valores de R²A, bajos EEE y SCE y una estimación de PP y TP similar a los reportados como normales para la raza, lo cual indicó que son las funciones que mejor predicen la producción de leche en vacas Holstein, en la cuenca lechera de Chiriquí en Panamá.

Estudios con vacas Holstein en Irán (Elahi-Torshizi et al 2011), donde utilizaron 7 modelos matemáticos para modelar la curva de lactancia de vacas Holstein, mediante criterios de seleccion R², R ²A, CME y DW, en Colombia (Cañas et al 2012), evaluaron los modelos de Wood, Brody, Wilmink y Papajcsik- Bordero, utilizando como criterio de selección el valor de R²A y en Brasil (Ferreira et al 2015) quienes evaluaron 4 modelos matemáticos para curvas de lactancia (Wood, Brody, Djikstra y Pollott) coincidieron en elegir como el modelo de mejor ajuste el de Wood. Diferentes resultados obtuvieron Strucken et al (2011) quienes reportaron el modelo Wilmink como el de mejor ajuste; por lo que se considera que la elección del mejor modelo varía entre autores, ya que depende de la muestra utilizada en el estudio y de las condiciones de manejo y alimentación a la estuvieron sujetos los animales en estudio

Figura 2. Análisis de los residuos para las funciones CLD, CL, PInv, W y GI.

Conclusiones

• Las funciones GI, W y CL predijeron mejor la producción de leche en vacas Holstein, bajo las condiciones de la cuenca lechera de Panamá, ya que presentaron los mejores ajustes a los criterios de selección considerados.

Agradecimientos

Esta investigación es resultado de la beca doctoralBIDP 2014-011 del programa EXCELENCIA PROFESIONAL – DOCTORADO EN AREAS ESPECIFICAS DEL CONOCIMIENTO, de la Secretaria Nacional de Ciencia, Tecnología e Innovación (SENACYT) en conjunto con el Instituto para la Formación y Aprovechamiento de los Recursos Humanos (IFARHU) de la República de Panamá con apoyo de la Universidad de Panamá – Facultad de Ciencias Agropecuarias. De igual manera a los productores que facilitaron los datos de las fincas para los análisis realizados.

Referencias

Bianchini Sobrinho E 1984 Estudo da curva de lactação de vacas da raçã Gir. Tesis de doctorado, USP, Ribeirão Preto, p. 84.

Bouallegue M, Hadad B, Aschi M and Ben Hamoud M 2013 Effect of environmental factor of lactation curves of milk production traits in Holstein-Friesian cows reared under North African conditions. Livestock Research for Rural Development 25(5): http://www.lrrd.org/lrrd25/5/boua25075.htm

Cañas J, Ceron-Muñoz M and Corrales J 2012 Modelación y parámetros genéticos de curvas de lactancia en bovinos Holstein en Colombia. Revista MVZ Cordoba 17 (2): 2998 - 3003

Cobby J and L Ledu 1978 On fitting curves to lactation data. Animal Production 26:127-134.

Dorneles C, Cobuci J, Rorato P, Weber T, Lopes J and Olivera H 2009 Estimacao de parámetros genéticos para producao de leite de vacas da raca Holando via regressao aleatória. Arquivo Brasileiro de Medicina Veterinária e Zootecnia 61: 407-412

Durbin J and Watson G 1971 Testing for serial correlation in least squares regression. III. Biometrika. 58: 19.

Elahi Torshizi M, Aslamenejad A, Nassiri M and Farhangfar H 2011 Comparison and evaluation of mathematical lactation curve functions of Iranian primiparous Holstein. South African Journal of Animal Science 41 (2):104 – 115

FAO 2017 Informe sobre la situación de los recursos zoogenéticos en Panamá. Documento electrónico, 57p: //ftp.fao.org/docrep/fao/011/a1250f/annexes/CountryReports/Panama.pdf

Ferreira A, Henrique D, Vieira R, Maeda E and Valotto A 2015 Fitting mathematical models to lactation curves from Holstein cows in the southwestern region of the state of Parana, Brasil. Anais da Academia Brasileira de Ciencias 87 (1).

Guerra W, Cabrera A and Fernández L 2003 Criteria for the selection of statistical models in scientific research. Cuban Journal of Agricultural Sciences 3: 37.

Guevara T 2012 Panorama del sector lechero nacional y avances de la cadena agroalimentaria de la leche. Panamá, 44p.

Hoffman R and Vieira S 1987 Análise de regressão. Uma introdução a econometria. São Paulo: Hucitec, 379p.

Kramer C Y 1956 Extension of multiple range tests to group means with unequal numbers of replications. Biometrics, 12, 307 - 310.

Mesen B 1999 Evaluación de la curva de lactancia de la raza Holstein. Tesis, Facultad de Agronomía, Universidad de Costa Rica, 80 p.

MIDA 2011 Plan de Acción Estratégico del Sector Agropecuario 2010-2014. http://www.bda.gob.pa/Templates/transparencia/PAE-MIDA.pdf

Nelder J A 1966 Inverse polynomials, a useful group of multifactor response functions. Biometrics 22:128-135.

SAS 2010 SAS User’s guide: Statistics. Version 9.3. De. SAS Institute. INC, Cary, N.C., USA.

Strucken E, de Koning D, Rahmatalla S and Brockmann G 2011 Lactation curve models for estimating gene effects over a timeline. Journal of Dairy Science 94:442–449.

Vargas O, Ramírez S, Gaitán J and Corrales J 2016 Modelación de curvas de lactancia para producción de leche por parto en vacas Holstein en Boyacá, Colombia. Revista Colombiana de Zootecnia 4: 13 – 20.

Villalobos A, Guerrero B, Hassán J and Herrera D 2016 Curvas de lactación de bovinos mestizos pardo suizo en la región de Azuero. Ciencia Agropecuaria 24: 103 – 110.

Wilmink JB 1987 Adjustment of lactation yield for age at calving in relation to level of production. Livestock Production Science 16:335–345.

Wolfinger R and O’Connell M 1993 Generalized linear models: a pseudo likelihood approach. Journal of Statistical Computation and Simulation. 48: 233-243.

Wood P 1967 Algebraic model of the lactation curve in cattle. Nature 216:164-165.

Received 16 December 2017; Accepted 24 January 2018; Published 1 April 2018

Go to top