Data
All data originated from purebred cattle registered in the American Angus Association (AAA) and commercial cattle enrolled in the AAA Breed Improvement Record program. Phenotypic data comprised hair shedding scores recorded by beef cattle producers enrolled in the Mizzou Hair Shedding Project (MU data) between 2016 and 2019 in combination with hair shedding scores collected by technicians in 2011, 2012, 2018, and 2019 as part of Angus Foundation-funded projects at Mississippi State University and North Carolina State University (AGI data). Across all years and both datasets, scores were recorded on 1 day between April 17th and June 30th in the late spring or early summer, with most scores recorded in mid- to late-May. Hair shedding was evaluated using a 1 to 5 visual appraisal scale, where 1 was 0% dead winter coat remaining and 5 was 100% winter coat remaining based on the systems developed by Turner and Schleger [21] and Gray et al. [22] (Fig. 1). While there is variation in the onset of hair shedding between individuals, cattle and other mammals tend to shed from the head towards the tail and from the topline towards the legs [2, 8, 23].
Records were removed when the breeder-reported sex of an animal did not match the sex recorded in the AAA pedigree. Hair shedding scores that originated from male animals comprised less than 1% of the dataset and only female records were retained. Age classifications were assigned to each record based on age in days determined by the AAA-recorded birth date and the date the hair shedding score was recorded. Similar to the system used in the Beef Improvement Federation (BIF) Guidelines for age-of-dam classification [24], age classifications were defined as \(\left( {n*365d} \right) - 90d\) to \(\left( {\left( {n + 1} \right)*365d} \right) - 90d\), where \(n\) is the age classification and \(d\) is days. Records where the breeder-reported age in years differed from the calculated age classification by more than 2 years and records from animals younger than 275 days-of-age were removed. When no calving season was reported, it was imputed using the most recent natural birth calving date available in the AAA database prior to the recorded score. When no natural birth calving dates were available, calving season was imputed using the animal’s own birth date. In the AGI data, some animals were scored by multiple scoring technicians on the same day. In these cases, phenotypes recorded on the same animal and the same day were averaged. In the MU data, participating producers were asked to report whether or not (yes or no) animals grazed toxic fescue during the spring of the recording year. Grazing status was not explicitly recorded in the AGI data, but animals scored in Texas were assumed not to have grazed toxic fescue. This resulted in 14,839 scores in the combined, filtered dataset. Among the 8619 individuals included, 49% had between 2 and 6 years of data. Most data came from herds in the Southeast and Fescue Belt (Fig. 2). The mean hair shedding score was slightly higher in the AGI data (\(\mu\) = 3.10; \(n\) = 6374) compared to the MU data (\(\mu\) = 2.65; \(n\) = 8465), but the standard deviation was identical in both datasets (\(\sigma\) = 1.15).
Genotypes and imputation
Genotypes for 3898 of the 8619 animals were imputed to a union marker set (\({\text{n}}\) = 233,246) of the GGP-F250 genotyping chip and various commercial assays using FImpute v.3.0 [25]. The commercial assays were those used in routine genotyping of Angus cattle for genomic selection purposes, which include ~ 50,000 markers or a lower density panel that can be imputed to ~ 50,000 with sufficient accuracy. Although FImpute provides the capacity to infer the genotypes of un-genotyped animals based on information from relatives, markers were imputed only for genotyped individuals. Prior to imputation, markers with a GenCall score lower than 0.15 were set to missing and individuals with Mendelian error rates higher than 2% had their parents set to missing in the pedigree. The GGP-F250 was designed to genotype functional variants and thus has more variants at low minor allele frequencies [26]. Therefore, no minor allele frequency filter was applied during or after imputation beyond the removal of monomorphic SNPs. Animals and markers with call rates lower than 85% were removed. The resulting marker set consisted of 174,894 autosomal variants.
Construction of the blended relationship matrix \({\mathbf{H}}^{ - 1}\)
In single-step genomic best linear unbiased prediction (BLUP) as used in the AAA National Cattle Evaluation (NCE), relationships between individuals are represented in the matrix \({\mathbf{H}}^{ - 1}\), which is a blended form of the genomic and additive relationship matrices [27], allowing information from both genotyped and non-genotyped animals to be used. \({\mathbf{H}}^{ - 1}\) is calculated as:
$${\mathbf{A}}^{ - 1} + \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & {{\mathbf{G}}_{\varvec{w}}^{ - 1} - {\mathbf{A}}_{22}^{ - 1} } \\ \end{array} } \right],$$
where \({\mathbf{A}}^{ - 1}\) represents the inverted pedigree relationship matrix traditionally used to represent relationships, \({\mathbf{A}}_{22}^{ - 1}\) represents the inverted pedigree relationship matrix for the subset of animals with genotypes available, and \({\mathbf{G}}_{\varvec{w}}^{ - 1}\) is the inverted genomic relationship matrix. The genomic relationship matrix was calculated using the VanRaden method [28] and was blended with \({\mathbf{A}}_{22}\) with the default weight of 0.05 using the preGSf90 program [29]. In all subsequent models including a random genetic effect, \({\mathbf{H}}^{ - 1}\) was constructed using the 3-generation pedigree (in total, 17,652 animals; 1987 distinct sires and 9509 distinct dams) in combination with imputed genotypes.
Effect of age on hair shedding score and contemporary group definition
Understanding how and which environmental factors shape phenotypic variation enables the development of more appropriate contemporary group definitions during genetic evaluation. In order to quantify the effect of animal age on hair shedding score, we fitted age as a categorical fixed effect in a repeated records animal model. Age categories were defined in three ways. First, age in years was fit (i.e. all possible values between 1 and 16). Second, ages were grouped as 1, 2, 3, or other (“four-class model”). Third, age groups were defined according to the guidelines set forth by the BIF for age-of-dam effects on birth weight and weaning weight (i.e., 2, 3, 4, 5-9, 10, 11, 12, 13+; [24]) plus yearlings (“BIF model”). The four-class model and the BIF model were each compared against a null model with no age effect included using Akaike’s Information Criterion (AIC) and likelihood ratio tests. In all three models with age classification fitted as a categorical fixed effect, classifications with fewer than five animals were excluded. These models are summarized below:
$${\mathbf{y}} = {\mathbf{X}}_{1} {\mathbf{b}} + {\mathbf{X}}_{2} {\mathbf{a}} + {\mathbf{Z}}_{1} {\mathbf{u}} + {\mathbf{Z}}_{2} {\mathbf{p}} + {\mathbf{e}},$$
where \({\mathbf{y}}\) is a vector of hair shedding scores; \({\mathbf{b}}\) is a vector of contemporary group effects for each hair shedding score, with contemporary group defined as farm ID, year scored, calving season, score group, and toxic fescue grazing status; \({\mathbf{a}}\) is a vector of age classification effects for each individual (based on age-in-years, BIF classifications, or the four age classes); \({\mathbf{u}}\) is the random additive genetic effect with \({\mathbf{u}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{H}}\sigma_{\text{a}}^{2} } \right)\); \({\mathbf{p}}\) is the random permanent environment effect with \({\mathbf{p}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{pe}}^{2} } \right)\); \({\mathbf{e}}\) is the random residual with \({\mathbf{e}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{e}}^{2} } \right)\); and \({\mathbf{X}}_{1}\), \({\mathbf{X}}_{2}\), \({\mathbf{Z}}_{1}\), and \({\mathbf{Z}}_{2}\) are incidence matrices relating the elements of \({\mathbf{y}}\) to \({\mathbf{b}}\), \({\mathbf{a}}\), \({\mathbf{u}}\), and \({\mathbf{p}}\), respectively.
Effect of toxic fescue grazing status on hair shedding
Cattle reared in heat-stressed regions but not exposed to endophyte-infected fescue demonstrate similar benefits from early summer hair shedding, but it is unclear if the biological mechanisms that govern hair shedding under fescue toxicosis and heat stress alone are the same. This could have implications for routine genetic evaluation, as it might require that some hair shedding score observations be treated as a separate trait. In order to clarify the relationship between hair shedding score while grazing toxic fescue versus while not grazing toxic fescue, we calculated the covariance and genetic correlation between hair shedding score grazing toxic fescue and not grazing toxic fescue using the bivariate repeated records animal model below:
$${\mathbf{y}}_{{\mathbf{t}}} = {\mathbf{X}}_{{\mathbf{t}}} {\mathbf{b}}_{{\mathbf{t}}} + {\mathbf{Z}}_{{1{\mathbf{t}}}} {\mathbf{u}}_{{\mathbf{t}}} + {\mathbf{Z}}_{{2{\mathbf{t}}}} {\mathbf{p}}_{{\mathbf{t}}} + {\mathbf{e}}_{{\mathbf{t}}} ,$$
where \({\mathbf{y}}\) is a vector hair shedding scores and \({\mathbf{t}}\) is toxic fescue grazing status (yes or no); \({\mathbf{b}}\) is a vector of contemporary group effects for each hair shedding score, with contemporary groups defined as farm ID, year scored, calving season, score group, and age class (yearling, 2-year-old, 3-year-old, or other; based on the results of the age classification analyses above); \({\mathbf{u}}\) is the additive genetic effect and \({\text{Var}}\left( {\mathbf{u}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{uYes}}^{2} } & {\sigma_{{{\text{uYes}},{\text{uNo}}}} } \\ {\sigma_{{{\text{uNo}},{\text{uYes}}}} } & {\sigma_{\text{uNo}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{H}}\); \({\mathbf{p}}\) is the permanent environment effect and \({\text{Var}}\left( {\mathbf{p}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{pYes}}^{2} } & 0 \\ 0 & {\sigma_{\text{pNo}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{I}}\); \({\mathbf{e}}\) is the random residual and \({\text{Var}}\left( {\mathbf{e}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{eYes}}^{2} } & {\sigma_{{{\text{eYes}},{\text{eNo}}}} } \\ {\sigma_{{{\text{eNo}},{\text{eYes}}}} } & {\sigma_{\text{eNo}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{I}}\); and \({\mathbf{X}}\), \({\mathbf{Z}}_{1}\), and \({\mathbf{Z}}_{2}\) are incidence matrices relating the elements of \({\mathbf{y}}\) to \({\mathbf{b}}\), \({\mathbf{u}}\), and \({\mathbf{p}}\), respectively.
In addition, we fitted a univariate model with toxic fescue grazing status included as a categorical fixed effect. The goal of this model was to quantify the effect of reported toxic fescue grazing status on hair shedding score:
$${\mathbf{y}} = {\mathbf{X}}_{1} {\mathbf{b}} + {\mathbf{X}}_{2} {\mathbf{f}} + {\mathbf{Z}}_{1} {\mathbf{u}} + {\mathbf{Z}}_{2} {\mathbf{p}} + {\mathbf{e}},$$
where \({\mathbf{y}}\) is a vector of hair shedding scores; \({\mathbf{b}}\) is a vector of contemporary group effects, defined in the same way as the univariate model above; \({\mathbf{f}}\) is the toxic fescue status effect (yes or no); \({\mathbf{u}}\) is the additive genetic effect with \({\mathbf{u}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{H}}\sigma_{\text{a}}^{2} } \right)\); \({\mathbf{p}}\) is the permanent environment effect with \({\mathbf{p}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{pe}}^{2} } \right)\); \({\mathbf{e}}\) is the random residual with \({\mathbf{e}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{e}}^{2} } \right)\); and \({\mathbf{X}}_{1}\), \({\mathbf{X}}_{2}\), \({\mathbf{Z}}_{1}\), and \({\mathbf{Z}}_{2}\) are incidence matrices relating the elements of \({\mathbf{y}}\) to \({\mathbf{b}}\), \({\mathbf{f}}\), \({\mathbf{u}}\), and \({\mathbf{p}}\), respectively.
In both models, only females with known toxic fescue grazing status were retained for analysis. Contemporary groups with fewer than five animals or no variation were discarded, resulting in 5832 observations from cattle grazing toxic fescue and 4197 observations from cattle not grazing toxic fescue. Three hundred ninety-six animals had observations over multiple years both grazing and not grazing toxic fescue.
Genetic parameters, breeding values, and estimated bias
Variance components, heritability, repeatability, and breeding values were estimated using the univariate repeated records animal model below implemented in AIREMLF90 [29].
$${\mathbf{y}} = {\mathbf{Xb}} + {\mathbf{Z}}_{1} {\mathbf{u}} + {\mathbf{Z}}_{2} {\mathbf{p}} + {\mathbf{e}},$$
where \({\mathbf{y}}\) is a vector of hair shedding scores; \({\mathbf{b}}\) is the contemporary group effect; \({\mathbf{u}}\) is the additive genetic effect with \({\mathbf{u}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{H}}\sigma_{\text{a}}^{2} } \right)\); \({\mathbf{p}}\) is the permanent environment effect with \({\mathbf{p}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{pe}}^{2} } \right)\); \({\mathbf{e}}\) is the random residual with \({\mathbf{e}} \sim {\text{N}}\left( {{\mathbf{0}},{\mathbf{I}}\sigma_{\text{e}}^{2} } \right)\); and \({\mathbf{X}}\), \({\mathbf{Z}}_{1}\), and \({\mathbf{Z}}_{2}\) are incidence matrices relating the elements of \({\mathbf{y}}\) to \({\mathbf{b}}\), \({\mathbf{u}}\), and \({\mathbf{p}}\), respectively.
The definition of contemporary groups used in this final prediction model was informed by the results of the age classification and toxic fescue grazing status analyses above. It included a combination of farm, year scored, calving season (spring or fall), toxic fescue grazing status (yes or no), age group (yearling, 2-year-old, 3-year-old, or other), and score group. In herds where cattle were scored for hair shedding over more than 1 day, the score group was determined using a 7-day sliding window to maximize the number of animals per contemporary group. In the future, it will be recommended that producers score all cattle for hair shedding within a week of one another to maximize the size of contemporary groups. Although yearling heifers have not yet experienced the stress of pregnancy, calving season/birth season is a good proxy for management group in the absence of breeder-reported codes. Therefore, “calving season” was included in the contemporary group definition for all animals regardless of reproductive status. Contemporary groups with fewer than five animals or no variation were dropped. This resulted in 14,438 total scores from 8449 animals in 395 contemporary groups.
In order to evaluate model bias, we estimated breeding values in ten separate iterations, excluding all phenotypes for a randomly selected 25% of animals. These “partial” breeding values were then compared to breeding values obtained via the “whole” model including all possible information using the “LR method” parameters proposed by Legarra and Reverter [30]. First, we calculated the absolute difference between whole breeding values and partial breeding values for the validation set, or animals whose phenotypes were excluded (\(d_{w,p}^{v}\)) and the reference set, or animals whose phenotypes were not excluded (\(d_{w,p}^{r}\)). The expectation of this value is zero in the absence of bias, where bias is introduced by incorrect estimation of the genetic trend. Next, we regressed whole breeding values on partial breeding values for both validation (\(b_{w,p}^{v}\)) and reference (\(b_{w,p}^{r}\)) sets. In this model, deviations of the slope from 1 are suggestive of dispersion. Finally, we calculated the correlation between partial and whole breeding values (\(\rho_{p,w} = \frac{{cov\left( {\hat{\mu }_{p} ,\hat{\mu }_{w} } \right)}}{{ \sqrt{var(\hat{\mu }_{p}) var (\hat{\mu }_{w})} }}\)) within the validation and reference sets, where the correlation within the validation set (\(\rho_{p,w}^{v}\)) is a metric of prediction accuracy.
Weaning weight
The effects of heat stress on pre-weaning growth are well characterized in cattle. Heat stress impacts calf performance most severely via reduced milk production in the dam [12]. Fescue toxicosis induces reduced milk production in a similar fashion [31]. Therefore, we quantified the phenotypic and genetic correlations between hair shedding score and weaning weight.
Weaning weight phenotypes and contemporary group designations came from the weekly growth run of the AAA national cattle evaluation (NCE). Prior to entering the NCE, phenotypes were adjusted for age-of-dam effects as used in the AAA weekly NCE and to 205 days-of-age. Weaning weight data were retrieved for: (1) own weaning weight of cows with at least one hair shedding score recorded, (2) all of cow’s recorded calves, (3) cow’s weaning weight contemporary group peers, and (4) all of their recorded calves’ weaning weight contemporary group peers. Weaning weights from animals born via embryo transfer and contemporary groups with fewer than five animals or no variation were excluded, resulting in 40,794 total weaning weights and 14,039 total hair shedding score records. Of the 45,420 phenotyped animals retained for analysis, 3850 had both a recorded weaning weight and at least one hair shedding score. Furthermore, 6448 dams had both hair shedding scores and calf weights recorded in at least 1 year (n = 9092 score/weight pairs).
Conceivably, environmental factors that affect a dam’s hair shedding performance could also affect the direct weaning weight of her calf and her maternal effect on the calf’s growth, creating a residual covariance between the two traits. In order to reflect this covariance, a bivariate model was fitted in which a direct hair shedding score effect was modeled for the cow, a direct weaning weight effect was modeled for the calf, and a maternal weaning weight effect was modeled for the cow. In practice, this model was implemented by fitting a maternal genetic effect for hair shedding, no direct genetic effect of hair shedding (no genetic effect of the calf on the hair shedding score of its dam), and direct and maternal genetic effects for weaning weight. This model created a direct tie between a dam’s hair shedding score and the calf she weaned that year, which reflects more accurately the relationship of interest and is similar to models used to assess the correlations between weaning weight and actual milk yield [32]. For cows with a hair shedding score but no calf weaning weight reported during the scoring year, a “dummy calf” with a weaning weight set to missing and unknown sire was created. This model was fitted three separate times: once including only dams explicitly recorded to have been grazing toxic fescue, once including only dams explicitly recorded to have not been grazing toxic fescue, and once with all available data.
$$\left[ {\begin{array}{*{20}c} {{\mathbf{y}}_{{{\mathbf{HS}}}} } \\ {{\mathbf{y}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {{\mathbf{X}}_{{{\mathbf{HS}}}} } & 0 \\ 0 & {{\mathbf{X}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\mathbf{b}}_{{{\mathbf{HS}}}} } \\ {{\mathbf{b}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & {{\mathbf{Z}}_{{1_{{{\mathbf{WW}}}} }} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} 0 \\ {{\mathbf{u}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right]$$
$$+ \left[ {\begin{array}{*{20}c} {{\mathbf{Z}}_{{2_{{{\mathbf{HS}}}} }} } & 0 \\ 0 & {{\mathbf{Z}}_{{2_{{{\mathbf{WW}}}} }} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\mathbf{m}}_{{{\mathbf{HS}}}} } \\ {{\mathbf{m}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {{\mathbf{Z}}_{{3_{{{\mathbf{HS}}}} }} } & 0 \\ 0 & {{\mathbf{Z}}_{{3_{{{\mathbf{WW}}}} }} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {{\mathbf{mpe}}_{{{\mathbf{HS}}}} } \\ {{\mathbf{mpe}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {{\mathbf{e}}_{{{\mathbf{HS}}}} } \\ {{\mathbf{e}}_{{{\mathbf{WW}}}} } \\ \end{array} } \right],$$
where \({\mathbf{y}}_{{\mathbf{t}}}\) is the phenotype and \({\mathbf{t}}\) is the trait (hair shedding score (HS) or weaning weight (WW)); \({\mathbf{b}}_{{\mathbf{t}}}\) is the contemporary group effect; \({\mathbf{u}}_{{\mathbf{t}}}\) is the calf genetic effect (fit only for weaning weight) and \({\text{Var}}\left( {\mathbf{u}} \right) = \left[ {\begin{array}{*{20}c} 0 & 0 \\ 0 & {\sigma_{\text{uWW}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{H}}\), where \(\sigma_{\text{uWW}}^{2}\) represents the genetic variance for the calf direct effect of weaning weight; \({\mathbf{m}}_{{\mathbf{t}}}\) is the cow genetic effect and \({\text{Var}}\left( {\mathbf{m}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{mHS}}^{2} } & {\sigma_{{{\text{mHS}},{\text{mWW}}}} } \\ {\sigma_{{{\text{mWW}},{\text{mHS}}}} } & {\sigma_{\text{mWW}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{H}}\), where \(\sigma_{\text{mHS}}^{2}\) represents the genetic variance for hair shedding and \(\sigma_{\text{mWW}}^{2}\) represents the genetic variance for the cow maternal effect of weaning weight; \({\mathbf{mpe}}_{{\mathbf{t}}}\) is the cow permanent environment effect and \({\text{Var}}\left( {{\mathbf{mpe}}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{mpeHS}}^{2} } & {\sigma_{{{\text{mpeHS}},{\text{mpeWW}}}} } \\ {\sigma_{{{\text{mpeWW}},{\text{mpeHS}}}} } & {\sigma_{\text{mpeWW}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{I}}\), where \(\sigma_{\text{mpeHS}}^{2}\) represents the permanent environmental variance for hair shedding and \(\sigma_{\text{mpeWW}}^{2}\) represents the permanent environmental variance for the maternal effect of weaning weight; \({\mathbf{e}}_{{\mathbf{t}}}\) is the random residual and \({\text{Var}}\left( {\mathbf{e}} \right) = \left[ {\begin{array}{*{20}c} {\sigma_{\text{eHS}}^{2} } & {\sigma_{{{\text{eHS}},{\text{eWW}}}} } \\ {\sigma_{{{\text{eWW}},{\text{eHS}}}} } & {\sigma_{\text{eWW}}^{2} } \\ \end{array} } \right] \otimes {\mathbf{I}}\); and \({\mathbf{X}}\), \({\mathbf{Z}}_{1}\), \({\mathbf{Z}}_{2}\), and \({\mathbf{Z}}_{3}\) are incidence matrices relating the elements of \({\mathbf{y}}\) to \({\mathbf{b}}\), \({\mathbf{u}}\), \({\mathbf{m}}\), and \({\mathbf{mpe}}\), respectively.
We also evaluated the phenotypic relationship between dam hair shedding score and calf weaning weight using the subset of 6448 dams with both hair shedding scores and calf weights recorded in at least 1 year. We did this by calculating the estimated change in calf weaning weight as a function of dam hair shedding score using four separate simple linear regression models. In the first two models, unadjusted calf weaning weight was regressed on unadjusted dam hair shedding score. Using weaning weight unadjusted for age in days captures increased gain from an earlier birth date (older when weighed), which might be an indicator of increased fertility for earlier shedding cows. In the other two models, 205-day, age-of-dam, and contemporary group solution adjusted calf weaning weight was regressed on unadjusted dam hair shedding score. Both the unadjusted weaning weight and adjusted weaning models were fitted separately for all available data, dams explicitly recorded as grazing toxic fescue, and dams explicitly recorded as not grazing toxic fescue.
Genome-wide association
In order to evaluate the genetic architecture of hair shedding and identify variants that contribute to hair shedding score breeding values, we performed a single-SNP genome-wide association analysis using the SNP1101 v.1 software [33]. The breeding values calculated above using AIREMLF90 were de-regressed and used as the phenotype such that each of the 3783 animals had one record. The de-regressed breeding values were weighted by their reliability \(1 - \frac{PEV}{{\sigma_{a}^{2} }}\), where \(PEV = \left( {SE^{2} } \right)*\sigma_{e}^{2}\) and \(\sigma_{a}^{2}\) and \(\sigma_{e}^{2}\) are the estimated additive genetic and residual variances for hair shedding score, respectively. Heritability was constrained to 0.40 and the genomic relatedness matrix used to control for family structure was calculated using the VanRaden method [28].
Using the UMD 3.1 bovine genome assembly [34] coordinates and annotations, we searched genes within 50 kb of SNPs with a genome-wide q-value lower than 0.05. The size of our search space was determined based on the density of our marker set, and the resulting gene list was used as input for cluster enrichment analysis within ClueGO v.2.5.6 [35]. KEGG pathways and biological process gene ontologies with at least four associated genes were considered for search terms. We also searched for protein–protein interaction between genes in our list using STRING v.10 [36], considering co-expression, experimental data, and curated databases as active interaction sources.