All animals sourced in this study belonged to Auburn University. All procedures with animals were performed in accordance with the protocols approved by Institutional Animal and Care and Use Committee in Auburn University.
Overall nutritional management of heifers
The dataset used in this study contained the first breeding season pregnancy outcome, phenotypic, and pedigree records for crossbred, beef heifers (Angus × Simmental cross; n = 259) born in the years 2010 to 2016 at three Auburn University Experimental Stations (Black Belt Research and Extension Center (n = 53); Gulf Coast Research and Extension Center (n = 136); Wiregrass Research and Extension Center (n = 70)). At weaning, a proportion of heifers born each year was retained as potential replacement heifers.
Heifers were managed to reach a target weight of 60% of their mature bodyweight (approximately 381 kg) by the start of their first breeding season, which began in early December of each year. Heifers at the Black Belt Research and Extension Center were weaned and developed on toxic endophyte infected tall fescue pastures and free-choice ryegrass hay. Heifers were supplemented as needed with a 50:50 mixture of corn gluten pellets and soyhull pellets. Heifers remained on tall fescue pastures from weaning through the winter grazing season and the time of pregnancy diagnosis. Heifers at the Gulf Coast Research and Extension Center were developed from weaning to breeding on bahiagrass pasture alongside free choice ryegrass hay. Heifers were supplemented as needed with a Nutrena NutreBeef® 13% protein pellet. After breeding, heifers were moved to a ryegrass pasture for the remainder of the winter grazing season. At the Wiregrass Research and Extension Center, weaned heifers were managed on bermudagrass pasture with supplementation of 50% pelleted soyhulls and 50% corn gluten feed that was provided as needed. As summer pastures entered dormancy, heifers were fed free-choice Tifton-85 bermudagrass hay and were allowed to graze pastures containing triticale, hairy vetch, and rape seed for the remainder of the winter grazing season.
Classification of heifers based on reproductive outcome
Approximately 30 d before the start of their first breeding season, heifers were evaluated for BCS (scale of 1–9 with 1 = emaciated and 9 = obese; [19, 20]) and assessed for RTS (scale of 1–5; 1 = pre-pubertal, 5 = pubertal, luteal phase; [21]) by a single, experienced veterinarian. At approximately 14 months of age (418.7 ± 22.6 days), heifers retained as replacements underwent estrous synchronization for fixed-time artificial insemination utilizing the 7-Day Co-Synch + CIDR protocol [22] to begin their first breeding season. Briefly, heifers received an injection of GnRH (i.m.; 100 μg; Cystorelin®; Merial, Duluth, GA) and insertion of a CIDR (intravaginal insert; 1.38 g progesterone; Eazi-Breed® CIDR®; Zoetis Inc., Kalamazoo, MI) on d −9, followed by CIDR removal and an injection of prostaglandin F2α (PGF; i.m.; 25 mg; Lutalyse®; Zoetis Inc., Kalamazoo, MI) on d −2. All heifers then received a second GnRH injection (i.m.; 100 μg; Cystorelin®; Merial, Duluth, GA) and were inseminated with a dose of semen of proven fertility on d 0, 54 ± 2 h after CIDR removal and PGF injection. Two professionals in random rotation were responsible for AI procedures at each experimental station for each year.
Fourteen days after insemination, heifers were exposed to two fertile bulls for natural breeding for the remainder of the breeding season. An experienced veterinarian performed initial pregnancy evaluation by transrectal palpation on d 62–89 post insemination, followed by final pregnancy evaluation on d 85–176 post insemination. Presence or absence of a conceptus, alongside morphological features indicating fetal age were recorded, and heifers were classified as “pregnant to AI’ (Preg AI), ‘pregnant to natural service’ (Preg NS), or ‘not pregnant’ (Not Preg).
Phenotypic dataset
All analytical procedures were carried out in R software [23]. A schematic diagram representing the phenotypic data, the reproductive data, the merging, and the curation procedures is depicted in (Additional file 1: Figure S1). We obtained performance records and pedigree information for all calves born at each station from 2000 to 2017 from the Alabama Beef Cattle Improvement Association (n = 2530). We then filtered this dataset to include only heifer calves (n = 1240) and merged the performance dataset with records for pre-breeding body condition score, reproductive tract score, artificial insemination date, and pregnancy outcome for all heifers exposed to breeding during their first breeding season. We computed age of dam by subtracting the year of birth for the dam from the year of birth for each heifer.
We curated the data and eliminated observations that appeared as abnormal data imputation or outliers. We retained records if weaning weight was recorded within 158.8–453.6 kg and adjusted weaning weight was less than 453.6 kg. For analyses of pregnancy outcome, we only retained records for heifers that conceived if the pregnancy was carried out to term and a healthy calf was born. Heifers experiencing pregnancy loss (n = 3) were removed from the dataset because conceptus losses were not the focus of this study and analyzing data from these heifers would create a confounding category between pregnant and not pregnant. In addition, heifers presenting RTS < 3 (n = 5) were removed from the dataset according to consistent data supporting the notion that heifers with an immature reproductive tract are significantly less likely to become pregnant [9, 16].
We assessed normality of the continuous traits by performing a Shapiro-Wilk test and by examining histograms, quantile-quantile, and density plots for each parameter. We utilized the data from all heifers to assess the normality of weaning weight, adjusted weaning weight, and age at weaning, regardless of whether we collected breeding data. We assessed normality of age at AI using the data from the heifers that were artificially inseminated. Amongst heifers included in this dataset, the variables weaning weight (WW) and adjusted weaning weight (adj WW) were normally distributed (P > 0.01, Shapiro-Wilk test, Additional file 1: Figure S2). The variables age at weaning and age at AI displayed a deviation when tested from normal distribution (P < 0.01, Shapiro-Wilk test, Additional file 1: Figure S2). Nonetheless, visual inspection of the data (Additional file 1: Figure S2) indicated strong resemblance of normal distribution and the skewness was likely an influence of the varied management strategies at each station.
Analysis of phenotypic parameters relative to pregnancy outcome
We analyzed the data using a mixed effect multinomial logistic regression model [24] because the heifers were categorized according to discrete reproductive outcomes. We modeled the phenotypic parameters relative to the reproductive outcomes according to two scenarios.
First, we accounted for the probability of three possible reproductive outcomes: Preg AI, Preg NS, or Not Preg. The variables station (Sj, j = 1, 2, 3), AI year (Yk, k = 2011, 2012, 2013, 2015, 2016, 2017), BCS (BCSl, l = 4, 5, 6, 7), RTS (RTSm, m = 3, 4, 5), age at weaning (AgeW), age at AI (AgeAI), dam age (AgeD), and weaning weight (WW) were considered for the model. Heifer’s sires, nor the bulls used in the breeding programs were included in the model as they were confounded with stations. The probabilities (Pr) of occurrence of each pregnancy outcome relative to the variables were estimated as follows:
$$ \mathit{\ln}\left(\frac{\Pr (PregAI)}{\Pr (NotPreg)}\right)={\beta}_{01}+{\beta}_{11}{S}_j+{\beta}_{21}{Y}_k+{\beta}_{31}{BCS}_l+{\beta}_{41}{RTS}_m+{\beta}_{51} AgeAI+{\beta}_{61} AgeD+{\beta}_{71} WW+{\varepsilon}_1 $$
(and)
$$ \mathit{\ln}\left(\frac{\Pr (PregNS)}{\Pr (NotPreg)}\right)={\beta}_{02}+{\beta}_{12}{S}_j+{\beta}_{22}{Y}_k+{\beta}_{31}{BCS}_l+{\beta}_{42}{RTS}_m+{\beta}_{52} AgeAI+{\beta}_{62} AgeD+{\beta}_{72} WW+{\varepsilon}_2 $$
Next, we accounted for the probability of two outcomes only: ‘pregnant’ or ‘not pregnant’. The binomial modeling followed the same structure as presented above with exception that the dependent variable was represented by \( \mathit{\ln}\left(\frac{\Pr (Preg)}{\Pr \left( Not\ preg\right)}\right) \).
We used the ‘nnet’ package [25] to fit the multinomial and binomial models. The likelihood of the ratios was calculated with a χ2 test using the ‘Anova’ function in the ‘car’ package. The model was assessed by the Akaike Information Criteria (AIC) [26] using the ‘MASS’ package. Statistical significance was inferred if P < 0.05.
The codes utilized for the analyses presented on this paper can be found on Additional file 2.