Introducción

La ganadería colombiana presenta bajos índices productivos que se han atribuido a bajas eficiencias reproductivas (ER), bajas tasas de crecimiento e incipientes programas de mejoramiento genético (Fedegan, 2006, citado por ^{Martínez-Rocha et al., 2011}). La baja ER aumenta los requerimientos de energía y los costos de alimentación para mantener la población parental (^{Elizalde & Riffel, 2012}); en consecuencia, compromete la competitividad del sistema, ya que no se obtiene un ternero por año y se aumenta la tasa de descarte de vientres (^{Walsh et al., 2011}).

Los factores que afectan la fertilidad son de orden nutricional, genético, ambiental, de gestión y las interacciones entre estos (^{Granja et al., 2012}; ^{Ossa et al., 2007}; ^{Vergara et al., 2008};). Por ello, su estudio es complejo y constituye un reto para los investigadores, quienes deben entender los procesos biológicos, los órganos y los tejidos del animal para encontrar las causas de la baja eficiencia y proponer estrategias de solución (^{Walsh et al., 2011}).

Los bovinos criollos Romosinuano provienen de un hato conformado por las razas gallega, berrenda y andaluz (^{Rouse, 1977}), que fueron introducidas a Colombia en el año 1525 (^{Ossa et al., 2013}) y se formaron gracias a algunos eventos biológicos como la selección natural (^{Ossa et al., 2011}). Este hecho permitió que esta raza adquiriera caracteres adaptativos de gran importancia económica, entre los que se destaca su ER (^{Martínez-Villate et al., 2009}).

En este sentido, el intervalo entre partos (IEP) constituye el carácter más importante en la evaluación de la ER individual de las hembras dentro de un hato (^{Ossa et al., 2006}). Actualmente, se reconoce la importancia del IEP para la economía y la eficiencia de la producción bovina al incorporarse en las evaluaciones genéticas de los programas de selección en los países donde se desarrollan planes de mejoramiento genético.

Es sabido que la forma más común de evaluar cómo una población responde a la selección es mediante la estimación de la heredabilidad del carácter, que es la relación entre la varianza genética aditiva y la varianza fenotípica (^{De Villemereuil, 2018}). Para estimar la varianza genética aditiva en los análisis modernos se suelen implementar los supuestos estadísticos del modelo infinitesimal desarrollado por ^{Ronald Fisher (1918)} mediante el modelo lineal mixto (^{De Villemereuil, 2020}).

Fisher desarrolló la teoría del modelo infinitesimal bajo un supuesto simplista que prevalece hasta nuestros días: la variación genética de una variable fenotípica está determinada por un gran número de factores mendelianos, cada uno de los cuales (de efecto infinitesimal) realiza una pequeña contribución a la variación total (^{Nelson et al., 2013}; ^{Turelli, 2017}). Esto generalmente resulta en la distribución gaussiana de los efectos genéticos de los individuos y de la variable fenotípica (^{Blasco, s. f.}; ^{Hill, 2014}; ^{De Villemereuil, 2018}).

Si bien esta recurrencia a la ley normal tiene muchas ventajas de cómputo e interpretación (^{Blasco, s. f.}), ¿qué sucedería si la distribución de la variable fenotípica no es de tipo gaussiana? Posiblemente, al usar herramientas estadísticas que no reflejen con precisión su verdadera distribución, se generarían conclusiones erróneas (^{De Villemereuil, 2018}). En respuesta a esto, durante el último decenio se han desarrollado metodologías compatibles con las suposiciones centrales del modelo infinitesimal, que se adaptan a una gama más amplia de propiedades de la variable fenotípica y no necesariamente están sujetas a este tipo de restricciones. Así, estas metodologías permiten estudiar variables con una distribución distinta a la normal (^{De Villemereuil et al., 2016}; ^{De Villemereuil, 2018}).

En este sentido, una práctica común en el mejoramiento genético animal consiste en emplear modelos lineales tradicionales diseñados a partir del modelo animal para estimar los parámetros genéticos (^{Galeano & Manrique, 2010}; ^{Martínez-Rocha et al., 2011}; ^{Vergara et al., 2016}). No obstante, en ocasiones, dichos modelos son asumidos sin un estricto cumplimiento del fundamento de que, si el número de observaciones es elevado, los datos pueden llegar a ajustarse a la normalidad (^{Bandera & Pérez, 2018}).

Así mismo, en estudios genéticos es frecuente el análisis de datos correlacionados provenientes de mediciones repetidas en los mismos individuos, como es el caso del IEP evaluado en los trabajos anteriormente citados. En este contexto no es recomendable usar modelos lineales clásicos, ya que al realizar mediciones repetidas se incumplen los supuestos requeridos para su utilización, entre ellos los de independencia, normalidad y homocedasticidad (^{Bandera & Pérez, 2018}).

La estimación de los componentes de varianza y heredabilidad de variables fenotípicas no gaussianas mediante nuevas metodologías será útil en estudios genéticos cuantitativos que pretendan evaluar dichos parámetros en una especie animal. Por lo tanto, el objetivo del presente estudio fue estimar los componentes de varianza y heredabilidad para la variable IEP en vacas de la raza Romosinuano mediante el enfoque del modelo lineal mixto generalizado (MLMG) y respetando la distribución en la que esta variable tiende a expresarse. Además, se buscó determinar los factores ambientales responsables de su variación para evaluar la efectividad de los programas de conservación o de mejoramiento genético.

Materiales y métodos

Descripción de los datos

Con el fin de obtener información veraz y consistente, se realizó un control de calidad inicial de 15.021 registros de partos. En el análisis se incluyeron los registros que presentaban información completa sobre la identificación de cada hembra y se observó de forma minuciosa que los partos fueran consecutivos, congruentes con la fecha de parto y los partos anteriores y posteriores. La depuración de la base de datos se realizó por medio del paquete dplyr (^{Wickham et al., 2019b}) del software R 3.6.0 (^{R Core Team, 2019}). Las figuras fueron realizadas usando el paquete ggplot2 (^{Wickham et al., 2019a}) del lenguaje de programación R.

Luego de editar la información, el número de registros de IEP útiles fue de 7.560. La información de pedigrí para el IEP consistió de 1.962 animales, incluyendo toros, vacas y sus crías. Por otra parte, 190 (73 %) del total de los toros tuvieron dos o más hijas, mientras que solo 68 (27 %) tuvieron una sola hija. El número total de vacas y sus crías fue de 925 y 1.190, respectivamente.

Las variables mes, año y orden de parto fueron clasificadas en variables tipo dummy o variables clasificadoras. Con el mes de parto se formaron cuatro clases de acuerdo con la variación de los registros de precipitación del Centro de Investigación Turipaná (lluvias, seca y transiciones), y el año de parto (periodo 1936-2016) se agrupó en periodos de diez años que conformaron ocho clases de año de parto. Esto se realizó con el fin de obtener una mejor representación del IEP en cada clase de mes y año de parto, pues al considerar cada variable sin formar las clases mencionadas la información era escasa, lo cual podría afectar la convergencia del modelo a ajustar. Siguiendo esta misma lógica, en cuanto al orden de IEP y exceptuando los partos del primero al sexto, fue necesario agrupar el orden de IEP siete y ocho en una misma clase (IEP 7) y los IEP iguales y mayores a nueve en una misma clase (IEP 8) (tabla 1, figura 1).

Tabla 1. Descripción de las variables de estudio usadas en la modelación 

Fuente: Elaboración propia

Con el fin de identificar la distribución para modelar la variable IEP, se calcularon las estadísticas descriptivas (asimetría y curtosis) de esta variable una vez editados los datos (figura 1A). Asimismo, se calculó el criterio de Cramér-von Mises, la prueba de Kolmogorov-Smirnov y la prueba de Anderson- Darling, usando el paquete fitdistrplus (^{Delignette-Muller et al., 2019}) del software R para el análisis estadístico de datos. Todo ello llevó a considerar el ajuste del modelo mediante el uso de una distribución asimétrica a la derecha: la distribución log-normal, al resultar, por ejemplo, en una asimetría igual a 1,26 y una curtosis de 5,50.

Fuente: Elaboración propia

Figura 1. Distribución de cada una de las variables de estudio.Nota: La variable para la cual se muestra una densidad (intervalo entre partos) se considera continua y las variables para las cuales se muestra un histograma (mes de parto, año de parto y orden de IEP) se consideran discretas. 

Modelo lineal mixto generalizado

Antes de ajustar el modelo estadístico, se verificó si existía colinealidad entre las variables de estudio (tabla 1), dado que esta puede causar problemas en la interpretación y la convergencia del modelo, y afectar la precisión de la heredabilidad a estimar (^{Harrison et al., 2018}). Para esto, se utilizó el factor de inflación de varianza generalizada (VIF, por su sigla en inglés), implementado en el paquete car (^{Fox & Weisberg, 2019}) del lenguaje de programación R. Como umbral se consideró un valor igual a 10; un valor por debajo de este es indicio de que el modelo a ajustar no implica problemas de colinealidad (^{Ringen et al., 2019}).

Dado que la variable es un conjunto de medidas repetidas tomadas a través del tiempo con distribución no gaussiana, esta se estudió mediante modelos estadísticos especiales que permitieran flexibilizar la distribución de la variable fenotípica. Por tanto, se utilizó el MLMG como variable fenotípica distribuida de forma log-normal y se permitió que el intercepto y la pendiente variaran entre los distintos individuos. De esta forma, al restringir que cada unidad no independiente (cada uno de los individuos con sus respectivos IEP) tuviera el mismo intercepto o pendiente (coeficientes de regresión asociados a cada individuo), fue posible tener en cuenta tanto las estructuras de asociación entre cada IEP del individuo (^{Harrison et al., 2018}) como la variabilidad aportada por su componente genético. Así, se logró representar el modo en que los datos fueron obtenidos y los aspectos biológicos de la variable fenotípica (^{Martínez Niño et al., 2012}).

La forma general del MLMG, que especifica las variables ambientales que pueden influir sobre la variable fenotípica a analizar (IEP), se describe en la ecuación 1.

Ecuación 1

Donde 𝑦_𝑖𝑗 y 𝑋_𝑖𝑗 son, respectivamente, el IEP y las variables ambientales (incluidas las clases mes de parto, año de parto y orden de IEP) para cada observación 𝑖 en el individuo 𝑗. El desvío aleatorio del individuo 𝑗 respecto al intercepto global fijo (𝛽₀) está representado por 𝛼_0𝑗, mientras que el desvío aleatorio de la pendiente global fija (𝛽₁) está representado por 𝛼1𝑗. Debido a que el intercepto y la pendiente están típicamente correlacionadas, se requieren parámetros adicionales para modelar el efecto aleatorio, representados por la matriz de covarianza 𝛴, que es igual a . Los elementos de la diagonal de 𝛴 son la varianza del intercepto aleatorio (𝜎_𝛼0²) y la varianza de la pendiente aleatoria (𝜎_𝛼1²), mientras que los elementos fuera de la diagonal corresponden a la covarianza (𝜎_𝛼0𝛼1²)entre el intercepto y la pendiente aleatoria. Por último, 𝜀𝑖𝑗 es el residuo de la observación 𝑖 del individuo 𝑗.

Para implementar el modelo anterior, se empleó el paquete brms (^{Bürkner, 2019}) del lenguaje de programación R, el cual utiliza el método Hamiltonian Monte Carlo y su extensión, el muestreador de No-U-Turn, que, además de alcanzar una buena convergencia del modelo, permite obtener resultados razonables incluso con unos pocos miles de muestras (^{Bürkner, 2017}). En este sentido, el modelo se ajustó usando dos cadenas de 2.000 iteraciones cada una, 500 de las cuales se utilizaron como calentamiento para un total de 3.000 muestras posteriores. En el proceso de ajuste del modelo se asumió que los parámetros a estimar presentaron distribución a priori uniforme; es decir, se optó por emplear distribuciones a priori no informativas.

Por último, se implementaron los denominados intervalos de credibilidad o intervalos bajo el paradigma bayesiano, para dar respuesta a la pregunta de si las diferencias observadas entre cada una de las clases formadas de determinada variable ambiental son relevantes para concluir si dicha variable es responsable de la variación en el valor de la variable IEP. Para ello, se construyó una hipótesis de equivalencia según la cual si el intervalo de credibilidad obtenido contiene el cero dentro de un conjunto diferente de valores plausibles para la diferencia en los coeficientes de regresión (𝛽), estimado en determinado nivel de la variable ambiental, se asume la hipótesis nula (𝐻₀) como cierta y se concluye, de este modo, que la variable ambiental no es una fuente de variación importante. En caso de no incluir el cero, se considera válida la hipótesis alterna (𝐻_𝑎), lo cual indica que la variable ambiental es una fuente importante de variación.

Modelo genético cuantitativo

Se calculó el parámetro de heredabilidad (ℎ²) para la variable IEP en vacas de la raza Romosinuano. Este se definió matemáticamente como ℎ² = 𝜎²⁄𝜎², donde 𝜎² es la varianza genética aditiva, 𝜎² es la 𝑎 𝑝 𝑎 𝑒 varianza residual o ambiental y 𝜎² = 𝜎² + 𝜎² corresponde a la varianza fenotípica total atribuible a la 𝑝 𝑎 𝑒 variación entre individuos. Para este cálculo, se empleó el paquete QGglmm (^{De Villemereuil, 2020}) del lenguaje de programación R, especializado en el análisis estadístico de datos. Como valores de entrada, se usaron los componentes de varianza estimados mediante el uso del MLMG anteriormente descrito, con una distribución log-normal.

La forma general (matricial) del MLMG usado para el análisis genético (^{Bennewitz et al., 2014}) se describe en la ecuación 2.

Ecuación 2a

Ecuación 2b

Ecuación 2c

En la ecuación 2a, la variable 𝑙 (escala latente de los datos) representa la variación del IEP atribuible a efectos genéticos aditivos y efectos ambientales; 𝜇 es el valor medio del IEP de la población evaluada; 𝛽 es un vector de efectos fijos (incluidas las clases mes de parto, año de parto y orden de IEP), que relaciona dichos efectos con las observaciones de IEP (𝑙) mediante la matriz de diseño de efectos fijos 𝑋; 𝛼 es un vector aleatorio de efectos genéticos aditivos individuales (valores genéticos), que relaciona dichos efectos con las observaciones de IEP (𝑙) mediante la matriz de diseño 𝑍, y 𝜀 es un vector de efectos residuales.

Se asume que la estructura de covarianza del efecto aleatorio es 𝑣𝑎𝑟(𝑎) = 𝐴𝜎², donde 𝜎² es la varianza _{𝑎 𝑎} genética aditiva y 𝐴 es la matriz de relación genética aditiva derivada del pedigrí.

En la ecuación 2b, la variable latente 𝑙 se transforma luego a una escala distinta 𝜂 (escala esperada de los datos) por medio de la función de enlace inversa 𝑔⁻¹ . Posteriormente, en la ecuación 2c se especifica que el IEP se distribuyó alrededor de los valores esperados 𝜂 de acuerdo con la distribución log-normal (𝐷). El parámetro genético de heredabilidad fue estimado con base en la ecuación 2c, dado que la distribución log-normal corresponde a la distribución en la que la variable IEP se expresó.

La precisión de la estimación de la heredabilidad (ℎ²) dada por el error estándar (𝐸𝐸) se calculó con base en la ecuación 3 (^{Swiger et al., 1964}).

Ecuación 3

Donde 𝑛 es el número total de crías; 𝑘 es el número ponderado de crías por reproductor dada la expresión 1⁄(𝑅 − 1) 𝑥[𝑛 − (∑𝑛² _i.⁄𝑛)] —siendo 𝑅 el número de reproductores, 𝑛_i el número total de crías y 𝑛_𝑖 el número de crías por reproductor—; 𝑠 es el número total de reproductores, y 𝑡 es el coeficiente de correlación intraclase dada la expresión (1⁄4)𝑥ℎ².

Resultados y discusión

En este estudio, el promedio general del IEP en vacas Romosinuano fue de 379,74 ± 49,22 días, con un coeficiente de variación de 12,96 %. Este valor fue menor a los obtenidos en el mismo hato por ^{Ossa et al. (2007)} (433 ± 139 días) y por ^{Ossa et al. (2013)} (422 ± 131 días), con información del IEP registrada en los periodos 1981-2003 y 1936-1979, respectivamente. Por otro lado, se han reportado mayores IEP en razas criollas colombianas distintas a la Romosinuano, con valores de 439,42 ± 92,24 días (^{Martínez- Villate et al., 2009}) y 487 ± 147 días (^{Martínez-Rocha et al., 2011}) en las razas Sanmartinero y Blanco Orejinegro, respectivamente. Además, en razas autóctonas españolas como la rubia gallega y la asturiana de los valles se han registrado valores promedios de IEP de 488 ± 177 días y 409 ± 62,17 días (^{Gutiérrez et al., 2003}; ^{Yagüe et al., 2007}), respectivamente.

Ante la diferencia entre el valor del IEP hallado en este estudio y los reportados en la literatura, se evidencia la importancia de describir los cruces genéticos y el comportamiento en las condiciones particulares en las que se desempeña cada raza (^{Quijano & Montoya, 2003}). Como se demuestra en los 81 años de este estudio, no es un mito la alta ER de las hembras Romosinuano.

Como toda variable fenotípica, el IEP es una medida cuyo valor puede variar por causas ajenas a la genética (factores ambientales), lo cual impide atribuir dicho comportamiento reproductivo al efecto exclusivo de los genes (^{Blasco, s. f.}; ^{Flint & Woolliams, 2008}). Identificar las variables ambientales responsables de la variación del IEP es una estrategia que permite establecer planes adecuados para el mejoramiento de la ER, ya que el IEP es su principal componente (^{Ossa et al., 2006}).

Dado que el intervalo de credibilidad entre las clases mes de parto 1 y 2, 2 y 4, y 3 y 4 se hallan en un entorno de cero (tabla 2), los valores del IEP en dichos niveles se consideran en esencia equivalentes. No obstante, al observar que el valor del intervalo de credibilidad entre las clases mes de parto 1 y 3 [0,02; 0,09], 1 y 4 [0,01; 0,06], y 2 y 3 [-0,08; -0,02] no contienen el cero, fue posible concluir que la variable clase mes de parto fue una causa importante de la variación presentada en el valor del IEP. A partir de estimaciones basadas en el método de cuadrados mínimos y en el análisis de varianza —el cual supone una distribución normal—, autores como ^{Bergmann et al. (1998)}, ^{Cavalcante et al. (2000)}, ^{Dias et al. (1994)}, ^{McManus et al. (2002)}, ^{Padua et al. (1994)} y ^{Pereira et al. (1980)} han encontrado, al igual que este estudio, un efecto significativo del mes de parto sobre el IEP.

Tabla 2. Intervalos de credibilidad del 95% para la diferencia en los coeficientes de regresión (𝛽) de los distintos niveles de la variable clase mes de parto 

^Nota:Intervalos de credibilidad obtenidos bajo hipótesis de equivalencia (𝐻₀ = 𝛽₁ − 𝛽₂ = 0 : ⇔ 𝛽₁ = 𝛽₂ vs 𝐻_𝐴 = 𝛽₁ − 𝛽₂ ≠ 0 : ⇔ 𝛽₁ ≠ 𝛽₂).

Fuente: Elaboración propia

En el análisis descriptivo del IEP en la clase mes de parto (figura 2), el mejor comportamiento reproductivo (menor valor de IEP) se presentó en las vacas cuyo parto ocurrió entre los meses abril-junio (clase 2), con un valor promedio IEP de 367,12 ± 46,06 días. Por el contrario, en las vacas que lograron parir en los tres últimos meses del año (clase 4), se obtuvo un valor promedio de 388,92 ± 57,96 días. Esto puede deberse al periodo de preñez; dado que el periodo de gestación en vacas Romosinuano es de aproximadamente 284 días (^{Ossa et al., 2006}), se deduce que estas quedaron preñadas entre finales de diciembre del año anterior y en el periodo de enero a mediados de marzo del mismo año del parto, tiempo en el cual se presenta un forraje con cantidad y calidad menores que se atribuye a la baja precipitación. A esto se suma la pérdida de peso que acompaña el inicio del periodo de amamantamiento de las crías, que puede incrementar los días abiertos y, por ende, el IEP asociado.

Fuente: Elaboración propia

Figura 2. Análisis descriptivo para el intervalo entre partos según la clase mes de parto. 1: enero a marzo; 2: abril-junio; 3: julio-septiembre; 4: octubre-diciembre. * Media general (en días) del intervalo entre partos en cada una de las clases formadas 

Se evidenciaron cambios en la media y varianza para el IEP según la clase año de parto. El valor más alto se observó entre los años 1946-1955 (clase 2), con una media de 387,69 ± 56,75 días, mientras que el menor valor correspondió a la media entre los años 1956-1965 (clase 3), con un valor de 370,48 ± 50,55 días. ^{Según Gómez et al. (1996)}, ^{Martínez-Villate et al. (2009)}, ^{Martínez-Rocha et al. (2011)}, ^{Ossa et al. (2006)}, ^{Tobón et al. (1996)} y ^{Quijano y Montoya et al. (2003)}, la diferencia en el comportamiento del IEP entre años podría atribuirse a los efectos ambientales (precipitación, cantidad y distribución de lluvias) que pueden incidir directamente en la cantidad y la calidad del forraje —favoreciendo o no el desempeño reproductivo de los bovinos— y a los cambios de alimentación y manejo a que los bovinos son sometidos. Las diferencias de comportamiento entre los decenios pueden obedecer también al grado de descarte anual de las vacas por distintas razones.

Como se puede observar a partir de los intervalos de credibilidad para la diferencia de los coeficientes de regresión estimados en cada clase de orden de IEP (figura 3), esta diferencia se puede considerar como una causa importante de la variación en la duración del IEP en vacas de la raza Romosinuano. Al comparar los coeficientes de regresión entre cada clase de orden de IEP, se identificaron diferencias entre la clase 1 (vacas con su primer IEP) y las restantes (a excepción de la clase orden de IEP 8), y entre las vacas de la clase 2 y las de las clases 3 y 4.

Fuente: Elaboración propia

Figura 3. Intervalos de credibilidad del 95 % para la diferencia en los coeficientes de regresión (𝛽) de los distintos niveles de la variable clase orden de IEP.Nota: Cuando el intervalo de credibilidad no contiene el cero, no hay diferencia entre un par de niveles. Los segmentos en amarillo son los niveles que, bajo hipótesis de equivalencia, resultaron ser diferentes. El significado de cada una de las clases formadas para el año de parto se describe en la tabla 1. 

Los mayores IEP observados (figura 4) se presentaron en las vacas de primero y segundo IEP, con valores de 389,70 ± 49,38 días y 382,18 ± 51,19 días, respectivamente. A partir de la clase orden de IEP 3, se evidenció una disminución en dicha variable, que alcanzó su valor más bajo en el cuarto IEP (373,93 ± 46,93 días). En general, la mayoría de los estudios indican que las vacas de primer IEP presentan un mayor valor, debido a que en vacas primerizas el inicio del ciclo reproductivo suele coincidir con su periodo de crecimiento, condición que afecta el proceso fisiológico de la hembra durante su primer parto (^{Martínez-Villate et al., 2009}). Adicionalmente, la etapa de desarrollo pospuberal y de primer parto de la hembra bovina coincide con la muda de sus dientes, lo que disminuye su capacidad de pastoreo y afecta de forma negativa su ER. Una vez que las vacas alcanzan su madurez fisiológica (cerca de los cuatro o cinco años), todos sus procesos se regulan y su IEP tiende a normalizarse (^{Martínez-Villate et al., 2009}).

En relación con este último resultado se recomiendan mayores cuidados, especialmente en el manejo y la alimentación de las vacas con menos de tres partos (uno y dos IEP). Esto con el fin de reducir su IEP al disminuir los parámetros de intervalo parto-concepción o días vacíos (^{Gutiérrez et al., 2012}; ^{Sánchez, 2010}) y, de esta forma, acercar dicha variable fenotípica a la media del estudio. Este valor demuestra la gran capacidad de adaptación de la raza Romosinuano a su nicho de desarrollo, como es el valle del Sinú en el departamento de Córdoba, Colombia.

Fuente: Elaboración propia

Figura 4. Análisis descriptivo para el intervalo entre partos según la clase orden de IEP. 1: vacas de primer IEP; 2: vacas de segundo IEP; 3: vacas de tercer IEP; 4: vacas de cuarto IEP; 5: vacas de quinto IEP; 6: vacas de sexto IEP; 7: vacas de sexto y séptimo IEP; 8: vacas con nueve o más IEP. * Media general (en días) del intervalo entre partos en cada una de las clases formadas. 

Como es sabido, la mayoría de los métodos estadísticos aplicados a las investigaciones en genética cuantitativa suponen que la distribución de la variable fenotípica es de tipo normal. Un ejemplo de esta clase de variables son el peso de la vaca a distintas edades y el índice de conversión. No obstante, existe otro tipo de variables que, si bien son de medida continua como las ya mencionadas, tienen una distribución que no sigue el supuesto clásico de normalidad, como es el caso del intervalo entre partos. Para estas variables es aconsejable usar métodos más adecuados que permitan estudiar su genética sin basarse en supuestos que pueden no ser ciertos.

Gracias a los avances tecnológicos recientes y a las herramientas estadísticas adaptadas a una gama más amplia de propiedades de la variable fenotípica, como el MLMG, las preguntas de investigación y el análisis de los datos ya no se limitan al supuesto clásico descrito con anterioridad (^{Bolker et al., 2009}; ^{Harrison et al., 2018}; ^{Kain et al., 2015}). La metodología del MLMG permite, además de identificar las variables ambientales responsables de la variación de la variable fenotípica, determinar si la diferencia entre individuos respecto a dicha variable se debe al conjunto de genes heredados o al ambiente en el que estos se desarrollan.

La heredabilidad estimada para el IEP en vacas de la raza Romosinuano mediante el MLMG fue de 0,0001 ± 0,002 (tabla 3), un valor bajo que indica que la variación de esta variable fenotípica en la población bovina analizada se debe primordialmente al componente ambiental y no al genético. Este valor fue menor y similar al estimado por ^{Ossa et al. (2007)} en la raza Romosinuano (0,05 ± 0,04), ^{López de Torre y Brinks (1990)} en ganado de carne de la raza Retinta (0,0), ^{Vergara (2001)} en ganado Cebú (0,00 ± 0,06), ^{Montes et al. (2009)} en vacas Cebú comercial (0,096 ± 0,08), ^{Wasike et al. (2009)} en vacas de la raza Boram (0,00), ^{Forni y De Albuquerque (2005)} (0,04 ± 0,6) y ^{Grossi et al. (2009)} (0,04 ± 0,06) en vacas de la raza Nelore, ^{Vergara et al. (2016)} en ganado cruzado Angus con Brahman (0,003 ± 0,046 y 0,006 ± 0,013), ^{Bekele et al. (2017)} en ganado Boram (0,00 ± 0,03), y ^{Lopez et al. (2019)} en ganado negro japonés Wagyu (0,01 ± 0,07). Todos estos resultados se obtuvieron mediante la metodología del MLMG.

Por otro lado, valores más altos de heredabilidad para el IEP han sido reportados en vacas de la raza Blanco Orejinegro (0,13 ± 0,031) por ^{Martínez-Rocha et al. (2011)} y en bovinos de esta última raza cruzados con Cebú y Angus (0,15 ± 0,07) por ^{Vergara et al. (2008)}. Así mismo, en un estudio realizado con animales puros y cruzados de las razas Cebú y Angus, ^{Vergara et al. (2009)} reportaron un valor de heredabilidad de 0,11. Según ^{Falconer y Mackay (2001)}, generalmente el valor de heredabilidad de las variables reproductivas es baja, como se evidenció en el presente estudio, y resulta menor en poblaciones bovinas del trópico, cuyas condiciones de manejo presentan un bajo control de los efectos ambientales.

Tabla 3. Estimados de componentes de varianza y heredabilidad para la variable IEP de vacas raza Romosinuano 

^Nota:El número promedio de crías por reproductor (𝑘) se estimó a partir de la fórmula 1⁄(𝑅 − 1) 𝑥[𝑛 − (∑𝑛²⁄𝑛)], propuesta por ^{Ossa (2017)}, donde 𝑅 es el número de reproductores, 𝑛 el número total de crías y 𝑛𝑖 el número de crías por reproductor.

Fuente: Elaboración propia

El intervalo de credibilidad del 95 % para el valor de heredabilidad de la variable IEP (figura 5), el cual indica el conjunto de valores plausibles para este parámetro genético, demostró, como en la estimación anterior, que la variación de esta variable fenotípica es causada en su mayoría por factores ambientales. Por lo tanto, el mejoramiento fenotípico de la variable dependería básicamente de las mejoras en las condiciones de alimentación y manejo de la población evaluada. No obstante, ^{Buzanskas et al. (2010)} sugirieron incluir el IEP en programas de mejora genética aun cuando presente baja heredabilidad, ya que una mejor reproducción genera grandes beneficios para los sistemas de producción de carne bovina.

Fuente: Elaboración propia

Figura 5. Distribución posterior del valor de heredabilidad obtenido para la variable IEP de vacas raza Romosinuano e incertidumbre asociada a su estimación. Nota: La incertidumbre está representada por el intervalo de credibilidad del 95 % obtenido mediante el enfoque bayesiano, cuyo valor esta entre 2,48 × 10-07 y 0,013. Por lo tanto, es posible afirmar que se tiene un 95 % de confianza de que el rango [2,48 × 10-07; 0,013] contiene el valor real de la heredabilidad. 

El uso del MLMG en este estudio no mostró mayores diferencias con otros estudios donde se aplicó la metodología del modelo lineal mixto, lo que sugiere su potencial para el estudio de otras poblaciones.

Conclusiones

El IEP de 379,74 ± 49,22 días es un indicativo de la adaptación de la raza Romosinuano a su lugar de origen, que se debe a su alta eficiencia reproductiva. La heredabilidad de la población evaluada fue baja, lo que demostró que esta variable fenotípica es afectada por factores ambientales y no genéticos. Por otra parte, los efectos ambientales (mes de parto, año de parto y orden de IEP) influyen sobre la variabilidad del IEP, por lo que es necesario establecer estrategias de alimentación y manejo que reduzcan su valor hasta lograr el ideal de una cría por año. El modelo lineal mixto generalizado para estimar los componentes de varianza y heredabilidad de la variable IEP permitió considerar la distribución presentada por dicha variable (distribución log-normal); se recomienda su uso en análisis que midan la variable fenotípica en múltiples ocasiones y no sigan el supuesto clásico de normalidad.