Genetic Markers for Improved Fertility and Thermotolerance in Holstein Dairy Cows Exposed to Environmental Heat Stress

1. Introduction

Climate change is intensifying environmental conditions that compromise livestock productivity, primarily by increasing the number of hot days and heat waves throughout the year due to rising global temperatures [1]. These conditions significantly increase the risk of HS in cattle, which is characterized by an elevated heart rate, increased respiratory rate, and higher body temperature. In particular, dairy cows are highly vulnerable to HS, as their elevated metabolic heat production amplifies thermal effects under high ambient temperatures [2,3]. In semi-arid regions of Mexico, extremely high summer temperatures induce HS, which markedly compromises reproductive efficiency in Holstein cows that calve in late spring. As a result, most producers postpone reproductive management until late autumn, resulting in a six-month delay in insemination and subsequent calving, which particularly impacts high-producing cows.

Heat stress in dairy cows negatively affects reproduction by disrupting the hypothalamic–pituitary–ovarian axis, reducing estradiol, LH, and progesterone levels, which compromises estrus expression, corpus luteum formation, and embryo viability [4]. High temperatures also damage oocytes, reducing embryonic development and requiring several estrous cycles for recovery. As a result, pregnancy rates decline markedly during summer, with reductions of up to 24% [5,6]. In addition to its direct effects on reproductive physiology, HS elicits systemic endocrine and metabolic adaptations that further exacerbate declines in productive and reproductive performance. Thermal stress reduces thyroid hormone concentrations as a physiological mechanism to lower metabolic rate [7]. This metabolic down-regulation is accompanied by decreased feed intake, resulting in limited systemic nutrient availability, which adversely affects both milk production and reproductive functions [8,9,10].

Holstein cows with favorable genotypes for genes within the PRL and GH/IGF1 pathways exhibited the ability to maintain reproductive success under HS conditions, suggesting a genetic basis for thermo-tolerance and fertility [11]. Similarly, a genomic marker association study reported SNPs in the genes AMH, IGFBP1, LGR5, and TLR4 as predictors for fertility traits and physiological variables indicative of HS tolerance [12]. The identification and validation of these genetic markers may offer novel opportunities to enhance reproductive efficiency during the summer months, thereby helping to maintain consistent annual milk production cycles in dairy cows. This approach could also serve as a viable strategy to mitigate the detrimental effects of HS in dairy herds and promote overall herd productivity and resilience.

To identify cows with greater reproductive potential in hot climates, a comprehensive understanding of their adaptive responses to thermal challenges is required. Significant individual variation in responses to HS has been reported, even among animals within the same group, breed, and environmental conditions, thereby highlighting opportunities for improvement through genetic selection [13]. With the advancement of bovine genome sequencing and high-density genotyping technologies for SNPs, the genome-wide association studies (GWAS) have become an essential tool for identifying genes associated with economically important traits [14]. In several dairy cattle breeds, GWAS have been successfully used to identify genomic regions involved in heat tolerance, fertility, and milk production [15,16].

Traditional selective breeding for heat tolerance in dairy cattle shows slow genetic progress due to long generation intervals, low heritability, and difficulty to measure HS-related phenotypes; however, in the face of increasing HS, faster genetic gains are required. Genomics and DNA marker technologies enable the selection of young animals and accelerate genetic gain for thermotolerance [17]. Therefore, genomic and marker-assisted technologies have been proposed as novel and effective molecular tools to improve genetic selection for relevant phenotypic traits [18]. Marker-assisted selection (MAS) approaches, based on GWAS findings, enable the validation of SNPs in independent populations to identify animals carrying favorable alleles at specific markers or genomic regions linked to economically important traits, including heat tolerance and fertility, in dairy cows exposed to HS conditions [19,20,21].

Therefore, our objective was to identify and validate genetic markers associated with fertility and thermotolerance traits in lactating Holstein cows exposed to environmental conditions that induce HS in northwest Mexico, using genome-wide and marker-assisted technologies.

2. Materials and Methods

The Institutional Animal Care and Use Committee of the Instituto Tecnológico de Sonora approved all procedures performed in animals (Approval code 2017-0079).

2.1. Location and Experimental Population

The study was conducted on three neighboring dairy farms located in the Yaqui Valley, Sonora, Mexico (27°21′ N, 109°54′ W). These farms operated under an intensive dairy production system characterized by confined housing, high stocking density, and mechanized milking. Climatic conditions throughout the year were characterized by wide variability in ambient temperature (−2.5 to 45.9 °C), relative humidity (25–90%), and solar radiation (0.05 to 761.7 W/m²). Annual precipitation averaged approximately 341 mm, with daily rainfall ranging from 3 to 90 mm on rainy days.

A total of 400 Holstein cows, aged 4 to 6 years, were included in the study. Cows had an average body weight of 645.1 ± 32.5 kg and an average body condition score of 3.8 ± 0.01. Only spring-calving cows (i.e., cows calving between March and June) were selected for the study because they were exposed to HS during summer reproductive management. All cows were managed under similar conditions, housed in shaded facilities with free access to water and a commercial trace mineral supplement. Flooring and shaded areas were adapted to meet the requirements of Holstein cows. To meet the nutritional requirements of lactating dairy cattle, cows were fed a total mixed ration twice daily, formulated for an average body weight of 650 kg and a milk yield of approximately 30 kg/day, with an average milk composition of 3.5% fat and 3.2% crude protein.

2.2. Reproductive Management

Cows were subjected to an ovulation synchronization protocol that included the insertion of an intravaginal progesterone-releasing device (CIDR®, Pfizer, Mexico) on day 0, along with the intramuscular (IM) administration of 0.01 mg of GnRH (Fertagyl®, Intervet, Mexico). The CIDR was removed 7 days later, and cows received an IM injection of 25 mg of prostaglandin F2α (PGF2α; Lutalyse®, Pfizer, Mexico). A second dose of 0.01 mg of GnRH was administered on day 9, and cows were inseminated by fixed-time artificial insemination (FTAI) 24 h later (day 10).

Pregnancy diagnosis was performed by transrectal ultrasonography 30 days after FTAI using a 7.5 MHz linear transducer (SonoSite MicroMaxx™, Bothell, WA, USA). Cows diagnosed as non-pregnant received an IM dose of 25 mg of PGF2α at that time, followed by 0.01 mg of GnRH 48 h later, and were inseminated again by FTAI 24 h later. Ultrasonographic pregnancy diagnosis was repeated 30 days after FTAI, and the protocol was repeated in non-pregnant cows until pregnancy was achieved. Once pregnancy was confirmed, the number of FTAI services required per cow was recorded to calculate services per conception (SPC). Days open (DO) were determined as the number of days from calving to confirmed conception, minus the 30 days corresponding to the interval between the onset of pregnancy and its ultrasonography confirmation.

2.3. Sampling and SNP Quality Control

A blood sample (3 ml) was drawn from each cow through venipuncture of the coccygeal vein using disposable sterile syringes. Five drops of the collected whole blood were spotted onto Fast Technology for Analysis of Nucleic Acids cards (FTA®), which were sent to Neogen AgriGenomics (Lincoln, NE) for DNA extraction and genotyping. The SNP panel BovineSNP50, which contains 53,218 highly informative SNPs evenly dispersed across the bovine genome, was used to obtain genotypes of each cow.

PLINK v1.07 software was used to perform quality control on the SNPs tested in a genotype to phenotype association analyses. Only SNPs that met the following criteria were included in analyses: 1) call rate of greater than 95% or a false discovery rate of less than 5%, 2) missing genotype frequency of less than 5%, 3) minor allele frequency (MAF) of more than 5%, and 4) no deviation from Hardy-Weinberg equilibrium (p-value of Chi-square goodness-of-fit test greater than 0.05, X² > 0.05). The analyses also excluded SNPs with unmapped loci or those on the sex chromosomes. After quality control, a total of 45,832 SNPs were retained for further analyses.

2.4. Genome-Wide Association Study (GWAS)

A principal component analysis (PCA) was applied to correct for batch effects and population stratification in the input test data. Subsequently, a single-locus mixed model (i.e., single-marker SNP-based GWAS) was used for genomic analysis to investigate associations between each SNP marker genotype and reproductive traits (i.e., DO and SPC) as the phenotypic observations. The GWAS was conducted on an SNP-by-SNP basis using the SNP & Variation Suite software, version 8.8.1 (SVSv8; Golden Helix, Inc., Bozeman, MT, USA; www.goldenhelix.com; accessed on 26 Jul 2025).

The additive mixed model applied was: y = Xβ + Za + e, where y represented the vector of phenotypic observations (i.e., DO and SPC); X was the incidence matrix of fixed effects; β was a vector of fixed effects, including herd, number of lactations, and the additive effect of the candidate SNP tested for association; Z was the design matrix for random additive genetic effects; a was the vector of random additive genetic effects; and e was the vector of residual effects. An animal-specific random effect was included in the model to account for population structure. Under this model, it was assumed that a ~ N(0, Gσ²ₐ) and e ~ N(0, Iσ²ₑ), where σ²ₐ represented the additive genetic variance, σ²ₑ the residual variance, G the genomic relationship matrix, and I the identity matrix.

2.5. Multiple Testing Adjustment

To correct for multiple comparisons, P-values derived from the genomic analyses were adjusted using the Bonferroni correction (b = α/n), which assumes independence among SNPs. The number of tests (n) was set equal to the total number of informative SNPs (n = 45,832), with an experiment-wise error rate of α = 0.05. A nominal P-value of 1.09 × 10⁻⁶ corresponded to the 5% genome-wide significance threshold, equivalent to a value of 5.96 on the −log₁₀ (p) scale.

All SNPs exceeding the Bonferroni-adjusted threshold were subsequently subjected to pairwise linkage disequilibrium (LD) analysis using PROC ALLELE (SAS software v9.4). Linkage disequilibrium was estimated using the coefficient of determination (R²), which represents the squared correlation between alleles at two loci. When a pair of SNPs showed a high correlation, a multicollinearity test was performed to assess the linear relationship between them.

2.6. SNP Validation Study

An independent Holstein cattle population, including 184 spring-calving cows from two neighboring commercial dairy farms, was selected. This study was conducted to validate the candidate SNP markers identified by the GWAS as being associated with the reproductive traits DO and SPC. The cows were lactating and had an average body weight of 651.3 ± 34.7 kg, an average body condition score of 3.0 ± 0.02, and were between 4 and 7 years of age. The animals were managed throughout the summer and were therefore exposed to environmental conditions associated with HS.

At 60 days postpartum, cows underwent a reproductive tract examination to confirm uterine involution. Subsequently, they were subjected to the previously described ovulation synchronization protocol, which included FTAI. Ovarian structures were monitored by transrectal ultrasonography using a 7.5 MHz transducer (SonoSite MicroMaxx™, Bothell, WA, USA). The diameter of the dominant follicle (FOL) was measured on the day of FTAI, and the size of the corpus luteum (CL) was recorded on day 18 after FTAI. Transrectal pregnancy diagnosis was performed by ultrasonography at 30 days after FTAI.

Cows diagnosed as non-pregnant received an intramuscular (IM) dose of 25 mg of PGF₂α, followed by an IM dose of 0.01 mg of GnRH 48 h later, and were inseminated again by FTAI 24 h later. Pregnancy diagnosis was repeated 30 days after each subsequent FTAI until pregnancy was achieved, allowing calculating DO and SPC.

2.7. Physiological Traits

Throughout the experimental period, rectal temperature (RT; °C) and respiratory rate (RR; breaths/min) were recorded twice weekly at 0700 h and 1700 h. RT was measured using a digital thermometer (TES-1310®) equipped with a type K contact probe (9 cm length), which was inserted to ensure contact with the rectal mucosa, whereas RR was assessed by visually counting flank movements. The same trained graduate student consistently collected all physiological measurements during the study.

2.8. SNP Selection and Genotyping

In the previously described genomic study, 21 SNPs with significant effects on days open (DO) and 20 SNPs associated with services per conception (SPC) were reported. Five SNPs were identified as common predictors of both reproductive variables, which constituted the primary criterion for their selection. These SNPs were rs133890206, rs136074864, rs109479159, rs43520457, and rs133770320, located within the genes IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1, respectively.

A blood sample was collected from each cow at the beginning of the study via jugular venipuncture using EDTA-coated Vacutainer tubes (EDTA-Na₂; Venoject®, Terumo, Lakewood, CA, USA). Samples were centrifuged at 3500 rpm for 15 min, after which 200 μL of the buffy coat was aspirated using a pipette and stored at −20 °C. Genomic DNA was subsequently extracted using a commercial kit (DNeasy Blood & Tissue Kit; QIAGEN, Hilden, Germany) following the manufacturer’s instructions. DNA concentration and purity were assessed using an automated NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

The five significant SNPs selected for the validation study were genotyped using the TaqMan allelic discrimination assay and RT-qPCR (StepOne™, Applied Biosystems, Foster City, CA, USA), according to the methodology described by Castillo-Salas et al. [22]. PCR amplification was performed using the StepOne Real-Time PCR System (Thermo Fisher Scientific, MA, USA). Finally, SNP genotyping and PCR data analysis were conducted using StepOne software (Life Technologies Corporation, version 2.3).

2.9. Statistics for the Validation Analysis

The MEANS procedure was used to calculate descriptive statistics for the reproductive traits FOL, CL, DO, and SPC, as well as for the thermo-tolerance traits RT and RR. The UNIVARIATE and GLM (Levene’s test) procedures were applied to assess the assumptions of normality of data distribution and homogeneity of variances, respectively. Allelic and genotypic frequencies were calculated using PROC ALLELE. Hardy–Weinberg equilibrium (HWE) was evaluated using a Chi-Square (χ²) test implemented in PROC FREQ. All statistical analyses were performed using SAS software (version 9.4; SAS Institute Inc., Cary, NC, USA).

A mixed-effects statistical model was used to evaluate the association between genotype and phenotype, employing the MIXED procedure for continuous variables. The association model was defined as: Y_ijklm = μ + A_i + B_j + C_k + D_l + E_m + e_ijklm, which included the response variable (Y_ijklm), the fixed effects of the SNP genotype (A_i), dairy herd (B_j), and number of lactations (C_k), the linear covariate of days in milk (D_l), the random effect of sire (E_m), and the residual error term (e_ijklm). When the association analysis identified the term genotype as a significant source of variation, preplanned pairwise comparisons of least squares means were performed using the PDIFF option. Mean separation tests were conducted using the LSMEANS statement with Bonferroni adjustment [23]. The MIXED procedure was also used to estimate average allele substitution effects by regressing the phenotype on the number of copies of a given SNP allele, included as a covariate. Additive and dominance (i.e., non-additive) genetic effects were also calculated following the methodology described by Falconer and Mackay [24].

2.10. mRNA Expression

The five genes corresponding to the candidate SNPs (i.e., IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1) were assessed for gene expression by quantitative polymerase chain reaction (qPCR) in a subset of Holstein cows (n = 30). Oligonucleotide primers were designed based on the bovine genome using Primer-BLAST. Optimal annealing temperatures and primer efficiencies were determined by PCR. Primer specificity was confirmed by nucleotide sequencing of PCR products cloned into the PCR 2.1-TOPO vector (Thermo Fisher Scientific Life Sciences, Waltham, MA, USA).

Complementary DNA (cDNA) was synthesized from 0.4 μg of total RNA per reaction using SuperScript III (Thermo Fisher Scientific). PCR amplification was performed using SYBR Green chemistry (Qiagen, Valencia, CA, USA) on an iQ5 Real-Time PCR Detection System (Bio-Rad Lab, Hercules, CA, USA). All qPCR results were normalized to the geometric mean of the reference gene ribosomal protein S15. Quantitative analysis was conducted using the comparative Ct method, and relative fold changes were calculated using the 2^−ΔΔCt method. Student’s t-test was used to compare results between cows with favorable (n = 16) and non-favorable SNP genotypes (n = 14).

2.11. Hormonal Assays

A 3 mL blood sample was collected using Vacutainer tubes with yellow caps and serum separator on the day of FTAI, and again 18 days later to determine serum concentrations of anti-Müllerian hormone (AMH) and progesterone, respectively. Additional blood samples were collected three times per week throughout the study to assess serum levels of cortisol and heat shock protein 70 (HSP70).

Serum AMH concentrations were measured using a bovine-specific AMH immunoassay (AMH Fertility Assay™, Minitube of America), which has a sensitivity of 0.04 pg/mL and an intra-assay coefficient of variation (CV) of 2.2%. Commercial ELISA kits were used to determine serum progesterone and cortisol concentrations (bovine progesterone and bovine cortisol; MyBioSource LLC). The progesterone assay exhibited a sensitivity of 0.2 ng/mL and an intra-assay CV of 11.3%, whereas the cortisol assay showed a sensitivity of 5 pg/mL and an intra-assay CV of 8.9%. Serum HSP70 levels were evaluated using a bovine HSP70 ELISA kit (Bioassay Technology Laboratory, China), with a sensitivity of 0.23 ng/mL.

Finally, a Pearson correlation analysis was performed to assess the associations between serum concentrations of reproductive hormones (AMH and progesterone) and stress-related biomarkers (cortisol and HSP70) with the expression levels of the genes IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1 in Holstein cows experiencing summer heat stress.

3. Results

3.1. Climate Data

Environmental conditions experienced by the cows in this study were evaluated through seasonal fluctuations in the temperature-humidity index (THI) (Figure 1). Average THI values were 70.1, 77.1, 80.3, and 82.2 units for May, June, July, and August, respectively.

3.2. Genome-Wide Association Study (GWAS)

After quality control, a total of 45,832 SNPs were retained for the GWAS. Using a single-locus mixed model, 72 SNPs were initially identified across the genome as associated with days open (DO). After correction for multiple testing, only 21 SNPs remained statistically significant, exceeding the Bonferroni-adjusted threshold of 1.09 × 10⁻⁶ (Figure 2). These results highlight genomic regions that may influence reproductive performance in Holstein cows and provide a basis for further functional studies.

Of the 21 SNPs significantly associated with DO, 17 were intronic variants located within known genes, whereas 4 were intergenic. The intronic SNPs were rs42416336, rs133573430, rs133770320, rs133890206, rs441125342, rs136074864, rs42468270, rs43520457, rs109740021, rs109479159, rs137336947, rs110893810, rs43613018, rs134587965, rs111023809, rs133073533, and rs137801036, which were mapped to the genes LONRF1, PARK7, TRHDE1, IQGAP1, LGR5, EXOC4, MIER2, GPR151, COL19A1, GLI3, SULF1, RSPH14, RMND1, HECW1, GUCYB1, NKX2-1, and HMGCLL1, respectively. The intergenic SNPs rs134756540, rs137461577, rs135454149, and rs109769851 were not associated with any annotated gene (Table 1).

Figure 2. Manhattan plot from single-marker GWAS for the trait days open (DO) in Holstein cows managed during summer in a semi-desert region.

Figure 2. Manhattan plot from single-marker GWAS for the trait days open (DO) in Holstein cows managed during summer in a semi-desert region.

Table 1. Significant SNPs (p < 1.09 × 10⁻⁶) from a single-marker genome-wide association study (GWAS) associated with DO in Holstein cows managed in a warm environment.

SNP ID 1	Variant 2	BTA 3	Position 4	Gene 5	Alleles 6	Var 7	p-Value 8
rs42416336	Intronic	27	24’273969	LONRF1	T/G	1.62	8.10 × 10−10
rs133573430	Intronic	16	45’413480	PARK7	A/G	1.62	1.64 × 10−9
rs133770320	Intronic	5	2′119998	TRHDE1	T/C	1.57	8.47 × 10−9
rs133890206	Intronic	21	22′177205	IQGAP1	A/C	1.52	1.06 × 10−8
rs441125342	Intronic	5	1′105051	LGR5	T/C	1.43	1.11 × 10−8
rs136074864	Intronic	4	97’452287	EXOC4	A/G	1.24	1.18 × 10−8
rs134756540	Intergenic	4	78′201917	--------	C/T	1.39	1.61 × 10−8
rs42468270	Intronic	7	43′047351	MIER2	A/T	1.31	7.08 × 10−8
rs43520457	Intronic	7	57′753560	GPR151	C/T	1.29	7.51 × 10−8
rs109740021	Intronic	9	91′141195	COL19A1	G/A	1.24	9.21 × 10−8
rs109479159	Intronic	4	78′759970	GLI3	A/G	1.23	1.12 × 10−7
rs137336947	Intronic	14	33’418670	SULF1	G/A	1.19	1.26 × 10−7
rs137461577	Intergenic	4	75′730827	--------	C/T	1.18	1.54 × 10−7
rs110893810	Intronic	17	71′821119	RSPH14	C/T	1.17	1.58 × 10−7
rs43613018	Intronic	9	88′450249	RMND1	T/G	1.16	2.43 × 10−7
rs134587965	Intronic	4	77′979197	HECW1	T/C	1.15	2.58 × 10−7
rs135454149	Intergenic	2	82′616155	--------	T/C	1.15	2.59 × 10−7
rs111023809	Intronic	17	43′600087	GUCYB1	G/A	1.15	3.04 × 10−7
rs133073533	Intronic	21	46′720726	NKX2-1	G/A	1.15	3.26 × 10−7
rs137801036	Intronic	23	4′831520	HMGCLL1	A/G	1.15	5.25 × 10−7
rs109769851	Intergenic	5	80′626601	--------	G/T	1.15	5.27 × 10−7

¹ SNP reference of the NCBI; ² type of SNP chromosome variant; ³ Bos taurus autosomal chromosome number; ⁴ SNP position within the chromosome; ⁵ candidate gene symbol; ⁶ alleles from the SNP; ⁷ percentage of trait variance explained by the SNP; ⁸ SNP statistical significance.

Table 1. Significant SNPs (p < 1.09 × 10⁻⁶) from a single-marker genome-wide association study (GWAS) associated with DO in Holstein cows managed in a warm environment.

SNP ID 1	Variant 2	BTA 3	Position 4	Gene 5	Alleles 6	Var 7	p-Value 8
rs42416336	Intronic	27	24’273969	LONRF1	T/G	1.62	8.10 × 10−10
rs133573430	Intronic	16	45’413480	PARK7	A/G	1.62	1.64 × 10−9
rs133770320	Intronic	5	2′119998	TRHDE1	T/C	1.57	8.47 × 10−9
rs133890206	Intronic	21	22′177205	IQGAP1	A/C	1.52	1.06 × 10−8
rs441125342	Intronic	5	1′105051	LGR5	T/C	1.43	1.11 × 10−8
rs136074864	Intronic	4	97’452287	EXOC4	A/G	1.24	1.18 × 10−8
rs134756540	Intergenic	4	78′201917	--------	C/T	1.39	1.61 × 10−8
rs42468270	Intronic	7	43′047351	MIER2	A/T	1.31	7.08 × 10−8
rs43520457	Intronic	7	57′753560	GPR151	C/T	1.29	7.51 × 10−8
rs109740021	Intronic	9	91′141195	COL19A1	G/A	1.24	9.21 × 10−8
rs109479159	Intronic	4	78′759970	GLI3	A/G	1.23	1.12 × 10−7
rs137336947	Intronic	14	33’418670	SULF1	G/A	1.19	1.26 × 10−7
rs137461577	Intergenic	4	75′730827	--------	C/T	1.18	1.54 × 10−7
rs110893810	Intronic	17	71′821119	RSPH14	C/T	1.17	1.58 × 10−7
rs43613018	Intronic	9	88′450249	RMND1	T/G	1.16	2.43 × 10−7
rs134587965	Intronic	4	77′979197	HECW1	T/C	1.15	2.58 × 10−7
rs135454149	Intergenic	2	82′616155	--------	T/C	1.15	2.59 × 10−7
rs111023809	Intronic	17	43′600087	GUCYB1	G/A	1.15	3.04 × 10−7
rs133073533	Intronic	21	46′720726	NKX2-1	G/A	1.15	3.26 × 10−7
rs137801036	Intronic	23	4′831520	HMGCLL1	A/G	1.15	5.25 × 10−7
rs109769851	Intergenic	5	80′626601	--------	G/T	1.15	5.27 × 10−7

¹ SNP reference of the NCBI; ² type of SNP chromosome variant; ³ Bos taurus autosomal chromosome number; ⁴ SNP position within the chromosome; ⁵ candidate gene symbol; ⁶ alleles from the SNP; ⁷ percentage of trait variance explained by the SNP; ⁸ SNP statistical significance.

The single-locus mixed model identified 86 SNPs across the genome associated with services per conception (SPC), of which only 20 remained significant after multiple testing correction, exceeding the Bonferroni-adjusted threshold of 1.09 × 10⁻⁶ (Figure 3).

Of the 20 SNPs associated with SPC, 15 were intronic variants, whereas 6 were intergenic. The intronic SNPs were rs110893810, rs133890206, rs43520457, rs133770320, rs109161550, rs137194049, rs136074864, rs109479159, rs110147189, rs110234465, rs134571905, rs43517716, rs42124106, rs1337126241, and rs8193046, mapped to RSPH14, IQGAP1, GPR151, TRHDE1, RASAL2, DISC1, EXOC4, GLI3, DENND5B, KDM7A, lncRNA, TCERG1, lncRNA, SCN8A, and TLR4, respectively. The intergenic SNPs were rs43475092, rs43768368, rs29024666, rs109769851, and rs136745124, which were not associated with any annotated gene (Table 2).

Figure 3. Manhattan plot from single-marker GWAS for services for the trait conception (SPC) in Holstein cows managed during summer in a semi-desert region.

Figure 3. Manhattan plot from single-marker GWAS for services for the trait conception (SPC) in Holstein cows managed during summer in a semi-desert region.

Table 2. Significant SNPs (P < 1.09 × 10⁻⁶) from a single-marker genome-wide association study (GWAS) associated with SPC in Holstein cows managed in a warm environment.

SNP ID 1	Variant 2	BTA 3	Position 4	Gene 5	Alleles 6	Var 7	p-Value 8
rs43475092	Intergenic	3	67’819418	--------	G/A	1.62	2.35 × 10−9
rs110893810	Intronic	17	71′821119	RSPH14	C/T	1.57	7.53 × 10−9
rs133890206	Intronic	21	22′177205	IQGAP1	A/C	1.52	8.16 × 10−9
rs43768368	Intergenic	9	4′831793	--------	A/G	1.43	2.38 × 10−8
rs43520457	Intronic	7	57′753560	GPR151	C/T	1.29	1.20 × 10−7
rs133770320	Intronic	5	2′119998	TRHDE1	T/C	1.43	2.38 × 10−8
rs29024666	Intergenic	17	20’625888	--------	A/C	00	2.38 × 10−8
rs109161550	Intronic	16	59’872943	RASAL2	A/G	00	2.38 × 10−8
rs109769851	Intergenic	5	80’626601	--------	A/C	1.31	9.86 × 10−8
rs137194049	Intronic	28	4′599045	DISC1	T/C	1.29	1.20 × 10−7
rs136074864	Intronic	4	97’452287	EXOC4	A/G	1.24	2.38 × 10−7
rs109479159	Intronic	4	78′759970	GLI3	A/G	1.23	2.46 × 10−7
rs110147189	UGeneV	5	78’338501	DENND5B	T/C	1.19	4.15 × 10−7
rs110234465	DGeneV	4	103′622830	KDM7A	A/G	1.18	4.24 × 10−7
rs134571905	Intronic	8	11′947725	lncRNA	C/T	1.17	5.22 × 10−7
rs43517716	Intronic	7	57′699520	TCERG1	T/C	1.16	5.24 × 10−7
rs42124106	Intronic	27	25′327107	lncRNA	A/C	1.15	5.49 × 10−7
rs1337126241	Intronic	5	28′142346	SCN8A	G/A	1.15	5.49 × 10−7
rs8193046	Intronic	8	108’062912	TLR4	G/A	1.15	5.49 × 10−7
rs136745124	Intergenic	6	11’486713	--------	G/A	1.15	5.49 × 10−7

¹ SNP reference of the NCBI; ² type of SNP chromosome variant; ³ Bos taurus autosomal chromosome number; ⁴ SNP position within the chromosome; ⁵ candidate gene symbol. ⁶ alleles from the SNP; ⁷ percentage of trait variance explained by the SNP; ⁸ SNP statistical significance.

Table 2. Significant SNPs (P < 1.09 × 10⁻⁶) from a single-marker genome-wide association study (GWAS) associated with SPC in Holstein cows managed in a warm environment.

SNP ID 1	Variant 2	BTA 3	Position 4	Gene 5	Alleles 6	Var 7	p-Value 8
rs43475092	Intergenic	3	67’819418	--------	G/A	1.62	2.35 × 10−9
rs110893810	Intronic	17	71′821119	RSPH14	C/T	1.57	7.53 × 10−9
rs133890206	Intronic	21	22′177205	IQGAP1	A/C	1.52	8.16 × 10−9
rs43768368	Intergenic	9	4′831793	--------	A/G	1.43	2.38 × 10−8
rs43520457	Intronic	7	57′753560	GPR151	C/T	1.29	1.20 × 10−7
rs133770320	Intronic	5	2′119998	TRHDE1	T/C	1.43	2.38 × 10−8
rs29024666	Intergenic	17	20’625888	--------	A/C	00	2.38 × 10−8
rs109161550	Intronic	16	59’872943	RASAL2	A/G	00	2.38 × 10−8
rs109769851	Intergenic	5	80’626601	--------	A/C	1.31	9.86 × 10−8
rs137194049	Intronic	28	4′599045	DISC1	T/C	1.29	1.20 × 10−7
rs136074864	Intronic	4	97’452287	EXOC4	A/G	1.24	2.38 × 10−7
rs109479159	Intronic	4	78′759970	GLI3	A/G	1.23	2.46 × 10−7
rs110147189	UGeneV	5	78’338501	DENND5B	T/C	1.19	4.15 × 10−7
rs110234465	DGeneV	4	103′622830	KDM7A	A/G	1.18	4.24 × 10−7
rs134571905	Intronic	8	11′947725	lncRNA	C/T	1.17	5.22 × 10−7
rs43517716	Intronic	7	57′699520	TCERG1	T/C	1.16	5.24 × 10−7
rs42124106	Intronic	27	25′327107	lncRNA	A/C	1.15	5.49 × 10−7
rs1337126241	Intronic	5	28′142346	SCN8A	G/A	1.15	5.49 × 10−7
rs8193046	Intronic	8	108’062912	TLR4	G/A	1.15	5.49 × 10−7
rs136745124	Intergenic	6	11’486713	--------	G/A	1.15	5.49 × 10−7

¹ SNP reference of the NCBI; ² type of SNP chromosome variant; ³ Bos taurus autosomal chromosome number; ⁴ SNP position within the chromosome; ⁵ candidate gene symbol. ⁶ alleles from the SNP; ⁷ percentage of trait variance explained by the SNP; ⁸ SNP statistical significance.

3.3. SNPs Selection

Out of the 21 and 20 SNPs significantly associated with the reproductive traits DO and SPC, respectively, five intragenic SNPs common to both GWAS analyses were selected for validation. These SNPs were in Hardy–Weinberg equilibrium (HWE, χ² > 0.05) and met the minor allele frequency criterion (i.e., MAF > 0.10). The selected SNPs rs133890206, rs136074864, rs109479159, rs43520457, and rs133770320 were located within the genes IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1, respectively (Table 3). These five SNPs met the criteria for inclusion in a genotype to phenotype association study, and were suitable candidates for validation as molecular markers of fertility in Holstein cows exposed to heat stress.

3.4. SNP Markers Validation

The least-squares means for ovulatory follicle diameter (FOL), corpus luteum size (CL), days open (DO), and services per conception (SPC), according to the genotypes of the five significant SNPs are presented in Table 4. The SNPs rs136074864, rs109479159, rs43520457, and rs133770320 were associated with FOL, CL, SPC, and DO (P < 0.05), with the favorable genotypes being GG, AA, TT, and CC, respectively, due to their association with improved reproductive performance. However, the SNP rs133890206 was not associated (P > 0.05) with any of the evaluated traits (i.e., FOL, CL, DO, and SPC). Allelic substitution analysis confirmed the positive contribution of the favorable allele of each SNP to the reproductive traits, while fixed-effects analysis corroborated the additive effect of the genes corresponding to each of the four significant SNPs.

The least-squares means for thermotolerance-related traits, including rectal temperature (RT) and respiratory rate (RR), are presented in Table 5. The SNPs rs136074864, rs109479159, rs43520457, and rs133770320 were significantly associated with both RT and RR (P < 0.05), with GG, AA, TT, and CC, respectively, identified as the favorable genotypes based on their association with a more favorable physiological response to heat stress. In contrast, SNP rs133890206 was not associated (P > 0.05) with either of the thermotolerance-related traits evaluated in this study (i.e., RT and RR). Allelic substitution and fixed-effects analyses confirmed the positive contribution of the favorable alleles and the additive effect of the corresponding genes, respectively.

3.5. Quantitative RT-PCR Validation

The expression patterns of the candidate genes associated with reproductive and thermo-tolerance traits were consistent across all measurements, as shown in Figure 4. The relative mRNA expression levels (Log₂ Fold Change) for the genes IQGAP1, EXOC4, GLI3, GPR151, and RSPH14 in cows with favorable genotypes were 1.22, 3.64, 4.12, 4.76, and 3.24, respectively. In cows with non-favorable genotypes, the corresponding values were 0.92, 1.76, 1.54, 2.12, and 0.95, respectively. Gene expression profiles differed significantly (P < 0.05) between cows carrying favorable and non-favorable genotypes for EXOC4, GLI3, GPR151, and TRHDE1, but not for IQGAP1.

3.6. Serum Hormone Profiles

Serum levels of the reproductive hormones AMH and progesterone are presented in Figure 5. AMH levels were 0.86 and 0.38 ng/mL, and progesterone levels were 1.74 and 1.16 ng/mL in cows with favorable and non-favorable genotypes, respectively. Significant differences were observed between groups for both hormones (P < 0.05).

Serum levels of the stress-related hormones cortisol and HSP70 are presented in Figure 6. Cortisol levels were 18.66 and 24.74 ng/mL, and HSP70 levels were 26.10 and 24.91 ng/mL in cows with favorable and non-favorable genotypes, respectively. Significant differences were observed between groups only for cortisol (P < 0.05).

3.7. Correlation Between Hormone Levels and Gene Expression

Pearson correlation coefficients between serum levels of reproductive hormones, stress-related biomarkers, and the expression levels of the candidate genes IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1, are presented in Table 6. Expression of EXOC4 and GPR151 showed moderate associations with serum levels of AMH, progesterone, cortisol, and HSP70 (P < 0.05). TRHDE1 was moderately correlated with AMH and progesterone (P < 0.05) but exhibited low correlations with cortisol and HSP70 (P < 0.05). GLI3 showed moderate correlations with progesterone, cortisol, and HSP70, whereas IQGAP1 was not significantly correlated with the reproductive hormones or with the stress-related biomarkers cortisol or HSP70 (P > 0.05).

4. Discussion

Genome-wide association studies (GWAS) and marker-assisted selection (MAS) were sequentially applied in the current study to identify SNPs from candidate genes as potential molecular markers for fertility and thermotolerance traits in Holstein cows exposed to HS. Following GWAS analysis, 21 and 20 SNPs were identified as significant predictors (P < 1.09 × 10⁻⁶) of days open (DO) and services per conception (SPC), respectively. Five SNPs were common predictors for both reproductive traits, and after meeting the established selection criteria were genotyped in an independent Holstein population using a TaqMan assay and qPCR. A genotype to phenotype association analysis validated four of these SNPs as molecular markers for fertility and thermo-tolerance, confirming their favorable additive effects on the evaluated traits (i.e., DO, SPC, RT, and RR). Additionally, cows carrying favorable SNP genotypes exhibited higher mRNA expression of the four candidate genes, along with corresponding reproductive and thermotolerance-related hormonal profiles. In summary, these findings confirm the utility of integrating genomic technologies, together with gene expression analyses, and hormonal profiling, to functionally validate genetic markers that may enhance fertility and heat-stress resilience in dairy cattle breeding programs.

Semiarid regions are characterized by intense summer heat, with elevated temperatures and humidity, leading to HS in animals. Heat stress refers to environmental conditions that disrupt the balance between heat gain and an animal’s capacity to dissipate excess body heat [25]. Dairy cattle are particularly susceptible to HS, as the metabolic demands of milk production increase internal heat load and further challenge thermoregulatory capacity [26]. Extensive research has focused on determining the onset of HS in dairy cows, demonstrating that exposure to elevated environmental temperatures triggers physiological and behavioral adaptations intended to preserve homeothermy [27,28]. These adaptive responses, including increased respiratory rate, elevated heart rate, reduced feed intake, greater water consumption, and a behavioral tendency to seek shaded areas, often compromise fertility and overall productive performance [4,29,30].

According to the environmental data analyzed, THI values confirmed that Holstein cows in the current study were exposed to increasing levels of moderate to severe HS (70-82 units) during summer reproductive management. Identifying HS in dairy cows can be challenging, therefore, behavioral and physiological indicators are commonly used to detect it, although their accurate assessment often requires trained personnel [3]. To enable a more objective evaluation, environmental indices such as the Temperature–Humidity Index (THI) have been developed and are widely applied [31]. THI shows a strong correlation with physiological responses, including respiration rate and body temperature, and is commonly used to assess the presence and severity of HS [32]. While cows can tolerate a certain THI range without major productivity losses, exceeding critical thresholds disrupts thermal balance, resulting in reduced milk yield and impaired reproductive performance [29,33].

Poor female fertility is a major global challenge in the modern dairy industry; largely due to unfavorable genetic correlations with milk yield traits that have been intensively selected over recent decades [34]. This slow improvement is additionally attributable to the low heritability of fertility traits and the fact that they are typically measured later in an animal’s productive life [35]. In addition, the climatic conditions in semi-arid regions that lead to HS further exacerbate reproductive inefficiency by disrupting endocrine function and reducing food intake, thereby intensifying the negative impact of high milk yield on fertility [12]. Faster genetic improvement in dairy cow fertility can be achieved by integrating quantitative trait loci (QTL) information into selection programs. Genome-wide association studies, enabled by high-density SNP panels, have been widely used since 2005 to identify genetic variants associated with female fertility across diverse Holstein populations [36]. However, the identification of genetic markers associated with high fertility in Holstein cows under intense HS remains scarce.

In the current study, single-marker GWAS analysis detected 21 and 20 SNPs associated with the reproductive traits DO and SPC, respectively. Five of these SNPs were common to both genomic analyses and met multiple-testing selection criteria, warranting their designation as genomic candidate SNPs for fertility traits in Holstein cows managed under HS environmental conditions. These SNPs were intron variants within the genes IQGAP1, EXOC4, GLI3, GPR151, and TRHDE1. Although located within intronic or other non-coding regions, these SNPs may still affect the expression of their associated genes [36]. Such effects are likely mediated by cis-regulatory elements, including transcription factor binding sites, enhancers, silencers, and insulators, which modulate gene activity. Furthermore, long non-coding RNAs residing in intronic regions may regulate protein expression through epigenetic, transcriptional, and post-transcriptional mechanisms [37,38].

Fertility in mammals is a complex trait shaped by multiple genetic loci and substantial environmental variation, making the identification of causal mutations particularly challenging [39]. In cattle, although numerous significant associations have been reported, limited overlap across populations highlights the difficulty of identifying consistent candidate genes [40]. Recently, multi-trait GWAS of several fertility traits revealed significant overlaps between cattle and human fertility-related genes, suggesting the presence of conserved genomic regions regulating fertility across mammals [41]. Following this rationale, we assume that detecting the same genetic variants across independent GWAS for reproductive traits may provide stronger evidence of shared functional genomic regions, thereby enabling inference of common genomic mechanisms underlying fertility traits in dairy cows exposed to HS.

GWAS for fertility traits performed in Canadian Holstein cattle detected eight genome-wide significant SNPs (1% FDR) on BTA21 associated with DO, located within the 53–59 Mb interval, a significant chromosomal region that overlaps the FAM181A gene [16]. In a similar GWAS analysis, the genes IL6R, SLC39A12, CACNB2, ZEB1, ZMIZ1, and FAM213A were reported as strong candidates for several female fertility traits in Holstein cows, including DO [42]. Recently, a longitudinal GWAS in Holstein cows analyzing DO and fertility records related to SPC identified the genes TMEM132C, IMPG1, DCHS2, CSMD1, and CSNK1A1 as potential fertility predictors [43]. GWAS of pregnancy and conception rates in Holstein heifers and cows identified genes known to affect reproduction, including GNRHR, SHBG, and ESR1 [44]. Likewise, a genomic study of fertility in Holstein cows revealed overlapping QTL on chromosomes 6 and 29, as well as positional candidate genes associated with oogenesis, pregnancy rate, conception rate, and reproductive longevity [45]

The five genomic SNPs selected in the current study met the candidate marker selection criteria (i.e., MAF > 0.10 and HWE > 0.05) and were considered appropriate for a SNP association study performed in an independent set of Holstein cows. Although GWAS statistical models incorporate multiple-testing correction procedures (Bonferroni adjustment) and false discovery rate (FDR) control to reduce false positives, these stringent thresholds may limit the detection of polygenic effects [46,47,48]. Therefore, despite applying these corrections, validating significant SNPs in independent animal populations remains the most reliable approach to confirm their true biological relevance, given the low probability of observing consistent associations across populations [49,50].

After our validation study, the SNPs rs136074864, rs109479159, rs43520457, and rs133770320 in the genes EXOC4, GLI3, GPR151, and TRHDE were confirmed as significant marker predictors for fertility (FOL, CL, SPC, DO) and thermo-tolerance traits (RT, RR). These results provide strong evidence of a genetic link between reproductive efficiency and heat-stress resilience in lactating Holstein cows, highlighting the potential role of pleiotropic mechanisms underlying fertility and thermo-tolerance phenotypes.

The exocyst complex component 4 (EXOC4) gene encodes a component of the exocyst complex, an eight-protein assembly (EXOC1–EXOC8) that directs exocytic vesicles to specific docking sites on the plasma membrane. It plays a key role in essential cellular processes, including cell growth, division, and polarization, while integrating signaling pathways and coordinating protein assembly and function [51]. EXOC4 protein is expressed in the placental syncytiotrophoblast, and its regulatory role in trophoblast differentiation and placental peptide hormone secretion may indicate some effect on reproductive function [52]. Additionally, EXOC4 showed a significant association with circulating triiodothyronine (T3) and thyroxine (T4) concentrations in Chinese Holstein cattle, providing insight into its potential functional role in bovine physiology [53]. A SNP within the EXOC4 gene was reported as associated with gestation length in Angus, Charolais, and Holstein-Friesian cattle populations [54]. Similarly, the EXOC4 gene was also reported as a candidate marker for age at first calving in Brazilian beef cattle [55] and calving-related fertility traits in Duroc sows [56].

In the current study, we hypothesized that the SNP rs136074864 in the EXOC4 gene might be associated with fertility in Holstein cows exposed to HS, given its central role in essential cellular processes and reproductive function. As a component of the exocyst complex, it is involved in vesicle trafficking, cell growth, and polarization, which are processes critical for embryonic development and placental function. Its expression and regulatory role in placental tissues and endocrine patterns suggest a direct influence on gestation maintenance. Moreover, its association with thyroid hormones (T3 and T4), which are key regulators of metabolism and thermal stress response, as well as with traits such as gestation length and age at first calving, further supports its relevance to reproductive efficiency under HS conditions.

The GLI Family Zinc Finger 3 (GLI3) gene encodes for a protein transcription factor considered as a major component of the Hedgehog signaling pathway. After phosphorylation and nuclear translocation, GLI3 can act either as an activator or as a repressor of the Hedgehog pathway [57]. This pathway is composed of 3 molecules: sonic (SHH), desert (DHH), and indian (IHH). SHH primarily regulates central nervous system development, DHH mainly regulates gonadal development, and IHH contributes to bone and joint formation [58], however, they are able to exhibit overlapping roles. DHH and IHH are expressed in ovarian granulosa with central role in regulating steroidogenesis and follicular development [59]. The HH pathway has been also associated with HS tolerance in sheep, together with glycerolipid and thyroid hormones metabolism, as these signaling mechanisms appeared to play an important role in mitigating the detrimental effects of thermal stress and enhancing adaptive thermotolerant responses [60].

In our study, we assumed that the SNP rs109479159 in the gene GLI3 is a potential candidate marker for fertility and thermotolerance because it regulates key biological processes involved in reproductive function and cellular adaptation to HS. Its role within the Hedgehog signaling pathway links GLI3 to mechanisms controlling follicular development, steroidogenesis, and gonadal function, while emerging evidence also associates this pathway with adaptive responses to HS. Therefore, genetic variation in GLI3 may contribute to individual differences in reproductive efficiency and resilience under thermal stress conditions, supporting its potential application as a genetic marker in dairy cattle breeding programs aimed at improving both fertility and thermotolerance.

The G Protein-Coupled Receptor 151 (GPR151) gene encodes an orphan receptor belonging to the class A rhodopsin-like group of G protein–coupled receptors (GPCRs), alongside galanin receptors (GalR1, GalR2, GalR3) and the kisspeptin receptor (Kiss1R- GPR54) [61,62]. Within this family, GPR151 is classified in the SOG subfamily, which comprises somatostatin, opioid, galanin, and kisspeptin receptors [63]. Because of its high sequence similarity to galanin receptors, it was initially proposed that GPR151 could be activated by the neuropeptides galanin or kisspeptin [64]. Galanin exerts diverse physiological effects on the neuroendocrine axis, including regulation of food intake, insulin secretion, and somatostatin release [65], whereas kisspeptins are essential regulators of the reproductive axis in humans and animals, as they stimulate GnRH secretion [66,67]. Additionally, GPR151 has been implicated in adipocyte differentiation and hepatic gluconeogenesis, processes that are considered essential components of the heat-stress response because of their key roles in energy metabolism [68,69].

Therefore, there is a strong reason to propose the SNP rs43520457 in the GPR151 gene as marker for fertility in cows exposed to HS, given its functional relationship with neuroendocrine and metabolic pathways. This gene is closely associated with galanin and kisspeptin receptors, both of which are involved in regulating the hypothalamic–pituitary–gonadal (HPG) axis. Under HS, cows reduce feed intake and mobilize adipose reserves as an energy source, leading to profound metabolic alterations that can disrupt reproductive hormonal signaling. Given that GPR151 has been implicated in adipocyte differentiation and hepatic gluconeogenesis, it may help integrate metabolic and neuroendocrine signals during negative energy balance. In this context, genetic variants in GPR151 could influence metabolic adaptation to thermal stress and, consequently, ovarian function and fertility in Holstein cows exposed to HS.

Thyrotropin-releasing hormone-degrading enzyme (TRHDE) gene encoding a peptidase belonging to the M1 family that regulates the availability and activity of thyrotropin-releasing hormone (TRH), a neuropeptide synthesized primarily in the hypothalamus [70]. TRH plays important roles in endocrine regulation by stimulating the secretion of thyroid-stimulating hormone (TSH) and prolactin (PRL), and therefore, TRHDE may influence processes related to metabolism, lactation, and reproduction by regulating TRH availability [71]. PRL, in addition to its essential role in the initiation and maintenance of lactation, is involved in the regulation of reproductive function and is a hormone responsive to stressful conditions [72]. In this regard, alterations in plasma PRL concentrations have been associated with heat stress in dairy cattle [73], while thermal stress can negatively affect the hypothalamic-pituitary-gonadal axis and consequently impair follicular dynamics and fertility [74].

In the current study, the association of the SNP rs133770320 in the TRHDE gene with fertility and thermotolerance in dairy cattle suggests that genetic variation in this gene may contribute to individual differences in the physiological response to heat stress. Given the role of TRHDE in endocrine regulation, the identified association may reflect differences in the ability of animals to maintain endocrine homeostasis under thermal challenges. Such regulation could influence reproductive function and, consequently, fertility in heat-stressed dairy cattle. Therefore, rs133770320 represents a promising candidate variant for further investigation of the genetic mechanisms underlying the interaction between HS, endocrine regulation, and reproductive performance.

The increased expression observed in our study in cows with favorable genotypes for the genes EXOC4, GLI3, GPR151, and TRHDE1 may be explained by their central role in integrating reproductive and metabolic processes under HS conditions. This association is supported by their involvement in essential cellular processes, including adipocyte differentiation, hepatic gluconeogenesis, and overall energy metabolism. During HS, reduced feed intake leads cows to enter a negative energy balance, which activates pathways that mobilize body reserves and optimize nutrient utilization. These coordinated metabolic responses help sustain cellular function, maintain endocrine balance, and ensure adequate nutrient availability, thereby supporting ovarian activity and reproductive performance under thermal stress.

Increased expression of the EXOC4 gene was reported in sows with SNPs favorably associated with the number of piglets born (NW) and litter weight at birth (LWW), which was due to the presence of variants located in regulatory regions of the promoter. This interaction likely enhances the gene’s transcriptional activity, leading to higher expression levels [56]. A GWAS analysis in Blackbelly sheep detected a chromosomal region in BTA5 harboring several SNPs within the GPR151 gene whose expression was associated with coat phenotypes [75]. Expression of GLI3 gene influences female fertility through its role in the developmental regulation of the hypothalamic–pituitary–gonadal axis. During embryonic development, GLI3 is required for the proper formation of olfactory unsheathing cells, which support the migration of gonadotropin-releasing hormone (GnRH-1) neurons from the vomeronasal region to the brain [76].

Gene expression studies are an essential complement to genome-wide association studies (GWAS) because they provide functional insights into the biological mechanisms underlying statistically associated loci. Expression analyses quantify transcribed RNA molecules that can be translated into proteins that contribute to phenotypic variation. While GWAS identify genomic regions or SNPs linked to specific traits, they do not necessarily reveal which gene is causative or how it influences the phenotype. Gene expression analysis helps determine whether candidate genes within associated regions are actively expressed and capable of influencing the phenotype [77,78].

In the current study, serum levels of AMH, progesterone, and HSP70 were higher in cows carrying favorable SNP genotypes for the genes EXOC4, GLI3, GPR151, and TRHDE1. In contrast, serum levels of cortisol, an important marker stress response, were significantly reduced in cows with favorable genotypes. These findings suggest that favorable genetic variants not only enhance reproductive endocrine function but also mitigate physiological stress responses, supporting improved fertility and greater adaptation to HS conditions. Interestingly, correlations observed between serum hormone profiles and gene expressions reinforce the functional role of these genes as key molecular links between reproductive efficiency and thermo-tolerance.

Heat stress in Holstein cows activates the hypothalamic–pituitary–adrenal (HPA) axis, the primary endocrine regulator of the stress response, thereby increasing cortisol secretion [79,80]. Cortisol is widely used as a biomarker of HS and can stimulate the production of heat shock proteins (HSPs). Among them, HSP70 is particularly sensitive to HS in ruminants, with expression levels rising in proportion to stress severity and the individual cow’s capacity for thermal stress response [81,82]. HSP70 plays a protective cellular role by stabilizing the cytoskeleton, regulating the cell cycle, modulating immune responses, and preventing apoptosis [83]. However, elevated cortisol can suppress luteinizing hormone secretion, negatively affecting reproductive function [84].

Therefore, cortisol and HSP70 are important biomarkers of HS that significantly influence the reproductive performance of Holstein cows. Interestingly, the negative correlation observed in our study between cortisol and the expression levels of the genes EXOC4, GLI3, GPR151, and TRHDE1, together with the positive association of these genes with reproductive hormones and HSP70, confirms their key regulatory role in physiological mechanisms that enable Holstein cows exposed to HS to maintain adequate reproductive performance. These findings further suggest that enhanced expression of these genes may contribute to improved endocrine balance and metabolic adaptation, thereby promoting resilience to thermal stress while sustaining fertility.

One limitation of this study is the relatively small number of Holstein cows included in the GWAS analysis, which may reduce statistical power and limit the detection of additional loci associated with fertility and thermotolerance traits. However, validation of the identified genomic markers in an independent Holstein population, together with functional gene expression analyses and measurement of serum hormone levels indicative of reproductive status and thermal balance, provided robust evidence supporting the key regulatory role of the genes EXOC4, GLI3, GPR151, and TRHDE1 in modulating fertility and thermotolerance in Holstein cows managed under conditions of intense heat stress.

5. Conclusions

Integration of GWAS and marker-assisted molecular technologies, further complemented by gene expression profiling and hormonal analyses, provided comprehensive evidence that SNPs within the genes EXOC4, GLI3, GPR151, and TRHDE1 play a central role in regulating fertility and thermo-tolerance in Holstein cows under HS conditions. These findings highlighted the functional relevance of these genes as molecular markers for reproductive efficiency and HS resilience. Collectively, this study contributes valuable insights for genetic selection strategies aimed at improving productivity and sustainability in dairy cattle exposed to challenging thermal environments. Nevertheless, further studies, including a larger number of Holstein cows, are necessary to identify additional genes that explain the variability observed in fertility under HS conditions.