Volumen 30(3):733-750. Septiembre-diciembre, 2019

e-ISSN 2215-3608, doi:10.15517/am.v30i3.35984

http://www.revistas.ucr.ac.cr/index.php/agromeso

Efecto del estrés calórico sobre la producción del ganado lechero en Costa Rica1

Heat stress effect on dairy cattle production in Costa Rica

Juan Ismael Ruiz-Jaramillo2, Bernardo Vargas-Leitón2, Sergio Abarca-Monge3, Hugo G. Hidalgo4

1 Recepción: 13 de febrero, 2019. Aceptación: 10 de abril, 2019. Esta Investigación formó parte de tesis de Maestría en Producción Animal Sostenible del primer autor. Posgrado Regional en Ciencias Veterinarias Tropicales, Universidad Nacional, Heredia, Costa Rica.

2 Universidad Nacional de Costa Rica, Posgrado Regional en Ciencias Veterinarias Tropicales, Campus Benjamín Núñez, Lagunilla, Heredia, Costa Rica. jirjjaramillo@gmail.com, bernardo.vargas.leiton@una.cr (https://orcid.org/0000-0002-1778-9672).

3 Instituto Nacional de Innovación y Transferencia en Tecnología Agropecuaria (INTA), San José, Costa Rica. sabarca@inta.go.cr

4 Universidad de Costa Rica, Centro de Investigaciones Geofísicas y Escuela de Física, San José, Costa Rica. hugo.hidalgo@ucr.ac.cr

Resumen

Introducción. La producción lechera en zonas tropicales se caracteriza por una alta dependencia del recurso forrajero, lo que la hace sensible a las variables climáticas, de las cuales existe información limitada sobre su efecto en el rendimiento de razas lecheras especializadas. Objetivo. El objetivo del presente estudio fue evaluar la asociación entre el índice de temperatura y humedad (ITH) y la producción diaria de leche (kgl) de vacas de razas Holstein, Jersey y sus cruces. Materiales y métodos. El estudio se efectuó en las regiones Norte y Central de Costa Rica, con datos recolectados entre 1990 y 2015. Los valores de ITH se obtuvieron mediante modelos predictivos lineales y modelos espaciales de autocorrelación, ajustados sobre 3547 registros mensuales de temperatura máxima y humedad relativa de diecisiete estaciones meteorológicas. Se analizaron 6 478 582 registros de kgl provenientes de 418 hatos lecheros, mediante un modelo lineal mixto generalizado (GLMM), que consideró los efectos fijos de la raza, el año y el mes calendario, número de parto, etapa de lactancia, raza×etapa de lactancia, ENOS (Efecto “El Niño”), ITH y raza×ITH, además del efecto aleatorio de la vaca. Resultados. Todos los factores presentaron efecto altamente significativo (p<0,0001) sobre la kgl. Para el rango de ITH entre 72-88 se estimaron reducciones lineales de 0,41 (r2=0,94), 0,36 (r2=0,95) y 0,29 (r2=0,82) kg día-1 para Holstein, Holstein×Jersey y Jersey, respectivamente. No se observaron tendencias significativas en los kgl para ITH <72. Las pérdidas económicas atribuibles a la mayor exposición al estrés en la región norte comparada con la región central se estimaron en $680, $587 y $477 por lactancia por vaca para Holstein, Holstein×Jersey y Jersey, respectivamente. Conclusión. Hubo relación inversa entre el índice de temperatura y humedad, y la producción de leche en vacas Holstein, Jersey y sus cruces, en esta región tropical.

Palabras clave: producción lechera, tolerancia al calor, cambio climático, estrés térmico.

Abstract

Introduction. Dairy production in tropical zones is characterized by a high dependence on the forage resource, which makes it sensitive to climatic variables of which there is limited information on their effect on the performance of specialized dairy breeds. Objective. The objective of this study was to evaluate the association between temperature and humidity index (THI) and daily milk production (kgl) of cows breeds Holstein, Jersey, and their crosses. Materials and methods. The study was conducted in the northern and central regions of Costa Rica, with data collected between 1990 and 2015. THI values were obtained using linear predictive models and spatial autocorrelation models, applied to 3,547 monthly records of maximum temperature and relative humidity from seventeen weather stations. 6,478,582 kgl records from 418 dairy herds were analyzed, using a generalized mixed linear model (GLMM), which considered the fixed effects of breed, year and month, birth number, lactation stage, breed × lactation stage, ENSO (Effect “El Niño”), ITH and breed×ITH, in addition to a random effect of the cow. Results. All factors had a highly significant effect (p<0.0001) on kgl. For the ITH range between 72-88, linear reductions of 0.41 (r2=0.94), 0.36 (r2=0.95) and 0.29 (r2=0.82) kg day-1 were estimated for Holstein, Holstein×Jersey and Jersey, respectively. No significant trends were observed for kgl when ITH <72. The economic losses attributable to the greater stress exposure in the north compared to the central region were estimated at $680, $587 and $477 per lactation and cow for Holstein, Holstein×Jersey and Jersey, respectively. Conclusion. There was an inverse relationship between temperature and humidity index and milk production in Holstein, Jersey and crossbreed cows for this tropical region.

Keywords: milk production, heat tolerance, climate change, heat stress.

Introducción

El estrés calórico se define como la sensación de malestar que experimentan los organismos cuando los mecanismos encargados de mantener la temperatura interna deben realizar esfuerzos ante la permanencia en un ambiente determinado (West et al., 2003). Una forma indirecta de cuantificar el estrés calórico al que se expone un organismo es mediante el índice de temperatura y humedad (ITH), que representa el efecto combinado de la temperatura y la humedad ambiental asociados al nivel de estrés térmico (Bohmanova et al., 2007).

El ITH fue propuesto originalmente para cuantificar el grado de estrés en humanos, pero ha sido modificado y adaptado para su uso en otras especies, entre ellas el ganado bovino (NRC, 1971; Bohmanova et al., 2007; Hammami et al., 2013). Actualmente, se utilizan distintas fórmulas de cálculo de ITH, que asignan diferentes pesos relativos a la temperatura y humedad, o que incluyen variables adicionales en el cálculo, tales como la radiación solar (Bohmanova et al., 2007; Zimbelman et al., 2009; Hammami et al., 2013).

En guías de manejo para hatos lecheros utilizadas en Estados Unidos, se ha indicado que la zona de confort para el bovino se encuentra en valores de ITH por debajo de 72 (Armstrong, 1994), lo que corresponde a una temperatura de 23 °C con una humedad relativa del 80 %. Valores de ITH entre 72 y 79 se asocian con estrés moderado, entre 80 y 89 son causantes de estrés severo y valores por encima de 89 ocasionan estrés grave (Armstrong, 1994). Sin embargo, se ha observado que estos rangos pueden variar en función de factores geográficos, de manejo o también factores biológicos como la raza, el nivel de actividad o el nivel de producción (Dunn et al., 2014).

Existe evidencia de que las vacas lecheras de alto valor genético, seleccionadas por su alta producción de leche, tienden a presentar mayor susceptibilidad al estrés calórico (Ravagnolo y Misztal, 2000; Hammami et al., 2013; Dunn et al., 2014; González, 2014). En las vacas Holstein el estrés térmico puede presentarse en niveles de ITH tan bajos como 65-69 (Bouraoui et al., 2002; Bryant et al., 2007; Zimbelman et al., 2009).

El estrés calórico se produce cuando el ITH es elevado y afecta negativamente la producción de leche (Bouraoui et al., 2002; West et al., 2003; Bohmanova et al., 2007; Gantner et al., 2011; Hammami et al., 2013; Zewdu et al., 2014; Garcia et al., 2015; Pragna et al., 2017). En el sur de Estados Unidos, la producción diaria de leche de vacas Holstein disminuyó cuando el ITH superó el valor de 72, con una tasa de reducción de 0,39 kg por cada unidad de incremento en ITH (Bohmanova et al., 2007). Otro estudio bajo condiciones similares informó sobre una reducción de 0,20 kg (Ravagnolo y Misztal, 2000). Mayores reducciones, en el orden de 0,88 y 0,60 kg, fueron reportadas en la producción diaria de vacas Holstein y Jersey, respectivamente, por cada unidad de incremento en ITH (West et al., 2003). En Etiopía, se reportó una reducción de hasta 190 kg en producción de leche por lactancia de vacas cruzadas Holstein×Deoni, por cada incremento de una unidad en el ITH máximo (Zewdu et al., 2014).

Entre los eventos fisiológicos desencadenados por el estrés calórico se mencionan el aumento de la temperatura rectal, la frecuencia respiratoria y el jadeo, que se manifiestan para buscar mantener la temperatura corporal (Pragna et al., 2017). Estos cambios conllevan una alteración en el patrón de alimentación y función del rumen, con reducción de la ingesta de materia seca y por consiguiente, de la productividad. El motivo de la reducción en la producción de leche es el balance energético negativo, ya que el animal intenta mantener la homeostasis para evitar la hipertermia. El estrés por calor también puede tener impacto negativo sobre la salud de la ubre, lo que en última instancia también conduce a la disminución de la leche (Pragna et al., 2017).

En conjunto, las pérdidas asociadas al efecto del estrés calórico sobre la producción lechera son cuantiosas. En algunos países de la franja tropical, se estima que las elevaciones en ITH causan pérdidas de 1,8 millones de toneladas de leche de la producción anual, lo cual representa más de 650 millones de dólares (Zewdu et al., 2014). En Estados Unidos, se calculó una pérdida anual de 897 millones de dólares por el mismo concepto (St-Pierre et al., 2003).

En Costa Rica, la ganadería lechera especializada se desarrolla principalmente en las regiones climáticas central y norte, donde los sistemas de producción son mayoritariamente basados en pastoreo con suplementación (72,7 %), y en menor escala semiestabulados (26,4 %) y estabulados (<1 %) (Vargas-Leitón et al., 2013). El clima de las regiones Central y Norte difiere en forma marcada; en la primera es seco en la parte baja, con temperatura promedio alrededor de 22 °C, templado en las partes medias y altas, con temperatura promedio alrededor de 15 °C (IMN, 2008). La precipitación promedio en esta región es de 2300 mm anuales en el valle occidental y 1700 mm en el valle oriental (IMN, 2008), con una marcada época seca entre los meses de diciembre hasta abril. La Región Norte presenta un clima tropical húmedo, con temperatura promedio de 20 °C en las partes altas y 26 °C en las partes bajas, y una precipitación de 3200 mm anuales (IMN, 2008); esta región es en su mayor parte plana y presenta un régimen de precipitación lluvioso todo el año, influenciado por la vertiente del Atlántico, con una disminución relativa de las lluvias en los meses de febrero, marzo y abril (IMN, 2008).

En Costa Rica y en zonas tropicales en general, la producción lechera y agrícola es altamente sensible al clima, por lo tanto, cambios en las variables climáticas pueden tener consecuencias directas sobre la producción y la disponibilidad de alimentos (Ríos e Ibrahim, 2008; Ordaz et al., 2010; Bouroncle et al., 2015). Eventos cíclicos tales como El Niño-Oscilación Sur (ENOS) o eventos climáticos extremos (frentes fríos, sequías, tormentas o huracanes) también han demostrado tener un alto impacto en la precipitación y temperatura de distintas épocas y años, afectando la producción agrícola y ganadera (Ordaz et al., 2010; Vallejos et al., 2012; Bonilla, 2014; Bouroncle et al., 2015).

En América Latina, los esfuerzos por determinar el impacto de variables climáticas sobre la producción agrícola han sido enfocados principalmente hacia los cultivos (Ordaz et al., 2010; Bouroncle et al., 2015), mientras que los estudios dedicados a ganadería son limitados (Garcia et al., 2015). Las zonas tropicales se caracterizan por estar expuestas a una mayor radiación solar y humedad, por lo que es importante evaluar y monitorear el impacto de las variables climáticas sobre las distintas razas lecheras, preferiblemente por periodos de tiempo prolongados. El objetivo del presente estudio fue evaluar la asociación entre el índice de temperatura y humedad (ITH) y la producción diaria de leche (kgl) de vacas de razas Holstein, Jersey y sus cruces.

Materiales y métodos

Fuente de datos climáticos

El cálculo del índice de temperatura y humedad (ITH) se realizó con base en series mensuales de temperatura y humedad relativa correspondientes al periodo comprendido entre 1990 y 2015 (veintiséis años), provenientes de diecisiete estaciones meteorológicas (diez automáticas y siete mecánicas), administradas por el Instituto Meteorológico Nacional (IMN), distribuidas en un rango de elevación entre los 15 y los 1850 msnm (Figura 1). Entre las estaciones, ocho se ubicaron en la Región Central y cinco en la Región Norte. Se agregaron además cuatro estaciones pertenecientes a la Región Caribe Norte, con el fin de lograr una mejor cobertura alrededor de los sitios de los hatos lecheros disponibles en la Región Norte.

Figura 1. Ubicación geográfica de los hatos (*) y las estaciones meteorológicas (+) incluidos en el estudio para evaluar la asociación entre el índice de temperatura y humedad (ITH) y la producción diaria de leche(kgl). Costa Rica, periodo 1990-2015.

Figure 1. Geographical location of herds (*) and meteorological stations (+) included in the study to evaluate the association between temperature and humidity index (THI) and daily milk production (kgl). Costa Rica, period 1990-2015.

Se contó con un total de 3547 pares de registros mensuales de temperatura máxima y humedad relativa en las diecisiete estaciones. Ninguna estación contó con las series mensuales completas para el periodo bajo estudio. El promedio de disponibilidad de información por estación para el cálculo del ITH fue 66,9 % con respecto al máximo posible de mediciones, que fue de 312 (26 años×12 meses), varió entre un mínimo de 14,7 % y un máximo de 85,9 %.

Inicialmente, se realizó un análisis de estadística descriptiva de las variables climáticas disponibles, incluyendo el cálculo de medidas de tendencia central y dispersión, así como una exploración gráfica de la distribución de las variables y sus tendencias en tiempo y espacio. Se identificaron y eliminaron valores inusuales posiblemente causados por errores de medición, el criterio empleado fue >3,5 desviaciones estándares dentro de un mes y una región, partiendo del supuesto de una distribución aproximadamente normal de las variables analizadas. No se realizó ningún tipo de imputación de datos climáticos, solo se incluyeron los meses para los cuales se contó con valores reales de temperatura y humedad relativa. Posteriormente, se procedió al cálculo del índice de temperatura y humedad (ITH), para lo cual se utilizó la fórmula propuesta por el National Research Council (NRC, 1971), definida como ITH=(1,8×T+32)-[(0,55-0,0055×HR)×(1,8×T-26)], donde T corresponde a la temperatura ambiental bulbo seco (°C) y HR es el porcentaje de humedad relativa. Esta fórmula ha sido sugerida para ganado bajo pastoreo (NRC, 1971). Diferentes estudios han utilizado tanto la temperatura promedio (West et al., 2003; Bomahnova et al., 2007) como la máxima (St-Pierre et al., 2003; West et al., 2003; Bryant et al., 2007) en el cálculo de ITH. En el presente estudio se calcularon ambos valores de ITH (promedio y máximo), pero la evaluación de asociación entre el ITH y producción se realizó con base en la temperatura máxima, con el fin de obtener una estimación de los valores más altos de ITH a que pueden estar expuestos los animales durante el día bajo condiciones locales.

Cálculo del índice de temperatura y humedad (ITH) en los sitios de los hatos

Con base en los datos de las diecisiete estaciones disponibles, se realizó una estimación de la variable ITH para los sitios específicos de los hatos bajo estudio, a lo largo de todos los meses y años, para los cuales existió información productiva en cada hato. Esta estimación se realizó en dos pasos, con base en el método descrito por Kolovos (2010), que combina predictores obtenidos a partir de un modelo lineal con estimadores estocásticos obtenidos a partir de un modelo de autocorrelación espacial. En el primer paso se obtuvieron predictores lineales de ITH a partir de un modelo lineal mixto generalizado, ajustado mediante el procedimiento GLIMMIX (SAS Institute Inc., 2013). Este modelo se ajustó con base en los datos reales de la variable ITH de las diecisiete estaciones meteorológicas. El modelo analizó la variable ITH en función de los efectos fijos de elevación (msnm), región, año y mes calendario, efecto climático denominado El Niño-Oscilación del sur o ENOS5, las interacciones región×mes, región×año y región×ENOS; así como el efecto aleatorio de variación dentro y entre estaciones meteorológicas. Se agregaron además los efectos de latitud, longitud y latitud×longitud, con el fin de eliminar posibles tendencias de superficie (Kolovos, 2010). Mediante el procedimiento PLM (SAS Institute Inc., 2013), la ecuación de predicción se utilizó para obtener predictores lineales del ITH en los sitios específicos de los hatos, para todos los años y meses en los cuales estos contaron con información de producción registrada.

En un segundo paso, los residuales libres de tendencia del modelo anterior se modelaron con el procedimiento VARIOGRAM (SAS Institute Inc., 2013) para explorar el ajuste de distintos modelos de autocorrelación espacial, entre ellos el gaussiano, el esférico, el exponencial y el sinusoidal. La selección del modelo de mejor ajuste se realizó con base en el criterio de minimización de la suma de cuadrados del error (Kolovos, 2010). En este paso se evaluó además el supuesto de anisotropía (Kolovos, 2010), que consiste en comparar el comportamiento de la semivarianza en cuatro diferentes direcciones (0, 45, 90 y 135 grados). El modelo seleccionado se utilizó para obtener estimadores aleatorios insesgados de los residuales del ITH en los sitios específicos de los hatos, utilizando el procedimiento KRIGE2D (SAS Institute Inc, 2013). En este procedimiento se utilizó el método de interpolación Kriging de tipo ordinario (Chilès y Desassis, 2018), que asigna ponderaciones basadas no solo en la distancia entre los puntos (estaciones meteorológicas) circundantes, sino también en la autocorrelación entre estos. La interpolación se realizó utilizando los datos reales de todas las estaciones circundantes con datos disponibles del ITH para el año y mes respectivo, ubicadas en un radio máximo de búsqueda de 0,30 grados decimales, equivalente a una distancia aproximada de 33,3 km con respecto al sitio de cada hato. En caso necesario, este radio se amplió hasta lograr ubicar un mínimo de tres estaciones con datos disponibles.

El predictor lineal del ITH obtenido (paso 1), se sumó con su correspondiente estimador aleatorio insesgado (paso 2) para obtener el valor predicho final del ITH en todos los sitios y meses de cada hato.

Fuente de información de producción

La información de la variable de producción diaria de leche (kgl) se obtuvo de la base de datos administrada por el proyecto CRIPAS (Romero et al., 2011) de la Universidad Nacional, Costa Rica. Se seleccionaron hatos de lechería especializada, con localización exacta conocida, ubicados en las regiones Central y Norte de Costa Rica, y que contaban con un mínimo de 1000 registros de producción de leche diaria (kgl) correspondientes a los primeros 360 días de lactancia de vacas de raza Holstein, Jersey o sus cruces, para el periodo comprendido entre el 01/01/1990 y el 31/12/2015. Estos criterios resultaron en un total de 418 hatos (171 en Región Norte y 247 en Región Central), ubicados en las provincias de Alajuela (n=272), Cartago (n=72), Heredia (n=38) y San José (n=34) (Figura 1). La gran mayoría de estos hatos se ubican en las faldas de los volcanes Poás, Barva e Irazú, a una elevación promedio de 1140 msnm (min=33, max=2965), contaban con un promedio de 14,3 años de seguimiento (min=1, max=26) y un promedio de 102 vacas con historial de producción registrada (min=6, max=777). Una vez realizada esta selección, se procedió al cálculo de estadísticos descriptivos de tendencia central (promedio, mediana) y medidas de dispersión (desviación estándar, intervalos de confianza), así como una exploración gráfica de la distribución de la variable kgl. Se realizó una revisión exhaustiva para detectar la existencia de posibles valores extremos (>3,5 desviaciones estándares dentro de raza y periodo de lactancia), que fueron eliminados de análisis subsecuentes.

Análisis de la asociación entre el índice de temperatura y humedad (ITH) y la producción de leche diaria (kgl)

Para explorar la asociación entre el ITH y la kgl, se ajustó un modelo lineal mixto generalizado con base en el procedimiento GLIMMIX (SAS Institute Inc., 2013). El modelo estadístico asumido fue el siguiente:

Y=β0+β1R+β2A+β3M+β4P+β5E+β6(ENOS)+β7(ITH)+β8(R×E)+β9(R×ITH)+v+ε

donde:

Y= variable dependiente, registros individuales de producción de leche diaria (kgl).

β0= intercepto.

β1R= efecto fijo ligado al grupo racial (clases: Holstein, Jersey y cruces Holstein×Jersey, que incluyó animales con diferentes proporciones de ambas razas, en su mayoría ½-½ y ¾-¼).

β2A= efecto fijo ligado al año calendario (clases: 1990…2015).

β3M= efecto fijo ligado al mes calendario (clases: enero, febrero,…diciembre).

β4P= efecto fijo ligado al número de parto (clases: 1, 2, 3, 4, 5, 6, ≥7).

β5E= efecto fijo ligado a la etapa de la lactancia (clasificado en categorías de amplitud 20 días: 20, 40…360).

β6(ENOS)= efecto fijo ligado a ENOS (clases: fase cálida-Niño, fase fría-Niña, normal).

β7(ITH)= efecto fijo ligado al índice de temperatura y humedad (ITH_máximo), el cual se obtuvo por predicción lineal y estocástica, según se describió en la sección anterior (clases: ≤60, 62,64. hasta ≥88, con categorías intermedias clasificadas en intervalos con amplitud de dos unidades).

β8(R×E)= efecto fijo ligado a la interacción de raza por etapa de lactancia (Ra_1/E_0….Ra_n/E_360).

β9(R×ITH)= efecto fijo ligado a la interacción de raza por ITH (Ra_1/ITH_60...Ra_n/ITH_88).

v= efecto aleatorio de la vaca, que considera la estructura de correlación existente entre registros diarios de producción para un mismo animal a lo largo de su vida productiva.

ε= error residual aleatorio.

El modelo anterior se ajustó de manera progresiva, con el fin de evaluar la consistencia de los estimados y evitar la sobreparametrización.

Resultados

Tendencias en índice de temperatura y humedad (ITH)

Los datos climáticos disponibles para el periodo analizado evidenciaron mayor temperatura, precipitación, humedad e ITH en la Región Norte, en comparación con la Región Central (p<0,0001) (Cuadro 1). Las diferencias fueron de 5,3, 7,7 y 6,0 °C en temperaturas promedio, mínima y máxima, respectivamente, 87 mms en precipitación mensual, 3,8% en humedad relativa, 9,1 en ITH promedio y 10,5 en ITH máximo. Por el contrario, el valor promedio de radiación fue mayor (1,9 W m-2) en la Región Central (p<0,0001). Las tendencias observadas para todas las variables climáticas disponibles fueron descritas en mayor detalle por Ruiz-Jaramillo (2017).

Cuadro 1. Promedio (), desviación estándar (D.E.), mínimo (Min) y máximo (Max) de variables climáticas para las regiones Central (ocho estaciones, 1630 registros) y Norte (cinco estaciones, 1117 registros). Costa Rica, periodo 1990-2015.

Table 1. Mean (), standard deviation (D.E.), minimum (Min) and maximum (Max) of climatic variables for the Central (eigth stations, 1630 records) and North (five stations, 1117 records) regions. Costa Rica, period 1990-2015.

La variable ITH (máximo) presentó una distribución normal, con dos picos pronunciados, claramente asociados a las dos regiones analizadas. El modelo estadístico ajustado sobre el ITH convergió a una solución satisfactoria, indicando que las variables con efectos altamente significativos (p<0,0001) fueron la elevación, el año y mes calendario, el efecto El Niño (ENOS) y la interacción región×mes. Por otra parte, las interacciones región×ENOS y región×año presentaron un efecto significativo (p<0,05). La variación entre y dentro de estaciones también fue altamente significativa (p<0,0001). La correlación entre valores observados de ITH y predichos por el modelo fue de 0,98 y el error estándar de los residuales fue de 1,51.

El mes calendario y la elevación (msnm) fueron las variables más estrechamente asociadas con el ITH. El modelo estimó una marcada relación lineal inversa entre la covariable de elevación y el ITH, con una reducción estimada de 0,92 puntos de ITH por cada incremento de 100 m en elevación (p<0,0001, r2= 0,78).

Tanto en la Región Norte como en la Central, las medias marginales del ITH a lo largo de los meses del año mostraron un patrón homogéneo, con dos curvas aproximadamente paralelas (Figura 2). En la Región Norte el ITH fue significativamente mayor (p<0,05) a lo largo de todos los meses, con medias marginales que variaron entre 85 y 89, mientras que en la Región Central variaron entre 75 y 79, para una diferencia de aproximadamente diez puntos entre regiones. Dentro de cada región, los meses más secos (diciembre a marzo), presentaron menor ITH (p<0,05) en comparación con los meses más húmedos (mayo/junio y agosto/setiembre). En ambas regiones, se observó un leve descenso del ITH en el mes de julio.

Figura 2. Medias marginales y bandas de confianza de 95 % para la variable índice de temperatura y humedad (ITH), en función de los meses del año y la región climática. Costa Rica, periodo 1990-2015.

Figure 2. Marginal means and 95 % confidence bands for the variable index of temperature and humidity (TIH) depending on the months of the year and the climatic region. Costa Rica, period 1990-2015.

En cuanto a la variación interanual también se observaron tendencias homogéneas entre regiones, pero con fluctuaciones irregulares, más evidentes alrededor de los años 1997/1998 y 2010/2011 (Figura 3). Un análisis de regresión lineal realizado sobre las medias marginales anuales del ITH por región no mostró tendencias consistentes, ya que se estimó un aumento anual no significativo de 0,018 puntos en ITH en la Región Norte (r2= 0,04, p=0,30), y una reducción anual no significativa de 0,017 en la Región Central (r2= 0,06, p=0,24).

Figura 3. Medias marginales y bandas de confianza de 95 % para la variable índice de temperatura y humedad (ITH) en función de los años y la región climática. Costa Rica, periodo 1990-2015.

Figure 3. Marginal means and 95 % confidence bands for the variable index of temperature and humidity(ITH) depending on the year and the climatic region. Costa Rica, period 1990-2015.

En cuanto al efecto ENOS, el modelo estimó medias marginales del ITH de 82,6, 82,5 y 81,9 para las fases cálida (El Niño), normal y fría (La Niña), respectivamente, con diferencias significativas (p<0,05) solo entre la fase fría y las dos restantes.

Predicciones del índice de temperatura y humedad (ITH) por hato/año/mes

Se obtuvieron predicciones lineales y estocásticas para un total de 47 934 clases de hato/año/mes. El modelo de autocorrelación espacial que mejor ajustó los residuales fue el sinusoidal (sine hole effect) cuyos parámetros optimizados fueron 2,20±0,003 (efecto pepita o estimado de la variabilidad entre mediciones realizadas en un mismo sitio), 0,19±0,003 (meseta) y 0,62±0,013 (rango). Se comprobó el cumplimiento adecuado del supuesto de anisotropía, ya que se obtuvieron patrones sinusoidales similares en las cuatro direcciones analizadas.

El promedio de estaciones circundantes utilizadas para obtener estas predicciones fue de 3,6 (Min=3, Max=8). La distancia promedio entre los sitios de los hatos y la estación meteorológica más próxima fue de 11,4 km (Min=0,32, Max=39,6). Del mismo modo, la distancia promedio a las tres estaciones más cercanas fue de 19,2 km (Min=4, Max=54,5). El 99 % de los hatos contaron con al menos una estación meteorológica dentro de un radio de 33,3 km del sitio.

Para un mismo año y mes, los valores del ITH predichos para los sitios de los hatos ubicados en la Región Central fueron altamente fluctuantes, oscilando entre 56 y 80, mientras que los hatos ubicados en la Región Norte presentaron valores de ITH más altos y estables, entre 80 y 84 (Figura 4).

Figura 4. Valores predichos de índice de temperatura y humedad (ITH) en las coordenadas geográficas (X= longitud, Y= latitud) de 220 hatos, ubicados en las regiones Central (azul) y Norte (rojo) de Costa Rica, contemporáneos al mes de enero, año 2015.

Figure 4. Predicted values of temperature index and humidity index (ITH) in the geographic coordinates (X= Longitude, Y= Latitude) of 220 herds located in the Central (blue) and North (red) regions of Costa Rica, contemporary to the month of January, year 2015.

En un mismo hato, los valores predichos del ITH de distintos meses y años también variaron marcadamente (Figura 5). El patrón intermensual concuerda con el observado previamente en la Figura 2, con diferentes magnitudes asociadas a la variación interanual.

Figura 5. Valores predichos de índice de temperatura y humedad (ITH) para un hato ubicado en la región Central de Costa Rica (elevación 2400 msnm), según mes (enero a diciembre) y año (1995, 2000, 2005, 2010 y 2015).

Figure 5. Predicted values of temperature index and humidity index (ITH) for a herd located in the Central Region of Costa Rica (elevation 2400 masl), according to month (January to December) and year (1995, 2000, 2005, 2010 and 2015).

Producción de leche diaria (kgl)

Se analizaron un total de 6 478 582 registros de producción diaria de leche (49 % Holstein, 16,7 % Holstein×Jersey, 34,3 % Jersey) procedentes de 418 hatos, 123 869 vacas y 362 815 lactancias. La variable kgl mostró una distribución aproximadamente normal, con un promedio general de 19,5 kg y una desviación estándar de 7,9 kg. De acuerdo con los resultados del modelo de regresión, todos los factores analizados presentaron un efecto altamente significativo (p<0,0001) sobre la variable kgl.

En cuanto al factor racial, se observaron diferencias significativas (p<0,05) entre todas las razas, con una mayor producción para la Holstein (20,1 kg), seguida por Holstein×Jersey (16,7 kg) y Jersey (14,9 kg). Las medias marginales por raza y etapa de la lactancia mostraron curvas de lactancia aproximadamente paralelas, con picos de producción diaria alrededor del día 40, en el orden de 26,5, 21,7 y 19,8 kg, respectivamente. Las medias marginales de producción según número de parto mostraron un patrón de incremento desde el primer parto (14 kg) hasta el cuarto parto (18,4 kg), para luego decaer hasta 16,8 kg en el parto 7 y posteriores.

El patrón intermensual indicó producciones marginales que variaron entre 16,2 kg en enero hasta 17,8 kg en junio, con producciones significativamente más altas (p<0,05) en los meses más lluviosos (Figura 6). En cuanto a la tendencia interanual de la producción, se observó un patrón de incremento escalonado desde 13 kg en 1990, hasta alrededor de 19,3 kg en los últimos años (Figura 7). Un análisis de regresión sobre las medias marginales indicó un incremento anual estimado de 0,23 kg (r2= 0,91, p<0,001), aunque se observaron periodos con fluctuaciones irregulares alrededor de los años 1996/1997, 2009/2010 y 2013/2014. En relación con el efecto ENOS, se obtuvieron diferencias significativas (p<0,05) en la producción diaria obtenida durante las fases cálida (17,5 kg), normal (17,2 kg) y fría (16,9 kg), respectivamente.

Figura 6. Medias marginales con intervalos de confianza de 95 %, para la variable producción diaria de leche (kgl), en función del mes calendario. Costa Rica. Periodo 1990-2015.

Figure 6. Marginal means with 95 % confidence intervals for the daily milk production variable (kgl) as a function of the calendar month. Costa Rica. Period 1990-2015.

Figura 7. Medias marginales con intervalos de confianza 95 %, para la variable producción diaria de leche (kgl), en función del año calendario. Costa Rica. Periodo 1990-2015.

Figure 7. Marginal means with 95 % confidence intervals for the daily milk production variable (kgl) as a function of the year. Costa Rica. Period 1990-2015.

En lo que se refiere al efecto de ITH sobre kgl, se pudo observar una tendencia global a la disminución en producción conforme aumentó el ITH (Figura 8). Esta reducción fue más pronunciada a partir de ITH superiores a 72. La reducción en el promedio global marginal de kgl fue desde 19,7 kg para ITH entre 60-64, hasta alrededor de 13,3 kg para ITH de 88. Un análisis de regresión lineal indicó una reducción de 0,26 kg (r2= 0,92, p<0,001) en kgl por cada unidad de aumento en ITH.

Figura 8. Medias marginales con intervalos de confianza 95 %, para la variable producción diaria de leche (kgl), en función del índice de temperatura y humedad (ITH). Costa Rica. Periodo 1990-2015.

Figure 8. Marginal means with 95 % confidence intervals for the daily milk production variable (kgl) as a function of the temperature and humidity index (ITH). Costa Rica. Period 1990-2015.

Al realizar el análisis dentro de cada grupo racial, se observó un patrón de reducción similar en producción asociado al incremento del ITH, con ligeras fluctuaciones (Figura 9). De acuerdo con un análisis de regresión lineal, en el rango de valores del ITH entre 60 y 72, se observó un impacto negativo pero no significativo (p>0,05) del ITH sobre la producción, con tasas de cambio de -0,14 (r2=0,37), +0,03 (r2=0,12) y -0,04 (r2=0,64) kg de leche diaria por cada unidad de incremento en el ITH para Holstein, Holstein×Jersey y Jersey, respectivamente. En el rango de ITH entre 74 y 88, el impacto negativo fue más consistente, con reducciones altamente significativas (p<0,0001) de 0,41 (r2=0,94), 0,36 (r2=0,95) y 0,29 (r2=0,82) kg de leche diaria por cada unidad de incremento en el ITH, respectivamente.

Figura 9. Medias marginales y tendencia lineal ajustada (con coeficiente de regresión β y coeficiente de determinación r2) para la variable producción diaria de leche (kgl) en función del índice de temperatura y humedad (ITH), según grupo racial. Costa Rica, periodo 1990-2015.

Figure 9. Marginal means and adjusted linear trend (with regression coefficient β and coefficient of determination r2) for the variable daily milk production (kgl) as a function of temperature index and humidity index (ITH), according to breed type. Costa Rica, period 1990-2015.

Discusión

En general, los valores observados para la variable índice de temperatura y humedad (ITH) son congruentes con lo esperado, en función de la ubicación geográfica y las condiciones topográfícas de las regiones incluidas en este estudio. La mayor temperatura, precipitación y humedad relativa presentes en la Región Norte ha sido ampliamente descrita en estudios previos de carácter meteorológico (IMN, 2008). Esto conlleva a diferencias de aproximadamente diez unidades en el ITH de las dos regiones estudiadas.

De acuerdo con la tabla de ITH vs. estrés calórico en ganado bovino, referida por Armstrong (1994), el ganado lechero de la Región Central podría estar sometido a un nivel de estrés moderado en gran parte del día. Bajo los niveles promedio de humedad relativa (82,7 %, Cuadro 1), esto sucedería cuando la temperatura supera los 23 °C. En la Región Norte, donde la humedad promedio es del 86,5 % (Cuadro 1), el nivel de estrés calórico podría llegar a considerarse como severo en gran parte del día, cuando las temperaturas exceden los 28 °C. Cabe aclarar que los valores obtenidos del ITH, al ser calculados con base en temperaturas máximas, reflejan el límite superior de estrés calórico al cual pueden estar sometidos los animales, lo cual sucede durante las horas de mayor temperatura, generalmente a principios de la tarde. Debido a la fórmula utilizada, el ITH estuvo influenciado principalmente por la temperatura, siendo menos sensible a la humedad relativa.

La Región Norte posee una gran influencia del litoral Atlántico (IMN, 2008), por lo que presenta mayor temperatura y humedad relativa, lo que incide en los mayores valores del ITH. Además, es una región que presenta menor variabilidad en elevación. La Región Central posee influencia de ambos litorales (IMN, 2008) y es una zona con topografía más irregular, por lo que presenta un promedio más bajo del ITH, pero un mayor rango de variación. La mayoría de los hatos lecheros incluidos en el presente estudio se ubican en las faldas montañosas, donde la temperatura e ITH tienden a ser menores.

Las fluctuaciones intermensuales observadas en el ITH (Figura 2) fueron altamente coincidentes con los patrones observados en temperatura y en menor grado con los mostrados para humedad relativa (Ruiz-Jaramillo, 2017). En el presente estudio, el descenso en el ITH en el mes de julio está ligado al fenómeno conocido como veranillo y la presencia de fuertes vientos alisios (IMN, 2008). Los resultados indicaron que la fluctuación intermensual en promedios del ITH puede ser de hasta cinco unidades y está ligada principalmente a las variaciones en temperatura. Estas fluctuaciones intermensuales son aún mayores a las interanuales, que generalmente se mantienen dentro del rango de dos unidades de variación.

Las fluctuaciones interanuales en ITH observadas en el presente estudio (Figura 3) parecen estar ligadas mayormente a la variabilidad climática normal, aunque se observaron variaciones más marcadas en los años 1997/1998 y 2010/2011, que concuerdan parcialmente con la presencia de fases acentuadas del evento ENOS (IMN, 2008). De acuerdo con registros históricos, los años con fases cálidas acentuadas fueron: 1987, 1991, 1992, 1997, 2002, 2009 y 2015, mientras que las fases frías acentuadas ocurrieron principalmente en los años: 1985, 1988, 1998, 1999, 2000, 2008, 2010 y 2011 (CPC, 2016).

Distintas investigaciones reportan que la magnitud y duración del evento ENOS no es la misma cada vez que se presenta, y su impacto depende de la susceptibilidad de la región al evento, ya sea a su fase fría o a su fase cálida, ya que este no se da aisladamente y puede modificarse por factores orográficos (IMN, 2008; Bonilla, 2014). En la Región Norte, se considera que la señal del efecto de ENOS no es robusta, por lo que su influencia ha sido descrita como muy variable (IMN, 2008).

Además del ENOS, en el periodo analizado ocurrieron varios eventos catalogados como extremos, los cuales también pueden contribuir a la variabilidad observada. Entre los eventos de mayor impacto en la región durante el periodo analizado se contabilizaron los huracanes César (1996) y Mitch (1998) y las tormentas tropicales Alma (2005) y Tomás (2010) (Vallejos et al., 2012). Estos eventos son generalmente de corta duración, por lo que su impacto en las tendencias interanuales es menos marcado.

Los datos disponibles no evidenciaron tendencias lineales significativas en el ITH interanual, dado que los estimados de regresión no fueron consistentes ni significativos en ambas regiones. Otros estudios tampoco han evidenciado tendencias significativas para la temperatura en estas dos regiones, aunque sí se han registrado aumentos significativos en la frecuencia de eventos extremos (IMN, 2008). La existencia de series climáticas incompletas y el uso de diferentes tecnologías en distintas estaciones hace difícil la comprobación de tendencias climáticas (IMN, 2008).

Los valores predichos del ITH para distintos hatos/año/mes (Figura 4), mostraron una mayor fluctuación en los hatos de la Región Central que en los de la Región Norte. Esto se debió principalmente a que los hatos lecheros de la Región Central se distribuyeron en un rango más amplio de elevación, desde 192 hasta 2965 msnm (promedio= 1663 msnm), mientras que, en la Región Norte, los hatos se ubicaron entre 33 y 1946 msnm (promedio= 383 msnm). La elevación, por su relación directa con la temperatura, fue la variable evaluada que tuvo mayor impacto sobre el ITH, por lo que su inclusión en el modelo predictivo demostró ser particularmente útil para lograr obtener estimados más precisos de esta variable a nivel de los sitios de los hatos.

En cuanto al comportamiento intermensual de la producción (Figura 6), la tendencia a mayor producción en los meses más lluviosos y cálidos, posiblemente estuvo ligada a una mayor disponibilidad de alimento. En una parte de los hatos, la disminución en precipitación en los meses más secos pudo limitar la disponibilidad del forraje y afectar la producción, principalmente en zonas bajas, con época seca más marcada y menor uso de suplementación. A nivel local, los niveles de suplementación de concentrado pueden variar considerablemente entre hatos. Un estudio realizado en 104 fincas de lechería especializada reportó que, del total de materia seca de la dieta consumida por las vacas lecheras, el porcentaje proveniente del concentrado osciló entre 23 % en el grupo de menor uso, hasta un 38 % en el grupo de mayor uso; mientras que el porcentaje proveniente del pasto varió desde un 18,8 % hasta un 71 % entre los diferentes grupos (Iñamagua-Uyaguari et al., 2016).

Las variaciones interanuales en producción (Figura 7) pueden estar asociadas con diferentes factores. La tendencia general de incremento en producción en el periodo analizado ha sido reportada por estudios previos (Vargas y Gamboa, 2008; Vargas, 2016), y se atribuye principalmente al mejoramiento genético y a las mejoras progresivas en las condiciones de manejo de los hatos lecheros. Las variaciones en la tendencia anual también podrían estar asociadas a eventos climáticos. En el presente estudio se evidenció reducción significativa en producción durante las fases frías del evento ENOS, no así durante las fases cálidas. La fase fría del ENOS es conocida por propiciar lluvias en el Pacífico, mientras que produce sequía en las regiones Central y Norte (IMN, 2008; Vallejos et al., 2012), lo que podría tener relación con este resultado.

El presente estudio evidenció una mayor producción de la raza Holstein por sobre el cruce Holstein×Jersey y la Jersey (Figura 9), mismo comportamiento que ha sido observado previamente a nivel local (Vargas y Gamboa, 2008; Vargas y Ulloa, 2008; Bolívar et al., 2009; Vargas-Leitón et al., 2012). Asimismo, la relación inversa observada entre ITH y producción en los tres grupos raciales (Figura 9) fue similar a lo reportado por otros estudios. El estrés calórico que se produce cuando el ITH es elevado afecta negativamente la producción de leche (Bouraoui et al., 2002; West et al., 2003; Bohmanova et al., 2007; Gantner et al., 2011; Hammami et al., 2013; Zewdu et al., 2014; Garcia et al., 2015; Pragna et al., 2017).

En un estudio poblacional similar al efectuado en el presente trabajo, en condiciones de clima subtropical, en vacas Holstein con niveles de producción mayores a 28 kg, se observó que la producción diaria de leche declinó a partir del ITH>72, con una tasa de reducción de 0,39 kg por cada unidad de incremento en el ITH (Bohmanova et al., 2007), muy similar al obtenido en el presente estudio para la misma raza. En condiciones de clima templado, un estudio similar con ganado Holstein de alta producción, reportó reducciones de 0,164 kg en producción diaria de leche para un ITH <62, 0,36 kg para un ITH entre 62-70 y de 0,96 kg para un ITH>70 (Hammami et al., 2013). En este caso, los autores atribuyeron la mayor susceptibilidad al ITH a una menor adaptabilidad de la raza Holstein al estrés calórico en climas templados. En Etiopía, otro estudio poblacional reportó una reducción de hasta 190 kg en producción de leche por lactancia de vacas cruzadas Holstein×Deoni por cada incremento de una unidad en el ITH máximo (Zewdu et al., 2014).

Otros estudios conducidos en condiciones semiexperimentales, con grupos reducidos de animales, también demostraron el efecto adverso del ITH sobre la producción (Bouraoui et al., 2002; West et al., 2003; Gantner et al., 2011). Un estudio experimental en vacas Holstein estabuladas con niveles de producción cercanos a 20 kg en condiciones de clima Mediterráneo, reportó una reducción de 0,41 kg de leche día-1 por cada unidad de incremento de ITH>69 (Bouraoui et al., 2002). En ganado estabulado Holstein con niveles de producción >16 kg se reportaron reducciones en producción diaria de 0,9 kg en vacas expuestas a un ITH >72 comparadas contra el grupo expuesto a un ITH≤72 (Gantner et al., 2011). En otro estudio, vacas estabuladas de razas Holstein (n=24, 24 kg leche día-1) y Jersey (n=8, 23 kg leche día-1), presentaron reducciones en el rango de 0,88 y 0,60 kg de producción diaria de leche, respectivamente, por cada unidad de incremento en el promedio del ITH de los dos días previos (West et al., 2003), valores que duplican los encontrados en el presente estudio.

Se ha reportado que la susceptibilidad de un individuo al estrés calórico depende de la raza, observándose mayores problemas en las vacas lecheras de alto valor genético, seleccionadas por su alta producción de leche (Ravagnolo y Misztal, 2000; Hammami et al., 2013; Dunn et al., 2014; González, 2014). La aparición de estrés térmico para las vacas Holstein de alta producción se encuentra a un menor umbral, quizás tan bajo como 65-69 de ITH (Bryant et al., 2007; Zimbelman et al., 2009; Hammami et al., 2013).

Los resultados del presente estudio sugieren efectos similares del ITH sobre los tres grupos raciales analizados. En términos económicos, la reducción en producción diaria de leche observada en las razas Holstein (β= 0,41), Holstein×Jersey (β= 0,36) y Jersey (β= 0,29), en el rango de ITH entre 75 (moderado) y 85 (severo), equivale a una diferencia de 1261, 1087 y 884 kg de leche en una lactancia estándar, con un costo local de $680, $577 y $477, respectivamente. Este es un costo que representa la pérdida por lactancia en que puede incurrir una vaca de los grupos raciales mencionados expuesta a estrés severo (p.e. Región Norte), en comparación con una expuesta a estrés moderado (p.e. región Central).

Todos los estudios coinciden en el efecto adverso del ITH, pese a que existen diferencias en cuanto a magnitud del efecto y los umbrales en que empiezan a manifestarse. Estas diferencias pueden estar ligadas a circunstancias particulares de cada estudio, tales como las razas implicadas y su nivel de producción, diferencias climáticas entre los sitios de estudio, diferentes sistemas de manejo y alojamiento, o inclusive detalles metodológicos, tales como el uso de valores diarios, medias móviles o promedios mensuales de ITH, así como la fórmula empleada para el cálculo del mismo.

Entre las ventajas del presente estudio están el uso de modelos predictivos para la estimación del ITH en los sitios exactos de los hatos durante todo el periodo bajo análisis, considerando todas las variables fijas y aleatorias de importancia, entre las cuales el mes y la elevación resultaron particularmente útiles. Es importante el uso de métodos de interpolación que ponderen la información proveniente de varias estaciones circundantes en función de su distancia al sitio, en vez de solo utilizar la más próxima. La predicción del ITH en el sitio de los hatos reviste importancia, sobretodo en el caso de regiones montañosas y de topografía muy variable, como las que se presentan en Costa Rica. Por otro lado, una limitante del presente estudio radica en el uso de promedios mensuales de ITH, ya que lo más adecuado sería el uso de valores correspondientes al mismo día del registro productivo, o medias móviles de los días inmediatos anteriores al registro de producción.

Conclusiones

El presente estudio aporta evidencia del efecto del estrés calórico sobre la producción del ganado lechero Bos taurus en las dos principales regiones lecheras de Costa Rica. Las tendencias observadas fueron consistentes y corroboran lo observado en otros estudios realizados en zonas subtropicales y templadas. El estudio evidencia una relación inversa entre la producción diaria de leche y el ITH para los tres grupos raciales analizados. Los efectos observados tienden a manifestarse a partir de ITH >72, como se evidencia en otros estudios. Las tendencias intermensuales observadas para el ITH indicaron que la Región Central se ubica mayormente en la franja de estrés moderado, mientras que la Región Norte se encuentra dentro de la franja de estrés severo.

El estudio muestra la factibilidad e importancia de integrar los datos meteorológicos para el estudio del impacto del clima sobre el rendimiento de las explotaciones ganaderas, y la necesidad de implementar de manera más dinámica esta información en la definición de políticas de manejo más adecuadas a nivel de hato, tendientes a prevenir o mitigar el efecto del estrés calórico.

Agradecimiento

El autor Hugo Hidalgo agradece al financiamiento ofrecido a su persona mediante los proyectos de la Universidad de Costa Rica: 805-B7-286 (UCREA), 805-B8-766 (Redes) y 805-B9-454 (Grupos).

Literatura citada

Armstrong, D. 1994. Heat stress interaction with shade and cooling. J. Dairy Sci. 77:2044-2050. doi:10.3168/jds.S0022-0302(94)77149-6

Bohmanova, J., I. Misztal, and J. Cole. 2007. Temperature-humidity indices as indicators of milk production losses due to heat stress. J. Dairy Sci. 90:1947-1956. doi:10.3168/jds.2006-513

Bolívar, D., J. Echevarry, L. Restrepo, y M. Muñoz. 2009. Productividad de vacas Jersey, Holstein y Jersey×Holstein en una zona de bosque húmedo montano bajo (Bh-MB). Livest. Res. Rural Dev. 21(6):80. http://www.lrrd.org/lrrd21/6/boli21080.htm (consultado 12 nov. 2016).

Bonilla, A. 2014. Patrones de sequía en Centroamérica. Su impacto en la producción de maíz y frijol y uso del Índice Normalizado de Precipitación para los Sistemas de Alerta Temprana. GWP Centroamérica, Tegucigalpa, HND. https://www.gwp.org/globalassets/global/gwp-cam_files/patrones-de-sequia_fin.pdf (consultado 11 jul. 2016).

Bouraoui, R., M. Lahmar, A. Majdoub, M. Djemali, and R. Belyea. 2002. The relationship of temperature-humidity index with milk production of dairy cows in a Mediterranean climate. Anim. Res. 51:479-491. doi:10.1051/animres:2002036

Bouroncle, C., P. Imbach, P. Läderach, B. Rodríguez, C. Medellín, E. Fung, M.R. Martínez, y C.I. Donatti. 2015. La agricultura de Costa Rica y el cambio climático: ¿Dónde están las prioridades para la adaptación? CGIAR, Montpellier, FRA. https://cgspace.cgiar.org/bitstream/handle/10568/45941/PB%20Costa%20Rica.pdf?sequence=7 (consultado 3 Jul. 2016).

Bryant, J.R., N. López-Villalobos, J.E. Pryce, C.W. Holmes, and D.L. Johnson. 2007. Quantifying the effect of thermal environment on production traits in three breeds of dairy cattle in New Zealand. N.Z.J. Agric. Res. 50:327-338. doi:10.1080/00288230709510301

Chilès, J.P., and N. Desassis. 2018. Fifty years of kriging. In: B.S. Daya et al., editors, Handbook of mathematical geosciences. Springer Nature, CHE. p. 589-612. doi:10.1007/978-3-319-78999-6_29

CPC (Climate Prediction Center). 2016. Cold & warm episodes by season. CPC, USA. http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml (accessed Aug. 3, 2016).

Dunn, R.J.H., N.E. Mead, K.M. Willet, and D.E. Parker. 2014. Analysis of heat stress in UK dairy cattle and impact on milk yields. Environ. Res. Lett. 9:064006. doi:10.1088/1748-9326/9/6/064006

Gantner, V., P. Mijić, K. Kuterovac, D. Solić, and R. Gantner. 2011. Temperature-humidity index values and their significance on the daily production of dairy cattle. Mljekarstvo 61(1):56-63.

Garcia, A., N. Angeli, L. Machado, F. de-Cardoso, and F. Gonzalez. 2015. Relationships between heat stress and metabolic and milk parameters in dairy cows in southern Brazil. Trop. Anim. Health Prod. 47:889-894. doi:10.1007/s11250-015-0804-9

González, J.M. 2014. El estrés calórico en bovinos. Universidad Especializada en Producción Animal (UDCA), ARG. http://www.produccion-animal.com.ar/etologia_y_bienestar/bienestar_en_bovinos/14-stres.pdf (consultado 13 nov. 2016).

Hammami, H., J. Bormann, N. M’hamdi, H. Montaldo, and N. Gengler. 2013. Evaluation of heat stress effects on production traits and somatic cell score of Holsteins in a temperate environment. J. Dairy Sci. 96:1844-1855. doi:10.3168/jds.2012-5947

IMN (Instituto Meteorológico Nacional). 2008. El clima, su variabilidad y cambio climático en Costa Rica. Segunda comunicación nacional. IMN, San José, CRI.

Iñamagua-Uyaguari, J.P., A. Jenet, L.G. Alarcón-Guerra, S.J. Vílchez-Mendoza, F. Casasola-Coto, y M.A. Wattiaux. 2016. Impactos económicos y ambientales de las estrategias de alimentación en lecherías de Costa Rica. Agron. Mesoam. 27:1-17. doi:10.15517/am.v27i1.21874

Kolovos, A. 2010. Everything in its place: efficient geostatistical analysis with SAS/STAT® spatial procedures. SAS Institute Inc., Cary, NC, USA. http://support.sas.com/resources/papers/proceedings10/337-2010.pdf (accessed Sep. 27, 2016).

NRC (National Research Council). 1971. A guide to environmental research on animals. National Academy of Sciences, WA, USA.

Ordaz, J., D. Ramírez, J. Mora, A. Acosta, y B. Serna. 2010. Costa Rica: efectos del cambio climático sobre la agricultura. CEPAL, México D.F., MEX.

Pragna, P., P. Archana, J. Aleena, V. Sejian, G. Krishnan, M. Bagath, A. Manimaran, V. Beena, E. Kurien, G. Varma, and R. Bhatta. 2017. Heat stress and dairy cow: Impact on both milk yield and composition. Int. J. Dairy Sci. 12:1-11. doi:10.3923/ijds.2017.1.11

Ravagnolo, O., and I. Misztal. 2000. Genetic component of heat stress in dairy cattle, parameter estimation. J. Dairy Sci. 83:2126-2130. doi:10.3168/jds.S0022-0302(00)75095-8

Ríos, N., y M. Ibrahim. 2008. Impactos del cambio climático sobre los recursos hídricos. Boletín Técnico No. 30. CATIE, Turrialba, CRI.

Romero, J., J. Rojas, y S. Estrada. 2011. El programa VAMPP Bovino, una herramienta para la toma de decisiones. Ventana Lechera 5:4-10.

Ruiz-Jaramillo, I. 2017. Variabilidad climática en las regiones Norte, Central y Caribe y su asociación con variables productivas en fincas lecheras costarricenses. Tesis M.Sc., Universidad Nacional, Heredia, CRI.

SAS Institute. 2013. SAS/STAT 9.4 User´s Guide. SAS Institute Inc., Cary, NC, USA.

St-Pierre, N.R., B. Cobanov, and G. Schnitkey. 2003. Economic losses from heat stress by US livestock industries. J. Dairy Sci. 86:E52-E77. doi:10.3168/jds.S0022-0302(03)74040-5

Vallejos, S., L. Esquivel, y M. Hidalgo. 2012. Histórico de desastres en Costa Rica (Febrero 1723 - Setiembre 2012). Biblioteca Virtual en Salud, San José, CRI. http://www.bvs.sa.cr/ambiente/4923.pdf (consultado 15 jul. 2016).

Vargas, B. 2016. Tendencias poblacionales por raza. Universidad Nacional, Heredia, CRI. http://www.medvet.una.ac.cr/posgrado/gen/tend/tendxraza1609.xlsx (consultado 15 dic. 2016).

Vargas, B., y G. Gamboa. 2008. Estimación de tendencias genéticas e interacción genotipo×ambiente en ganado lechero de Costa Rica. Téc. Pecu. Méx. 46:371-386.

Vargas, B., y J. Ulloa. 2008. Relación entre curvas de crecimiento y parámetros reproductivos en grupos raciales lecheros de distintas zonas agroecológicas de Costa Rica. Livest. Res. Rural Dev. 20(7):122. http://www.lrrd.org/lrrd20/7/varg20103.htm (consultado 18 oct. 2015).

Vargas-Leitón, B., Y. Marín-Marín, y J.J. Romero-Zúñiga. 2012. Comparación bioeconómica de grupos raciales Holstein, Jersey y Holstein×Jersey en Costa Rica. Agron. Mesoam. 23:329-342. doi:10.15517/am.v23i2.6533

Vargas-Leitón, B., O. Solís-Guzmán, F. Saénz-Segura, y H. León-Hidalgo. 2013. Caracterización y clasificación de hatos lecheros en Costa Rica mediante análisis multivariado. Agron. Mesoam. 24:257-275. doi:10.15517/am.v24i2.12525

West, J., B. Mullinix, and J. Bernard. 2003. Effects of hot, humid weather on milk temperature, dry matter intake, and milk yield of lactating dairy cows. J. Dairy Sci. 86:232-242. doi:10.3168/jds.S0022-0302(03)73602-9

Zewdu, W., B. Thombre, and D. Bainwad. 2014. Effect of macroclimatic factors on milk production and reproductive efficiency of Holstein Friesian×Deoni crossbred cows. J. Cell Anim. Biol. 8:51-60. doi:10.5897/JCAB2014.0408

Zimbelman, R., R. Rhoads, M. Rhoads, L. Baumgard, and R. Collier. 2009. A re-evaluation of the impact of temperature humidity index (THI) and black globe humidity index (BGHI) on milk production in high producing dairy cows. In: R.J. Collier, editor, Proc. 24th Annual Southwest Nutrition and Management Conference. University of Arizona, Tempe, AZ, USA. p. 158-169.

---