Gabriel A. Zayas1*
Camila Santos Rojas1
Eduardo E. Rodriguez1Aakilah S. Hernandez2Ashley M. Beard1
Fahad Rafiq1
Kaitlyn M. Sarlo Davila3
Raluca G. Mateescu1- 1Department of Animal Sciences, University of Florida, Gainesville, FL, United States
- 2Department of Animal Science, North Carolina State University, Raleigh, NC, United States
- 3Virus and Prion Research Unit, National Animal Disease Center, Agricultural Research Service, United States Department of Agriculture, Ames, IA, United States
Background: Rising temperatures increasingly expose beef cattle to heat stress, reducing productivity and welfare, especially in tropical climates. Crossbreeding Bos t. taurus and Bos t. indicus has emerged as a critical strategy to balance the production efficiency of taurine breeds with the superior thermotolerance of indicine breeds. Understanding the genetic architecture of thermotolerance traits is essential for improving heat resilience in beef cattle populations.
Methods: Phenotypes for short hair length (SHL, undercoat) and long hair length (LHL, topcoat), sweat gland area (SGA), and thermal stress slope (TSS), a measure of body temperature fluctuations under heat stress, were collected from 3,962 crossbred Angus-Brahman heifers. Heifers were genotyped, and breed-of-origin (BOA) for each marker was determined using LAMP-LD. Genome-wide association studies were conducted using SNP-only, BOA-only, and integrated SNP + BOA models to identify quantitative trait loci (QTLs) associated with thermotolerance traits. Genes in QTL regions were used for functional enrichment analysis using Gene Ontology (GO) and KEGG pathways.
Results: Significant QTLs for SHL and LHL were identified on BTA20, overlapping the PRLR gene. A QTL on BTA19 for SHL and LHL was driven solely by BOA effects, with Brahman BOA associated with shorter hair lengths. For SGA, six suggestive QTLs were detected, predominantly linked to Angus-derived alleles associated with reduced sweat gland area. For TSS, a significant QTL on BTA1 exhibited a strong BOA effect, with Angus BOA associated with higher TSS values, indicative of reduced thermoregulatory efficiency. Integrated SNP + BOA models provided greater resolution and revealed novel QTLs compared to single-effect models. Functional enrichment using GO and KEGG identified MAPK and estrogen signaling pathways in both LHL and TSS, indicating potential overlap in the biological processes influencing hair length and thermoregulation.
Conclusion: This study demonstrates the value of integrating BOA with SNP-based models to uncover the genetic architecture of thermotolerance traits in beef cattle. By better capturing breed-specific contributions, these findings enhance our understanding of thermoregulation and provide actionable insights for improving heat resilience in cattle.
1 Background
Rising global temperatures have introduced challenges for beef production, as approximately 80% of cattle globally experience heat stress conditions for at least 30 days per year, with the severity being greater in tropical regions (Carvajal et al., 2021). Heat stress adversely affects productivity, health, and welfare in cattle, necessitating strategies to enhance thermotolerance. Thermotolerance refers to the cattle’s ability to regulate body temperatures whilst maintaining production in the presence of heat stress conditions. One promising approach to improving thermotolerance in cattle is crossbreeding Bos t. taurus and Bos t. indicus breeds. Bos t. taurus breeds are valued for their production efficiency, while Bos t. indicus breeds are generally very heat tolerant. By combining these attributes, crossbreeding offers a practical strategy to improve both productivity and thermotolerance in beef cattle. The evolutionary history of Bos t. indicus and Bos t. taurus cattle are quite different, as each originated from distinct domestication events and experienced different selection pressures and adaptations (Bolormaa et al., 2013; Koufariotis et al., 2018). Bos t. taurus breeds have been selectively bred for a variety of traits that optimize production, leading to the development of numerous specialized breeds, each selected for specific purposes, such as meat quality, growth rate, or milk production. In contrast, Bos t. indicus cattle have undergone selection for resilience in tropical climates, with distinct physical adaptations for regulating body temperature in high heat and humidity (Utsunomiya et al., 2019; Verdugo et al., 2019; Mateescu et al., 2020). These adaptations include a larger skin surface area, longer ears (Utsunomiya et al., 2019), slicker hair coats (Sarlo Davila et al., 2019), and larger sweat glands for efficient heat dissipation (Mateescu et al., 2023; Hernandez et al., 2024). The observable morphological differences between Bos t. taurus and Bos t. indicus cattle are a direct reflection of underlying variations in their genetic architecture (Koufariotis et al., 2018). Crossing breeds from these two different origins can result in cattle that truly embody the “best of both worlds.”
Understanding the genetic basis of thermotolerance in beef cattle is essential for improving performance in hot and humid environments. The mixing of Bos t. taurus and Bos t. indicus genetic backgrounds often results in greater variation in thermotolerance traits, creating opportunities to identify and select individuals with superior heat tolerance. Previous studies have described the heritability of hair length traits (Sarlo Davila et al., 2019), sweat gland area (Mateescu et al., 2023; Hernandez et al., 2024) and core body temperature traits (Sarlo Davila et al., 2019) in crossbred Angus-Brahman cattle, reporting high to moderate heritabilities. This variability enhances the potential for identifying genetic regions influencing thermotolerance, as greater genetic and phenotypic variation increases the power to detect associations between genetic markers and traits of interest. However, differences in the evolutionary histories of Bos t. taurus and Bos t. indicus breeds result in differences in linkage disequilibrium (LD) patterns, which present challenges for genomic analyses. In crossbred and composite populations, the LD structure is often less stable, as it represents a mixture of the parental breeds’ LD patterns. This can reduce the accuracy of identifying quantitative trait loci (QTL), as these analyses rely on LD between genetic markers and causal variants. Addressing these challenges requires advanced genomic approaches that take into account crossbred populations unique genetic structure.
Breed of Origin of Alleles (BOA) offers an innovative solution by capturing breed-specific genomic contributions (Vandenplas et al., 2016; Claudia et al., 2019). BOA assigns each genetic marker to its ancestral breed, enabling simultaneous modeling of SNP and BOA effects. This approach enhances the precision of genomic evaluations by accounting for the distinct genetic architectures of the contributing breeds, thereby facilitating more accurate predictions of traits (Guillenea et al., 2023). By accurately assigning each allele to its breed of origin, BOA enables researchers to discern how breed-specific genetic contributions influence various traits in crossbred animals. This joint approach ensures that QTL identification accounts for both genomic variation and breed-specific contributions to trait expression. Previous studies in crossbred Angus-Brahman cattle have identified BOA effects associated with economically relevant traits such as marbling and carcass weight (Zayas et al., 2024b). However, these studies modeled SNP and BOA effects separately, which can limit the ability to detect meaningful genomic associations. Fitting both SNP and BOA information simultaneously, can improve the robustness and accuracy of QTL detection in crossbred cattle.
This study aimed to identify QTLs associated with thermotolerance traits in beef cattle, including short hair length (SHL), long hair length (LHL), sweat gland area (SGA) and thermal stress slope (TSS). Genome-wide association studies (GWAS) were conducted modeling SNP and BOA effects simultaneously. This approach was compared to traditional SNP-only and BOA-only models to identify genomic regions essential for thermotolerance traits.
2 Methods
2.1 Animals and phenotypic measurements
The research protocol for this study was approved by the University of Florida Institutional Animal Care and Use Committee (IACUC #202200000019). The study population consisted of 3,962 commercial Brangus heifers and Angus-Brahman crossbred heifers from the University of Florida’s multibreed Angus-Brahman (MAB) herd. A total of 2,772 commercial Brangus heifers were sampled from the Seminole Tribe of Florida, Inc., located in Okeechobee, FL, over the years 2016–2019. Additionally, 669 commercial Brangus heifers were sourced from Williamson Cattle Company in Chiefland, FL, during 2021 and 2022. The UF MAB heifers comprised 521 individuals collected in 2017, 2018, 2023, and 2024.
Phenotypic measurements and DNA samples, collected via tissue or blood, were obtained during the summer months. Heifers were sampled in groups of 70–200, with a total of 22 groups. All heifers within a given group and year were from the same cohort and were of similar ages, reflecting uniform management practices. Sampling occurred during the summer months to capture phenotypic data relevant to heat stress under natural environmental conditions. Four traits were measured to assess thermotolerance traits: SHL, LHL, SGA and TSS.
2.1.1 Short and long hair length
Hair samples were collected from the shoulder, 10 cm below the spine and halfway along the horizontal axis of each animal, as described in previous studies (Sarlo Davila et al., 2019; Davila et al., 2020). Samples were spread on white paper and photographed alongside a ruler, which was used to convert pixels to millimeters. Hair lengths were measured using ImageJ software (Schindelin et al., 2012). The undercoat (SHL) and topcoat (LHL) were assessed by averaging five short hairs and five long hairs from each animal. Outliers, defined as hairs with lengths exceeding three standard deviations from the mean, were excluded to ensure accurate representation of typical hair length.
2.1.2 Sweat gland area
Skin biopsies were collected from the shoulder, 10 cm below the spine and midway along the horizontal axis of each animal, following previously described methods (Mateescu et al., 2023; Hernandez et al., 2024). Briefly, the skin was disinfected with 70% ethanol and a 2% chlorhexidine solution, then sprayed with 4% lidocaine topical anesthetic. A 6-mm biopsy punch (Miltex Inc., York, PA, USA) was used to obtain the sample, which was then fixed in 10% formalin for 16–24 h at room temperature. After fixation, biopsies were dehydrated in 70% ethanol, embedded in paraffin, and sectioned at a thickness of 7 μm. Sections were mounted on slides, stained with Harris-Eosin Hematoxylin, and digitized using a Nikon T3000 inverted phase microscope equipped with NIS Image Elements software. Images were captured at ×40 magnification and analyzed using ImageJ software (Schindelin et al., 2012). SGA was measured as the total area occupied by sweat glands within a 1,100 × 1,100-pixel image section. Measurements were converted from pixels to millimeters using a scale of 1,000 pixels = 2.145 mm.
2.1.3 Thermal stress slope
Core body temperature was measured as vaginal temperature using iButton data loggers (type DS1922L, temperature range −40°C to 85°C, accuracy ±0.5°C, resolution ±0.0625°C; Maxim Integrated, San Jose, CA) during a 5-day period during the summer months across each year of collection. iButtons were attached to blank controlled internal drug-release (CIDR) devices and inserted into the vagina of each heifer. Before deployment, all iButtons were calibrated to ensure measurement accuracy and preprogrammed to record body temperature at 15-min intervals on a 24-h cycle. Data from the first and last day, corresponding to CIDR insertion and removal, were excluded.
Environmental conditions, including temperature, humidity, dew-point temperature, and black globe temperature, were continuously monitored using HOBO data loggers (Onset Computer Corp., Bourne, MA). Hourly averages of body temperature were matched with the corresponding hourly averages of THI recorded during the same time. The temperature-humidity index (THI) was calculated using the formula:
where T is air temperature (°C) and RH is relative humidity (%). This equation has been validated as an indicator of heat stress in livestock (Animals, 1971; Dikmen and Hansen, 2009).
To quantify the impact of daily thermal fluctuations on beef cattle under heat stress conditions a novel metric was used. TSS was calculated for each individual as the ratio of the diurnal variation in body temperature—determined by the difference between maximum and minimum body temperatures within a 24-h period—to the THI load for the same day. The THI load is defined as the cumulative deviation of the THI above a baseline threshold of 70, reflecting the intensity and duration of thermal stress experienced by the animals. To ensure consistency in assessing heat stress, a single day was selected for each group of animals. These days were chosen based on high THI load and the absence of rain events, which could otherwise confound the effects of thermal stress. The following formula was used to calculate TSS:
where VTmax and VTmin represent the maximum and minimum vaginal temperatures recorded during the day, and THIi represents the hourly THI. This approach quantifies an individual’s thermoregulatory capacity under naturally occurring environmental heat stress, offering a precise metric for genetic and phenotypic analysis.
2.2 Genotyping and quality control
Genomic DNA was extracted from either tissue or blood samples using the QIAamp DNA Mini DNA kit (Qiagen, Valencia, CA, United States) following the manufacturer’s protocol and stored at −20°C. Genotyping was conducted using the Bovine GGP F250 array (GeneSeek, Inc., Lincoln, NE) containing 221,115 SNPs enriched with functional variants including non-synonymous, frameshift, and stop mutations. Only autosomal SNP were mapped to the ARS-UCD1.2 assembly and retained for further analysis. Quality control (QC) filtering was performed with PLINK2 (Chang et al., 2015). QC filters excluded animals with a genotype completion rate below 90%, markers with a minor allele frequency below 5% and a genotype call rate below 99%.
2.3 Breed of origin assignment
Local Ancestry in adMixed Populations (LAMP-LD) was used to infer percentages of local ancestry of each animal (Paşaniuc et al., 2009; Baran et al., 2012). LAMP-LD uses hidden Markov models of haplotype diversity of the ancestral/purebred populations within a window-based framework to trace the origin of alleles in the admixed population (Paşaniuc et al., 2009; Baran et al., 2012). Purebred Angus and Brahman cattle from the University of Florida’s MAB herd were used to represent the purebred populations for the LAMP-LD Analysis. A total of 123 purebred Angus cattle and 406 purebred Brahman cattle were used as the reference population. Only markers with an allele frequency difference (AFD) ≥5% between purebred population were used, to ensure sufficient differentiation between breeds. After QC filtering and AFD filtering a total of 107,153 markers remained for further analysis. The local ancestry results from LAMP-LD were subsequently used to infer the BOA. The BOA of the resulting 107,153 SNPs were then converted into a pseudo-genotype format using in-house scripts, where 0 represented homozygous Angus (AA), 1 represented the heterozygote state (AB/BA) and 2 represented homozygous Brahman (BB). The pipeline utilized is available under https://github.com/gzayasPR/BOA_Estimation (Zayas et al., 2024b), along with further details on data processing and management (Tange, 2011; Knaus and Grünwald, 2017; R Core Team, 2023; Barrett et al., 2024; Wickham et al., 2024).
2.4 Genome-wide association studies
GWAS was conducted to assess SNP and BOA effects on the phenotypic trait of interest using the following linear mixed models:
where: y denotes a vector of phenotypic records, SNP represents the SNP effect, modeled as a fixed effect encoded as a linear covariable as 0, 1, 2 (0 = homozygous for the major allele, 1 = heterozygous and 2 = homozygous for the minor allele). BOA is the Breed of Origin of Allele effect, as a fixed effect encoded as a linear covariable as 0, 1, 2 (0 = homozygous for the Angus BOA, 1 = heterozygous BOA and 2 = homozygous for the Brahman BOA). X is the incidence matrix, connecting phenotypic records to fixed effects, b is a vector of the fixed effects, including a categorical collection group, Z is the incidence matrix linking random genetic effects to phenotypes,
The genomic relationship matrix G was constructed based on the method proposed by VanRaden (VanRaden, 2008)
To minimize bias from proximal markers, a Leave-One-Chromosome-Out (LOCO) approach was implemented, where the
The single-marker GWAS was performed using the lmm.diago function from the gaston package (Perdry and Dandine-Roulland, 2023) in Rv4.4 (R Core Team, 2023). The lmm.diago function in the gaston package was chosen for its computational efficiency in handling large genomic datasets and its robust implementation of linear mixed models for GWAS. Marker significance was determined via single or joint Wald tests, and results were represented as -log(p-values).
To correct for multiple testing, a Bonferroni correction was implemented with an alpha level of 0.1 to control for family-wise error rate. In addition, a suggestive threshold was set at 1/number of tests to identify markers that showed potential associations but did not meet the stringent genome-wide significance criteria. GWAS results were processed using custom R scripts (R Core Team, 2023) and plotted using ggplot2 (Wickham, 2009). For mapping SNPs to specific genes, the Ensembl version 107 (Cunningham et al., 2022) and the ARS-UCD 1.2.0 genome assembly (Zimin et al., 2009) were utilized. Pipeline for GWAS can be found in https://github.com/gzayasPR/BOAxSNP-GWAS.git.
2.5 Functional enrichment analysis
To investigate the biological relevance of genomic regions associated with thermotolerance traits, functional enrichment analyses were performed using gene sets located near significant and suggestive SNPs. For each trait, genomic regions ±50 kb around the lead SNPs were annotated using Ensembl’s Bos taurus gene build (release July 2023) via the biomaRt package in R (Durinck et al., 2005; Durinck et al., 2009; R Core Team, 2023).
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the clusterProfiler package in R, with Bos taurus as the reference organism (Yu et al., 2012; Xu et al., 2024; Yu, 2024). GO enrichment was conducted across all three categories: Biological Process, Molecular Function, and Cellular Component. KEGG enrichment was carried out to identify overrepresented signaling and metabolic pathways. In both analyses, terms were considered significant at a false discovery rate (FDR) threshold of <0.1, and only gene sets with a minimum of three genes were included. All analyses were implemented in R (R Core Team, 2023), and trait-specific GO and KEGG results were saved as CSV files for reproducibility.
3 Results and discussion
The findings of this study highlight the value of BOA in understanding the genetic architecture of thermotolerance traits, in crossbred beef cattle. Unlike single-breed systems, crossbred cattle are characterized by unique patterns of LD that arise from the combination of parental breeds with distinct evolutionary histories. These complex LD patterns often present challenges for traditional SNP-based methods in detecting QTLs. BOA overcomes these limitations by leveraging ancestral information, enabling more precise identification of genomic regions and their associated effects.
3.1 Phenotypes
Table 1 provides descriptive statistics for all traits analyzed. The traits exhibited considerable variability, as indicated by their coefficient of variation, which ranged from 28.24% (LHL) to 43.57% (SGA). The wide range of values for each trait underscores the biological variation within the dataset, highlighting its potential to provide meaningful insights into association analyses and potential selection.
Table 1. Summary of phenotypic data.
3.2 GWAS long and short hair length
Figure 1 illustrates the GWAS Manhattan plot for the LHL trait, with the corresponding QTL intervals and their significance levels summarized in Table 2. For detailed information of SNP and BOA effects see Supplementary Table S1. In total, thirteen QTL were identified for LHL across various chromosomes (BTA). Of these, eight QTL, located on BTA 1 (42.7–57.2 Mb and 118.6–154.0 Mb), BTA 5, BTA 6 (55.7–60.8 Mb), BTA 19, BTA 20, and BTA 24, exceeded the Bonferroni significance threshold. Additional QTL on BTA 2, 6 (21.7–22.5 Mb), 7, 10, and 16 reached suggestive levels of significance.
Figure 1. Manhattan plots from GWAS for long hair length using (a) SNP-based, (b) BOA-based, and (c) SNP + BOA combined analyses. The -log10(p-values) are plotted across chromosomes, with green and red dashed lines representing the suggestive threshold (1/number of SNPs) and genome-wide significance threshold (Bonferroni correction, α = 0.1), respectively.
Table 2. Summary of significant and suggestive QTLs for traits related to short hair length, long hair length, sweat gland area, and thermal stress slope.
Figure 2 displays the GWAS Manhattan plot for the SHL trait, with associated QTL intervals and significance levels provided in Table 2. For detailed information of SNP and BOA effects see Supplementary Table S2. Seven QTL were identified for SHL. Three QTL, located on BTA 3 (19.8–22.4 Mb), BTA 12, and BTA 20 surpassed the Bonferroni threshold and were classified as significant. The remaining four QTL located on BTA 3 (26.1 Mb region), BTA 18, BTA 19 and BTA 24 reached suggestive levels of significance.
Figure 2. Manhattan plots from GWAS for short hair length using (a) SNP-based, (b) BOA-based, and (c) SNP + BOA combined analyses. The -log10(p-values) are plotted across chromosomes, with green and red dashed lines representing the suggestive threshold (1/number of SNPs) and genome-wide significance threshold (Bonferroni correction, α = 0.1), respectively.
The QTL on BTA20 was by far the most significant peak across both hair length traits. The alleles resulting in shorter hairs for both LHL and SHL were predominantly of Brahman origin, indicating a strong genetic contribution from this breed. However, the same allele was also observed in some Angus haplotypes, suggesting shared genetic variation between the breeds in this region of the genome. This QTL on BTA20, encompasses the prolactin receptor (PRLR) gene, which is strongly associated with hair length traits critical for thermoregulation (Sarlo Davila et al., 2019; Sosa et al., 2022). Shorter hair minimizes insulation and enhances heat dissipation, making this region a valuable target for breeding programs in hot climates (Olson et al., 2003; Dikmen et al., 2014; Sarlo Davila et al., 2019; Evangelista Façanha et al., 2020; Hansen, 2020). PRLR plays a role in the slick-hair phenotype observed in Criollo cattle, known to enhance thermotolerance (Dikmen et al., 2008; Huson et al., 2014; Porto-Neto et al., 2018; Sosa et al., 2021). The SLICK mutation in PRLR has been successfully introgressed into Holstein cattle, with the goal to improve temperature regulation in Holstein cattle in Florida, Unites States (Dikmen et al., 2008; Zayas et al., 2024a). Additionally, the PRLR gene has recently been successfully gene edited into Angus calves in an experimental setting, with the goal of improving thermotolerance in Angus cattle (Cuellar et al., 2024). Although the mutation identified by our study differs from the classic SLICK mutations, it likely performs a similar function in reducing hair length. This mutation was previously reported in a subset of this population, but its ancestral origin remained unclear (Sarlo Davila et al., 2020). Our model incorporating both BOA and SNP information found that a significant allele substitution effect −1.02 ± 0.13 whilst the BOA effect was 0.00 ± 0.11. Indicating that this effect for BOA was due to the higher prevalence of the beneficial allele rather than a BOA effect.
A unique QTL of interest for SHL is the QTL on BTA3 (19.7–22.3 Mb), overlapping ECM1, a biologically relevant gene encoding an extracellular matrix glycoprotein (Sercu et al., 2008). Mutations in ECM1 are linked to alopecia, attributed to compromised skin basement membrane integrity and impaired hair follicle health (Vahidnezhad et al., 1993).
A suggestive SHL QTL was identified on BTA19 (59.87–63.43 Mb), where the BOA model revealed a significant Angus substitution effect (0.22 ± 0.05, rs42969983). Supplementary Figure S1 shows a more detailed Manhattan plot of this region and boxplots showing the effect of BOA. This region contains several functionally relevant genes related to hair development and cycle, including ABCA5, PRKCA, and WIPI1. ABCA5 and its paralogs (ABCA6, 9, 10) are highly expressed in skin and hair follicles, with mutations linked to hypertrichosis and elongated follicles (DeStefano et al., 2014). Mutations in ABCA5 are known to cause excessive hair overgrowth, with affected individuals exhibiting significantly longer and thicker hair follicles (DeStefano et al., 2014). PRKCA encodes PKC-α, a kinase upregulated during hair growth and expressed in the outer root sheath of anagen follicles, its overexpression in mouse models has been associated with accelerated keratinocyte proliferation and enhanced hair follicle growth (Li et al., 2003). WIPI1, an autophagy regulator, is critical for hair follicle cell survival and timing of hair cycle transitions (Parodi et al., 2018; Sun et al., 2024).
The QTL on BTA19 for LHL demonstrates the unique value of BOA in capturing genetic signals that are otherwise missed by traditional SNP-based analyses. This QTL exhibited significant Angus substitution effects 0.56 ± 0.12 for LHL (rs41922675), with Angus showing considerably longer hair relative to Brahman. Supplementary Figure S1 shows a more detailed Manhattan plot of this region and boxplots showing the effect of BOA. These findings align with well-documented breed differences in hair length (Sarlo Davila et al., 2019, 2020). Several genes in this region are involved in hair development and thermoregulation. These genes include ACE and MAP3K3 which are involved in blood flow regulation and signaling pathways in response to stress. ACE plays a critical role in the renin-angiotensin system (RAS), which regulates blood flow, vascular tone, and inflammatory responses (Pontremoli et al., 1996; Karnik et al., 2009). Its association with hair follicle biology is underscored by its involvement in androgenetic alopecia (Karnik et al., 2009; Heymann, 2022). During stress, ACE redirects blood flow to vital organs by changing vaso- dilation and constriction, reducing blood supply to hair follicles thus resulting in shorter hair length (Yang et al., 1996; Saavedra and Benicky, 2007). Similarly, MAP3K3 acts as an upstream regulator of stress-activated protein kinases, such as the c-Jun N-terminal kinase (JNK) and p38 MAPK pathways. These pathways influence hair follicle biology by regulating keratinocyte function and inflammatory responses, which are crucial for hair growth and cycling (Chuang et al., 2000; Botchkareva et al., 2006; Quan et al., 2019). Dysregulated MAP3K3 activity could lead to premature hair follicle regression, impacting hair growth and length. This is further supported by the fact that the KEGG enrichment analysis on suggestive and significant genes for LHL, found enrichment for the MAPK signaling pathway. Another gene in this region, CCDC47 has been linked to wooly hair in humans, highlighting its potential role in hair structure and morphology (Morimoto et al., 2018). While the precise mechanism remains unclear, its involvement in calcium ion homeostasis may influence keratinocyte function and hair follicle health (Morimoto et al., 2018). Given that Brahman origin is beneficial in this region and there is a presence of multiple biologically relevant genes, we speculate that the observed effect arises from several advantageous variants distributed across these genes rather than a single causative mutation This complexity explains why the BOA analysis, which captures the combined effects of ancestral haplotypes could identified this region. Whereas the traditional SNP analysis which relies on individual marker effects or markers in LD (often disrupted in crossbred populations) failed to detect the effect.
3.3 GWAS sweat gland area
The GWAS results for SGA are presented in Figure 3, with putative QTLs summarized in Table 2 and significant markers detailed in Supplementary Table S3. A total of ten QTLs were identified, with three significant QTLs on BTA 5,13 and 20 and 7 QTLs across BTA 5,7,15,19,2 and 26. These QTLs were predominantly associated with BOA effects, with some regions also showing combined SNP and BOA contributions. Across all identified QTLs, Angus BOA inheritance was consistently associated with a decrease in SGA, as shown in Supplementary Table S3. This pattern underscores the ancestral influence of Angus genetics on sweat gland morphology, suggesting a potential evolutionary adaptation or breed-specific trait contributing to this characteristic.
Figure 3. Manhattan plots from GWAS for sweat gland area using (a) SNP-based, (b) BOA-based, and (c) SNP + BOA combined analyses. The -log10(p-values) are plotted across chromosomes, with green and red dashed lines representing the suggestive threshold (1/number of SNPs) and genome-wide significance threshold (Bonferroni correction, α = 0.1), respectively.
Modeling the SNP and BOA effects simultaneously, demonstrated superior resolution in identifying QTLs compared to the SNP-only and the BOA-only models. This integrated approach revealed combined effects in genomic regions that were less distinct when analyzed using either SNP or BOA effects independently. By accounting for both additive and breed-specific contributions, we were able to gain deeper insights into the genetic architecture of complex traits such as SGA.
The significant QTL on BTA 5 includes several genes with known roles in skin and sweat gland development, notably WIF1, GRIP1, MSRB3, and LEMD3. WIF1 encodes a protein that antagonizes Wnt signaling by binding to Wnt proteins (Hsieh et al., 1999). As Wnt signaling is essential for the formation and differentiation of ectodermal appendages, including sweat glands (Lu et al., 2012), dysregulation of WIF1 may impact sweat gland morphogenesis, size, and function (Coates et al., 2019). GRIP1 encodes a scaffolding protein that organizes transmembrane receptors at the cell membrane. It plays a crucial role in epidermal basement membrane attachment. Mutations in GRIP1 cause Fraser syndrome, characterized by skin syndactyly and dermo-epidermal junction defects; in mouse models, loss of GRIP1 results in epidermal blistering due to impaired dermal-epidermal junction integrity (Vogel et al., 2012).
MSRB3 encodes a methionine sulfoxide reductase involved in repairing oxidatively damaged proteins, thereby protecting skin cells from oxidative stress. It is expressed in vascular endothelial cells of the skin and contributes to dermal integrity and aging processes (Taungjaruwinai et al., 2009). Lastly, LEMD3 encodes the inner nuclear membrane protein MAN1, which antagonizes TGF-β/BMP signaling through Smad binding. This regulation is critical for connective tissue and bone homeostasis. Loss-of-function mutations in LEMD3 cause Buschke-Ollendorff syndrome, characterized by connective tissue nevi and osteopoikilosis, indicating a requirement for LEMD3 in normal dermal connective tissue development (Hellemans et al., 2004). BOA effects showed that Angus ancestry significantly reduced SGA (−23.28 ± 4.37,rs210687570), consistent with previous findings linking Angus genetics to smaller sweat glands and reduced heat dissipation efficiency (Mateescu et al., 2023; Hernandez et al., 2024). This suggests that targeting Brahman-derived alleles could be advantageous for improving thermotolerance in cattle.
3.4 GWAS thermal stress slope
TSS reflects an animal’s ability to maintain stable core body temperature under heat stress, serving as a critical metric for thermotolerance. Animals with lower TSS values demonstrate better thermoregulation, experiencing smaller fluctuations in body temperature under heat stress. In contrast, higher TSS values indicate reduced thermoregulatory capacity due to less efficient responses to heat load. This trait, which accounts for daily body temperature variations and cumulative heat stress exposure, is pivotal for understanding the genetic basis of thermoregulation and guiding breeding strategies aimed at enhancing heat resilience.
The GWAS results for TSS, presented in Figure 4, identified nine QTLs, of which 4 significant QTLs were found on BTA1. Five suggestive QTLs were found on BTA 2,19,18 and 28. Supplementary Table S4 provides details of suggestive and significant markers.
Figure 4. Manhattan plots from GWAS for thermal stress slope (TSS) using (a) SNP-based, (b) BOA-based, and (c) SNP + BOA combined analyses. The -log10(p-values) are plotted across chromosomes, with green and red dashed lines representing the suggestive threshold (1/number of SNPs) and genome-wide significance threshold (Bonferroni correction, α = 0.1), respectively.
The QTL on BTA1(44.0–58.5 Mb) spans a large genomic region (Table 2) containing multiple genes involved in stress response and adaptation. Among these, FILIP1L and ALCAM are particularly notable for their interactions with heat shock proteins and pathways related to immune and heat stress responses. The FILIP1L gene encode for a Heat shock factor (Hsf1) interacting protein, which reduces Hsf1-mediated transcription Hsf1 by promoting Hsf1 ubiquitination and degradation (Hu and Mivechi, 2011). Heats shock proteins are vital to heat stress response, these results indicate that Brahman BOA may have variants that better regulate the expression of heat shock proteins. The Activated Leukocyte Cell Adhesion Molecule (ALCAM) gene plays a role in immune cell adhesion and migration, which is critical for immune function, particularly under stress conditions like heat exposure (Zimmerman et al., 2006; te Riet et al., 2007; Chan et al., 2015). Heat stress can upregulate integrin adhesion molecules, which interact closely with ALCAM to enhance leukocyte adhesion and migration (Chan et al., 2015; Tukaj, 2020). Heat shock protein 70 (HSP70), a molecular chaperone, plays a critical role in stabilizing cell adhesion proteins during thermal stress, enhancing leukocyte function and adhesion (te Riet et al., 2007; Chan et al., 2015; McRobb et al., 2019). HSP70 stabilizes cell adhesion proteins, potentially enhancing leukocyte functionality and adhesion under heat stress conditions (Chan et al., 2015; Rigg et al., 2016). The interaction between ALCAM, integrins, and HSP70 supports immune function and helps mitigate thermal damage, contributing to immune resilience during heat stress. Other genes related to immune and stress response include the NFKBIZ and DCBLD2 genes. NFKBIZ functions as a regulator in the NF-κB signaling pathway, which is crucial for immune and stress responses (Wong et al., 1997; Malhotra et al., 2002). Alterations in NF-κB signaling have been linked to heat stress responses in livestock. DCBLD2, a transmembrane protein with roles in cellular proliferation and stress responses, has been implicated in oxidative stress pathways a consequence of heat stress (Chen et al., 2021; Wang et al., 2024). The Angus substitution effect of this QTL was estimated at 0.27 ± 0.05 (rs3423104800), indicating that Angus inheritance has a detrimental impact on thermotolerance. A more detailed Manhattan plot of the QTL alongside other QTLs on BTA 1 can be found in Supplementary Figure S2. The observed effect of Angus alleles in this region, associated with an unfavorable increase in TSS, suggests that Brahman BOA may offer genetic advantages for heat resilience. Integrating the SNP and BOA effects, demonstrated superior resolution in disentangling the overlapping contributions of these factors. This integrated approach not only improves the precision of QTL mapping but also offers deeper insights into how breed-specific and additive genetic factors interact to influence complex traits like thermoregulatory capacity under heat stress.
3.5 Functional enrichment analysis
To further explore the biological relevance of the genomic regions associated with thermotolerance traits, gene ontology (GO) and KEGG pathway enrichment analyses were conducted for the set of candidate genes identified within significant and suggestive QTLs for each trait.
For LHL, enriched GO biological processes included keratinocyte differentiation (GO:0030216, FDR = 0.001), skin development (GO:0043588, FDR = 0.001), and keratinization (GO:0031424, FDR = 0.003). These terms reflect the biological underpinnings of hair growth and structure. Additionally, KEGG pathway analysis revealed enrichment for the MAPK signaling pathway (bta04010, FDR = 0.047) and estrogen signaling pathway (bta04915, FDR = 0.045), both implicated in skin and hair follicle regulation (Oh and Smart, 1996; Ohnemus et al., 2005; Akilli Öztürk et al., 2015; Lu et al., 2021). For SHL, GO analysis highlighted terms related to epidermis morphogenesis (GO:0048730), cell-cell adhesion (GO:0098742), and adherens junction organization (GO:0034332), alongside molecular functions such as DNA-binding transcription repressor activity and calcium-activated chloride channel activity. KEGG results showed significant enrichment for adherens junctions (bta04520), consistent with the involvement of adhesion molecules and cytoskeletal elements in follicular structure and integrity (Young et al., 2003; Wu et al., 2014).
The SGA trait exhibited highly significant enrichment (FDR < 0.001) for multiple GO terms associated with epidermal and epithelial development, including epidermis development (GO:0008544), keratinocyte differentiation (GO:0030216), and intermediate filament organization (GO:0045109). These findings align with the morphological nature of sweat gland area and its dependence on dermal structure.
For TSS, GO terms associated with skin development (GO:0043588, FDR = 0.013), keratinization (GO:0031424, FDR = 0.013), and innate immune responses (GO:0002220, FDR = 0.038) were enriched. These results underscore the dual role of skin in thermoregulation and immune response under heat stress. KEGG pathways significantly enriched for TSS included the MAPK signaling pathway (bta04010, FDR = 0.047) and estrogen signaling pathway (bta04915, FDR = 0.045), both relevant to heat shock response and cellular adaptation (Stice and Knowlton, 2008; Kumar et al., 2020; Hu et al., 2023).
Overall, enrichment results corroborate the trait-specific biological processes and pathways highlighted in our GWAS findings. The consistent enrichment of skin- and keratin-related processes across LHL, SHL, and SGA, and immune/stress response pathways for TSS, reinforce the validity of the detected QTLs and support the functional relevance of the implicated genes in thermotolerance phenotypes.
4 Broader implications and limitations
This study builds upon previous research utilizing BOA, which demonstrated the value of ancestral information for genomic evaluations and QTL identification (Guillenea et al., 2023; Warburton et al., 2023; Hernandez et al., 2024; Zayas et al., 2024b). The integrated modeling strategy provides a more comprehensive understanding of genetic contributions, as demonstrated by the identification of suggestive and significant QTLs for thermotolerance traits. This approach differs from studies that have modeled SNPs in a breed-specific manner without directly incorporating BOA effects (Eiríksson et al., 2021; Guillenea et al., 2023; Warburton et al., 2023). By leveraging BOA, we identified significant QTLs for thermotolerance traits that were undetectable using SNP-only approaches. These findings highlight the critical role of breed-specific genetic contributions in shaping thermotolerance traits. For instance, Brahman-derived alleles were associated with improved thermoregulation, while Angus-derived alleles were linked to traits such as longer hair and reduced sweat gland area. This detailed understanding of ancestral effects provides actionable insights for breeding programs. Specifically, targeted selection for Brahman-derived alleles represents a promising strategy to enhance heat resilience in cattle in tropical climates.
Future studies are needed to investigate the mechanisms behind the BOA effects observed for these traits. We hypothesize that BOA may provide a more accurate model for QTLs in crossbred cattle compared to SNP-based approaches alone. Another possibility is that BOA has the potential to capture non-SNP factors, such as regulatory or epigenetic mechanisms, that influence gene expression and phenotypes. Because BOA assigns variants to their ancestral origin, it provides a framework for disentangling breed-specific regulatory effects in crossbred populations. Integrating BOA analyses with functional genomic tools such as transcriptomics, eQTL mapping, RNA-Seq, Chip-Seq or ATAC-Seq could enable the identification of regulatory elements or expression patterns that are unique to a particular ancestral background. This breed-aware regulatory mapping could uncover ancestry-specific transcription factor binding sites, epigenetic marks, or chromatin accessibility profiles that underlie observed phenotypic differences. Such analyses may offer a deeper understanding of the genetic and regulatory networks driving complex traits, ultimately improving biological interpretation and improving the precision of breeding strategies.
Despite its potential, the application of BOA has certain limitations. The resolution of BOA assignments is constrained by the density of SNP arrays, which may miss finer-scale ancestral contributions. Future studies utilizing whole-genome sequencing or high-density genotyping arrays could enhance the precision and accuracy of BOA analyses, providing a more detailed understanding of ancestral genetic influences. Furthermore, the computational complexity of integrating BOA with SNP-based approaches underscores the need for efficient analytical pipelines to support wider adoption.
Incorporating BOA into genomic studies represents a powerful tool for advancing our understanding of the genetic architecture of complex traits using crossbred and composite populations. By leveraging breed-specific ancestry information, this approach offers valuable insights into the complex interactions between ancestry and phenotype. It paves the way for a more nuanced understanding of the genetic dynamics within these populations, ultimately supporting more effective breeding and selection strategies.
5 Conclusion
This study demonstrates the power of BOA for dissecting the genetic architecture of thermotolerance using crossbred cattle. By integrating BOA with SNP effects, we uncovered novel QTL and ancestral contributions that are critical for improving thermoregulation in cattle. For LHL and SHL, a major QTL on BTA20 overlapping PRLR showed strong associations with shorter hair, with beneficial alleles predominantly of Brahman origin. For SGA, Angus BOA consistently decreased gland size, with a significant BTA5 QTL containing genes such as WIF1, GRIP1, and LEMD3, all involved in skin structure and development. The most significant TSS QTL on BTA1 harbored FILIP1L and ALCAM, which regulate immune response and interact with heat shock proteins during thermal stress. Functional enrichment revealed MAPK and estrogen signaling pathways were common to both LHL and TSS, underscoring shared regulatory mechanisms. These results highlight BOA’s ability to uncover breed-specific genetic effects missed by SNP-only models. These findings provide valuable insights for breeding programs aiming to enhance heat resilience of beef cattle in tropical and subtropical climates. Future research should aim to validate these findings and explore functional mechanisms underlying the observed QTL, maximizing their potential for use in genetic improvement strategies.
Data availability statement
The genotypic data is available in the European Variation Archive repository, PRJEB75981. Scripts used for analysis are available on Github https://github.com/gzayasPR/BOA_GWAS and https://github.com/gzayasPR/BOAxSNP-GWAS.git.
Ethics statement
The animal study was approved by University of Florida Institutional Animal Care and Use Committee number 202200000019. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
GZ: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing. CS: Data curation, Investigation, Writing – review and editing. ER: Data curation, Formal Analysis, Writing – review and editing, Conceptualization. AH: Conceptualization, Data curation, Writing – review and editing. AB: Data curation, Writing – review and editing. FR: Writing – review and editing, Data curation. KS: Conceptualization, Data curation, Investigation, Methodology, Writing – review and editing. RM: Conceptualization, Data curation, Formal Analysis, Methodology, Supervision, Writing – review and editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This research was supported by UF Agricultural Experimental Station Hatch FLA-ANS-005548, USDA-NIFA Grants 2017-67007-26143, 2020-67015-30820 and 2021-67015-33414. Gabriel Zayas is a USDA National Needs Fellow supported by USDA NIFA grant 2019-38420-28977 and Town Creek Farms.
Acknowledgments
Authors thank the owners and staff of the Seminole Tribe of Florida for their participation in this study. Additionally, thanks are extended to Zoetis Animal Health (Kalamazoo, MI) and Mitchell Blanding for providing blank controlled internal drug releases.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1576966/full#supplementary-material
References
Akilli Öztürk, Ö., Pakula, H., Chmielowiec, J., Qi, J., Stein, S., Lan, L., et al. (2015). Gab1 and mapk signaling are essential in the hair cycle and hair follicle stem cell quiescence. Cell Rep. 13, 561–572. doi:10.1016/j.celrep.2015.09.015
Animals, N. R. C. C. P. E. E. F. (1971). A Guide to environmental research on animals. National Academy of Sciences. Washington D.C. Available online at: https://books.google.com/books?id=gzsrAAAAYAAJ.
Baran, Y., Pasaniuc, B., Sankararaman, S., Torgerson, D. G., Gignoux, C., Eng, C., et al. (2012). Fast and accurate inference of local ancestry in Latino populations. Bioinformatics 28, 1359–1367. doi:10.1093/bioinformatics/bts144
Barrett, T., Dowle, M., Srinivasan, A., Gorecki, J., Chirico, M., Hocking, T., et al. (2024). data.table: extension of `data.frame`. Available online at: https://CRAN.R-project.org/package=data.table.
Bolormaa, S., Pryce, J. E., Kemper, K. E., Hayes, B. J., Zhang, Y., Tier, B., et al. (2013). Detection of quantitative trait loci in bos indicus and bos taurus cattle using genome-wide association studies. Genet. Sel. Evol. 45, 43. doi:10.1186/1297-9686-45-43
Botchkareva, N. V., Ahluwalia, G., and Shander, D. (2006). Apoptosis in the hair follicle. J. Investigative Dermatology 126, 258–264. doi:10.1038/sj.jid.5700007
Carvajal, M. A., Alaniz, A. J., Gutiérrez-Gómez, C., Vergara, P. M., Sejian, V., and Bozinovic, F. (2021). Increasing importance of heat stress for cattle farming under future global climate scenarios. Sci. Total Environ. 801, 149661. doi:10.1016/j.scitotenv.2021.149661
Chan, Y.-C., Greenwood, D. R., Yang, Y., Leung, E., and Krissansen, G. W. (2015). Leukocyte integrin α4β7 associates with heat shock protein 70. Mol. Cell Biochem. 409, 263–269. doi:10.1007/s11010-015-2530-z
Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaSci 4, 7. doi:10.1186/s13742-015-0047-8
Chen, X., Lv, Y., Xu, K., Wang, X., Zhao, Y., Li, J., et al. (2021). DCBLD2 mediates epithelial-mesenchymal transition-induced metastasis by cisplatin in lung adenocarcinoma. Cancers 13, 1403. doi:10.3390/cancers13061403
Chuang, S.-M., Wang, I.-C., and Yang, J.-L. (2000). Roles of JNK, p38 and ERK mitogen-activated protein kinases in the growth inhibition and apoptosis induced by cadmium. Carcinogenesis 21, 1423–1432. doi:10.1093/carcin/21.7.1423
Claudia, A. S., Sevillano, C. A., Bovenhuis, H., Bovenhuis, H., Bovenhuis, H., Calus, M. P. L., et al. (2019). Genomic evaluation for a crossbreeding system implementing breed-of-origin for targeted markers. Front. Genet. 10, 418. doi:10.3389/fgene.2019.00418
Coates, M., Mariottoni, P., Corcoran, D. L., Kirshner, H. F., Jaleel, T., Brown, D. A., et al. (2019). The skin transcriptome in hidradenitis suppurativa uncovers an antimicrobial and sweat gland gene signature which has distinct overlap with wounded skin. PLOS ONE 14, e0216249. doi:10.1371/journal.pone.0216249
Cuellar, C. J., Amaral, T. F., Rodriguez-Villami, P., Ongaratto, F. L., Martinez, D. O., Labrecque, R., et al. (2024). Consequences of gene editing of PRLR on thermotolerance, growth, and male reproduction in cattle. FASEB BioAdvances 6, 223–234. doi:10.1096/fba.2024-00029
Cunningham, F., Allen, J. E., Allen, J., Alvarez-Jarreta, J., Amode, M. R., Armean, I. M., et al. (2022). Ensembl 2022. Nucleic Acids Res. 50, D988–D995. doi:10.1093/nar/gkab1049
DeStefano, G. M., Kurban, M., Anyane-Yeboa, K., Dall’Armi, C., Paolo, G. D., Feenstra, H., et al. (2014). Mutations in the cholesterol transporter gene ABCA5 are associated with excessive hair overgrowth. PLOS Genet. 10, e1004333. doi:10.1371/journal.pgen.1004333
Dikmen, S., Alava, E., Pontes, E., Fear, J. M., Dikmen, B. Y., Olson, T. A., et al. (2008). Differences in thermoregulatory ability between slick-haired and wild-type lactating Holstein cows in response to acute heat stress. J. Dairy Sci. 91, 3395–3402. doi:10.3168/jds.2008-1072
Dikmen, S., and Hansen, P. J. (2009). Is the temperature-humidity index the best indicator of heat stress in lactating dairy cows in a subtropical environment? J. Dairy Sci. 92, 109–116. doi:10.3168/jds.2008-1370
Dikmen, S., Khan, F. A., Huson, H. J., Sonstegard, T. S., Moss, J. I., Dahl, G. E., et al. (2014). The SLICK hair locus derived from Senepol cattle confers thermotolerance to intensively managed lactating Holstein cows. J. Dairy Sci. 97, 5508–5520. doi:10.3168/jds.2014-8087
Durinck, S., Moreau, Y., Kasprzyk, A., Davis, S., De Moor, B., Brazma, A., et al. (2005). BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics 21, 3439–3440. doi:10.1093/bioinformatics/bti525
Durinck, S., Spellman, P. T., Birney, E., and Huber, W. (2009). Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat. Protoc. 4, 1184–1191. doi:10.1038/nprot.2009.97
Eiríksson, J. H., Karaman, E., Su, G., and Christensen, O. F. (2021). Breed of origin of alleles and genomic predictions for crossbred dairy cows. Genet. Sel. Evol. 53, 84. doi:10.1186/s12711-021-00678-3
Evangelista Façanha, D. A., Leite, J. H. G. M., Rufino Sousa, J. E., Ferreira, J. B., Costa, W. P., Jannuzzi Valente, F. L., et al. (2020). Phenotypic selection of Brahman bulls: adaptive typification for breeding programs. J. Therm. Biol. 93, 102736. doi:10.1016/j.jtherbio.2020.102736
Guillenea, A., Lund, M. S., Evans, R., Boerner, V., and Karaman, E. (2023). A breed-of-origin of alleles model that includes crossbred data improves predictive ability for crossbred animals in a multi-breed population. Genet. Sel. Evol. 55, 34. doi:10.1186/s12711-023-00806-1
Hansen, P. J. (2020). Prospects for gene introgression or gene editing as a strategy for reduction of the impact of heat stress on production and reproduction in cattle. Theriogenology 154, 190–202. doi:10.1016/j.theriogenology.2020.05.010
Hellemans, J., Preobrazhenska, O., Willaert, A., Debeer, P., Verdonk, P. C. M., Costa, T., et al. (2004). Loss-of-function mutations in LEMD3 result in osteopoikilosis, Buschke-Ollendorff syndrome and melorheostosis. Nat. Genet. 36, 1213–1218. doi:10.1038/ng1453
Hernandez, A. S., Zayas, G. A., Rodriguez, E. E., Sarlo Davila, K. M., Rafiq, F., Nunez, A. N., et al. (2024). Exploring the genetic control of sweat gland characteristics in beef cattle for enhanced heat tolerance. J. Animal Sci. Biotechnol. 15, 66. doi:10.1186/s40104-024-01025-4
Heymann, W. R. (2022). The inflammatory component of androgenetic alopecia. J. Am. Acad. Dermatology 86, 301–302. doi:10.1016/j.jaad.2021.11.013
Hsieh, J. C., Kodjabachian, L., Rebbert, M. L., Rattner, A., Smallwood, P. M., Samos, C. H., et al. (1999). A new secreted protein that binds to Wnt proteins and inhibits their activities. Nature 398, 431–436. doi:10.1038/18899
Hu, Y., Lin, L., Liu, K., Liu, E., Han, S., Gong, Z., et al. (2023). l -Theanine alleviates heat stress-induced impairment of immune function by regulating the p38 MAPK signalling pathway in mice. Food Funct. 14, 335–343. doi:10.1039/D2FO02775E
Hu, Y., and Mivechi, N. F. (2011). Promotion of heat shock factor Hsf1 degradation via adaptor protein filamin A-interacting protein 1-like (FILIP-1L). J. Biol. Chem. 286, 31397–31408. doi:10.1074/jbc.M111.255851
Huson, H. J., Kim, E.-S., Godfrey, R. W., Olson, T. A., McClure, M. C., Chase, C. C., et al. (2014). Genome-wide association study and ancestral origins of the slick-hair coat in tropically adapted cattle. Front. Genet. 5, 101. doi:10.3389/fgene.2014.00101
Karnik, P., Tekeste, Z., McCormick, T. S., Gilliam, A. C., Price, V. H., Cooper, K. D., et al. (2009). Hair follicle stem cell-specific PPARγ deletion causes scarring alopecia. J. Investigative Dermatology 129, 1243–1257. doi:10.1038/jid.2008.369
Knaus, B. J., and Grünwald, N. J. (2017). VCFR: a package to manipulate and visualize variant call format data in R. Mol. Ecol. Resour. 17, 44–53. doi:10.1111/1755-0998.12549
Koufariotis, L., Hayes, B. J., Kelly, M., Burns, B. M., Lyons, R., Stothard, P., et al. (2018). Sequencing the mosaic genome of Brahman cattle identifies historic and recent introgression including polled. Sci. Rep. 8, 17761. doi:10.1038/s41598-018-35698-5
Kumar, R. R., Arora, K., Goswami, S., Sakhare, A., Singh, B., Chinnusamy, V., et al. (2020). MAPK enzymes: a ROS activated signaling sensors involved in modulating heat stress response, tolerance and grain stability of wheat under heat stress. 3 Biotech. 10, 380. doi:10.1007/s13205-020-02377-0
Li, L.-F., Guo, J., and Gao, Z.-F. (2003). Overexpression of skin protein kinase C-alpha in anagen hair follicles during induced growth of mouse hair. Clin. Exp. Dermatol 28, 429–433. doi:10.1046/j.1365-2230.2003.01328.x
Lu, C. P., Polak, L., Rocha, A. S., Pasolli, H. A., Chen, S.-C., Sharma, N., et al. (2012). Identification of stem cell populations in sweat glands and ducts reveals roles in homeostasis and wound repair. Cell 150, 136–150. doi:10.1016/j.cell.2012.04.045
Lu, Q., Gao, Y., Fan, Z., Xiao, X., Chen, Y., Si, Y., et al. (2021). Amphiregulin promotes hair regeneration of skin-derived precursors via the PI3K and MAPK pathways. Cell Prolif. 54, e13106. doi:10.1111/cpr.13106
Malhotra, V., Eaves-Pyles, T., Odoms, K., Quaid, G., Shanley, T. P., and Wong, H. R. (2002). Heat shock inhibits activation of NF-kappaB in the absence of heat shock factor-1. Biochem. Biophysical Res. Commun. 291, 453–457. doi:10.1006/bbrc.2002.6470
Mateescu, R. G., Sarlo-Davila, K. M., Dikmen, S., Rodriguez, E., and Oltenacu, P. A. (2020). The effect of Brahman genes on body temperature plasticity of heifers on pasture under heat stress. J. Animal Sci. 98, skaa126. doi:10.1093/jas/skaa126
Mateescu, R. G., Sarlo Davila, K. M., Hernandez, A. S., Andrade, A. N., Zayas, G. A., Rodriguez, E. E., et al. (2023). Impact of Brahman genetics on skin histology characteristics with implications for heat tolerance in cattle. Front. Genet. 14, 1107468. doi:10.3389/fgene.2023.1107468
McRobb, L. S., McKay, M. J., Gauden, A. J., Lee, V. S., Subramanian, S., Thomas, S. G., et al. (2019). Radiation-stimulated translocation of CD166 and CRYAB to the endothelial surface provides potential vascular targets on irradiated brain arteriovenous malformations. Int. J. Mol. Sci. 20, 5830. doi:10.3390/ijms20235830
Morimoto, M., Waller-Evans, H., Ammous, Z., Song, X., Strauss, K. A., Pehlivan, D., et al. (2018). Bi-Allelic CCDC47 variants cause a disorder characterized by woolly hair, liver dysfunction, dysmorphic features, and global developmental delay. Am. J. Hum. Genet. 103, 794–807. doi:10.1016/j.ajhg.2018.09.014
Oh, H. S., and Smart, R. C. (1996). An estrogen receptor pathway regulates the telogen-anagen hair follicle transition and influences epidermal cell proliferation. Proc. Natl. Acad. Sci. 93, 12525–12530. doi:10.1073/pnas.93.22.12525
Ohnemus, U., Uenalan, M., Conrad, F., Handjiski, B., Mecklenburg, L., Nakamura, M., et al. (2005). Hair cycle control by estrogens: catagen induction via estrogen receptor (ER)-α is checked by ERβ signaling. Endocrinology 146, 1214–1225. doi:10.1210/en.2004-1219
Olson, T. A., Lucena, C., Chase, C. C., and Hammond, A. C. (2003). Evidence of a major gene influencing hair length and heat tolerance in Bos taurus cattle. J. Animal Sci. 81, 80–90. doi:10.2527/2003.81180x
Parodi, C., Hardman, J. A., Allavena, G., Marotta, R., Catelani, T., Bertolini, M., et al. (2018). Autophagy is essential for maintaining the growth of a human (mini-)organ: evidence from scalp hair follicle organ culture. PLOS Biol. 16, e2002864. doi:10.1371/journal.pbio.2002864
Paşaniuc, B., Sankararaman, S., Kimmel, G., and Halperin, E. (2009). Inference of locus-specific ancestry in closely related populations. Bioinformatics 25, i213–i221. doi:10.1093/bioinformatics/btp197
Perdry, H., and Dandine-Roulland, C. (2023). Class, D. B. (C++ gzstream, and class), L. K. (C++ gzstream gaston): genetic data handling (QC, GRM, LD, PCA) and linear mixed models. Available online at: https://cran.r-project.org/web/packages/gaston/index.html (Accessed December 4, 2024).
Pontremoli, R., Sofia, A., Tirotta, A., Ravera, M., Nicolella, C., Viazzi, F., et al. (1996). The deletion polymorphism of the angiotensin I-converting enzyme gene is associated with target organ damage in essential hypertension. J. Am. Soc. Nephrol. 7, 2550–2558. doi:10.1681/ASN.V7122550
Porto-Neto, L. R., Bickhart, D. M., Landaeta-Hernandez, A. J., Utsunomiya, Y. T., Pagan, M., Jimenez, E., et al. (2018). Convergent evolution of slick coat in cattle through truncation mutations in the prolactin receptor. Front. Genet. 9, 57. doi:10.3389/fgene.2018.00057
Quan, V. L., Zhang, B., Mohan, L. S., Shi, K., Isales, M. C., Panah, E., et al. (2019). Activating structural alterations in MAPK genes are distinct genetic drivers in a unique subgroup of spitzoid neoplasms. Am. J. Surg. Pathology 43, 538–548. doi:10.1097/PAS.0000000000001213
R Core Team (2023). R: a Language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/.
Rigg, R. A., Healy, L. D., Nowak, M. S., Mallet, J., Thierheimer, M. L. D., Pang, J., et al. (2016). Heat shock protein 70 regulates platelet integrin activation, granule secretion and aggregation. Am. J. Physiology-Cell Physiology 310, C568–C575. doi:10.1152/ajpcell.00362.2015
Saavedra, J. M., and Benicky, J. (2007). Brain and peripheral angiotensin II play a major role in stress. Stress 10, 185–193. doi:10.1080/10253890701350735
Sarlo Davila, K. M., Hamblen, H., Hansen, P. J., Dikmen, S., Oltenacu, P. A., and Mateescu, R. G. (2019). Genetic parameters for hair characteristics and core body temperature in a multibreed Brahman–Angus herd1. J. Animal Sci. 97, 3246–3252. doi:10.1093/jas/skz188
Sarlo Davila, K. M., Howell, A., Nunez, A., Orelien, A., Roe, V., Rodriguez, E., et al. (2020). Genome-wide association study identifies variants associated with hair length in Brangus cattle. Anim. Genet. 51, 811–814. doi:10.1111/age.12970
Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M., Pietzsch, T., et al. (2012). Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682. doi:10.1038/nmeth.2019
Sercu, S., Zhang, M., Oyama, N., Hansen, U., Ghalbzouri, A. E. L., Jun, G., et al. (2008). Interaction of extracellular matrix protein 1 with extracellular matrix components: ECM1 is a basement membrane protein of the skin. J. Invest Dermatol 128, 1397–1408. doi:10.1038/sj.jid.5701231
Sosa, F., Carmickle, A. T., Jiménez-Cabán, E., Ortega, M. S., Dikmen, S., Negrón-Pérez, V., et al. (2021). Inheritance of the SLICK1 allele of PRLR in cattle. Anim. Genet. 52, 887–890. doi:10.1111/age.13145
Sosa, F., Santos, J. E. P., Rae, D. O., Larson, C. C., Macchietto, M., Abrahante, J. E., et al. (2022). Effects of the SLICK1 mutation in PRLR on regulation of core body temperature and global gene expression in liver in cattle. animal 16, 100523. doi:10.1016/j.animal.2022.100523
Stice, J. P., and Knowlton, A. A. (2008). Estrogen, NFkappaB, and the heat shock response. Mol. Med. 14, 517–527. doi:10.2119/2008-00026.Stice
Sun, P., Wang, Z., Li, S., Yin, J., Gan, Y., Liu, S., et al. (2024). Autophagy induces hair follicle stem cell activation and hair follicle regeneration by regulating glycolysis. Cell Biosci. 14, 6. doi:10.1186/s13578-023-01177-2
Tange, O. (2011). GNU parallel - the command-line power tool. login USENIX Mag. 36, 42–47. doi:10.5281/zenodo.16303
Taungjaruwinai, W. M., Bhawan, J., Keady, M., and Thiele, J. J. (2009). Differential expression of the antioxidant repair enzyme methionine sulfoxide reductase (MSRA and MSRB) in human skin. Am. J. Dermatopathol. 31, 427–431. doi:10.1097/DAD.0b013e3181882c21
te Riet, J., Zimmerman, A. W., Cambi, A., Joosten, B., Speller, S., Torensma, R., et al. (2007). Distinct kinetic and mechanical properties govern ALCAM-mediated interactions as shown by single-molecule force spectroscopy. J. Cell Sci. 120, 3965–3976. doi:10.1242/jcs.004010
Tukaj, S. (2020). Heat shock protein 70 as a double agent acting inside and outside the cell: insights into autoimmunity. Int. J. Mol. Sci. 21, 5298. doi:10.3390/ijms21155298
Utsunomiya, Y. T., Milanesi, M., Fortes, M. R. S., Porto-Neto, L. R., Utsunomiya, A. T. H., Silva, M. V. G. B., et al. (2019). Genomic clues of the evolutionary history of Bos indicus cattle. Anim. Genet. 50, 557–568. doi:10.1111/age.12836
Vahidnezhad, H., Youssefian, L., and Uitto, J. (1993). “Lipoid Proteinosis,” in GeneReviews®. Editors M. P. Adam, J. Feldman, G. M. Mirzaa, R. A. Pagon, S. E. Wallace, and A. Amemiya (Seattle (WA): University of Washington, Seattle). Available online at: http://www.ncbi.nlm.nih.gov/books/NBK338540/(Accessed March 30, 2025).
Vandenplas, J., Calus, M. P. L., Sevillano, C. A., Windig, J. J., and Bastiaansen, J. W. M. (2016). Assigning breed origin to alleles in crossbred animals. Genet. Sel. Evol. 48, 61. doi:10.1186/s12711-016-0240-y
VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. Dairy Sci. 91, 4414–4423. doi:10.3168/jds.2007-0980
Verdugo, M. P., Mullin, V. E., Scheu, A., Mattiangeli, V., Daly, K. G., Maisano Delser, P., et al. (2019). Ancient cattle genomics, origins, and rapid turnover in the Fertile Crescent. Science 365, 173–176. doi:10.1126/science.aav1002
Vogel, M. J., van Zon, P., Brueton, L., Gijzen, M., van Tuil, M. C., Cox, P., et al. (2012). Mutations in GRIP1 cause Fraser syndrome. J. Med. Genet. 49, 303–306. doi:10.1136/jmedgenet-2011-100590
Wang, J., Zhang, H., Chen, L., Fu, K., Yan, Y., and Liu, Z. (2024). CircDCBLD2 alleviates liver fibrosis by regulating ferroptosis via facilitating STUB1-mediated PARK7 ubiquitination degradation. J. Gastroenterol. 59, 229–249. doi:10.1007/s00535-023-02068-6
Warburton, C. L., Costilla, R., Engle, B. N., Moore, S. S., Corbet, N. J., Fordyce, G., et al. (2023). Concurrently mapping quantitative trait loci associations from multiple subspecies within hybrid populations. Heredity 131, 350–360. doi:10.1038/s41437-023-00651-4
Wickham, H. (2009). ggplot2: elegant graphics for data analysis. New York, NY: Springer New York. doi:10.1007/978-0-387-98141-3
Wickham, H., Vaughan, D., and Girlich, M. (2024). Tidyr: tidy messy data. Available online at: https://CRAN.R-project.org/package=tidyr.
Wong, H. R., Ryan, M., and Wispé, J. R. (1997). The heat shock response inhibits inducible nitric oxide synthase gene expression by blocking iκ-B degradation and NF-κB nuclear translocation. Biochem. Biophysical Res. Commun. 231, 257–263. doi:10.1006/bbrc.1997.6076
Wu, M., Guo, X., Yang, L., Wang, Y., Tang, Y., Yang, Y., et al. (2014). Mesenchymal stem cells with modification of junctional adhesion molecule A induce hair formation. Stem Cells Transl. Med. 3, 481–488. doi:10.5966/sctm.2013-0165
Xu, S., Hu, E., Cai, Y., Xie, Z., Luo, X., Zhan, L., et al. (2024). Using clusterProfiler to characterize multiomics data. Nat. Protoc. 19, 3292–3320. doi:10.1038/s41596-024-01020-z
Yang, G., Wan, Y., and Zhu, Y. (1996). Angiotensin II- an important stress hormone. Neurosignals 5, 1–8. doi:10.1159/000109168
Young, P., Boussadia, O., Halfter, H., Grose, R., Berger, P., Leone, D. P., et al. (2003). E-cadherin controls adherens junctions in the epidermis and the renewal of hair follicles. EMBO J. 22, 5723–5733. doi:10.1093/emboj/cdg560
Yu, G. (2024). Thirteen years of clusterProfiler. Innovation 5, 100722. doi:10.1016/j.xinn.2024.100722
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for comparing biological themes among gene clusters. OMICS A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zayas, G. A., Dikmen, S., Mateescu, R. G., and Hansen, P. J. (2024a). Maintaining breed integrity: successful introgression of the SLICK1 allele into the Holstein breed. J. Hered., esae057. doi:10.1093/jhered/esae057
Zayas, G. A., Rodriguez, E., Hernandez, A., Rezende, F. M., and Mateescu, R. G. (2024b). Breed of origin analysis in genome-wide association studies: enhancing SNP-based insights into production traits in a commercial Brangus population. BMC Genomics 25, 654. doi:10.1186/s12864-024-10465-1
Zimin, A. V., Delcher, A. L., Florea, L., Kelley, D. R., Schatz, M. C., Puiu, D., et al. (2009). A whole-genome assembly of the domestic cow, Bos taurus. Genome Biol. 10, R42. doi:10.1186/gb-2009-10-4-r42
Keywords: thermotolerance, PRLR, heat stress, GWAS, crossbreeding
Citation: Zayas GA, Santos Rojas C, Rodriguez EE, Hernandez AS, Beard AM, Rafiq F, Sarlo Davila KM and Mateescu RG (2025) Genetic architecture of thermotolerance traits in beef cattle: a novel integration of SNP and breed-of-origin effects. Front. Genet. 16:1576966. doi: 10.3389/fgene.2025.1576966
Received: 14 February 2025; Accepted: 22 April 2025;
Published: 30 April 2025.
Edited by:
Aline Silva Mello Cesar, University of São Paulo, BrazilReviewed by:
Bárbara Silva-Vignato, University of São Paulo, BrazilShamik Polley, West Bengal University of Animal and Fishery Sciences, India
Copyright © 2025 Zayas, Santos Rojas, Rodriguez, Hernandez, Beard, Rafiq, Sarlo Davila and Mateescu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Gabriel A. Zayas, Z3pheWFzOTdAdWZsLmVkdQ==