Genetic dissection of reproductive performance of dairy cows under heat stress

Summary Heat stress negatively impacts the reproductive performance of dairy cows. The main objective of this study was to dissect the genetic basis underlying dairy cow fertility under heat stress conditions. Our first goal was to estimate genetic components of cow conception across lactations considering heat stress. Our second goal was to reveal individual genes and functional gene‐sets that explain a cow’s ability to conceive under thermal stress. Data consisted of 74 221 insemination records on 13 704 Holstein cows. Multitrait linear repeatability test‐day models with random regressions on a function of temperature–humidity index values were used for the analyses. Heritability estimates for cow conception under heat stress were around 2–3%, whereas genetic correlations between general and thermotolerance additive genetic effects were negative and ranged between −0.35 and −0.82, indicating an unfavorable relationship between cows’ ability to conceive under thermo‐neutral vs. thermo‐stress conditions. Whole‐genome scans identified at least six genomic regions on BTA1, BTA10, BTA11, BTA17, BTA21 and BTA23 associated with conception under thermal stress. These regions harbor candidate genes such as BRWD1, EXD2, ADAM20, EPAS1, TAOK3, and NOS1, which are directly implicated in reproductive functions and cellular response to heat stress. The gene‐set enrichment analysis revealed functional terms related to fertilization, developmental biology, heat shock proteins and oxidative stress, among others. Overall, our findings contribute to a better understanding of the genetics underlying the reproductive performance of dairy cattle under heat stress conditions and point out novel genomic strategies for improving thermotolerance and fertility via marker‐assisted breeding.


Introduction
Heat stress negatively impacts dairy cattle performance. It is now known that the intense selection for productive traits in recent decades has compromised the thermoregulatory competence of dairy cows (Aguilar et al. 2009;Santana et al. 2017;Sigdel et al. 2019). In particular, heat stress is a major cause of low fertility in dairy cattle. Indeed, heat stress disrupts several reproductive processes, including follicular growth, steroid production, estrous expression, oocyte competence, uterine endometrial response and embryonic growth, leading to conception failure, early embryonic mortality and pregnancy loss (Jordan 2003;Hansen 2009;Silva et al. 2013). Several studies have reported that there is a marked decrease in pregnancy rate during the summer (Hansen & Arechiga 1999;Sartori et al. 2002;de Vries & Risco 2005). Summer heat stress affects reproductive performance because cows undergo adaptive mechanisms during periods of thermal stress by redirecting blood flow from core to periphery, thereby impacting cyclicity, pregnancy establishment and fetal development (Al-Katanani et al. 2002;Hansen 2009). In dairy herds in the southeast USA, there is a strong seasonal variation in conception per insemination. For instance, de Vries & Risco (2005) investigated the seasonality of annual pregnancy rate (between 71 and 364 days since last calving) in dairy herds in Florida and Georgia, reporting a 15.8% pregnancy rate during winter compared with 5.6% during summer. This decline in pregnancy rate during periods of heat stress has a huge economic cost to the dairy industry (De Vries 2006). Various strategies, including cooling, shading and nutrition, have been used to alleviate the negative effects of thermal stress on dairy cow fertility. Nonetheless, there is a clear decline in conception in hot and tropical climates (Schuller et al. 2014;Hagiya et al. 2017). In this context, genetic selection for improved thermotolerance is an attractive complementary approach for reducing the effects of heat stress and subsequently improving reproductive performance. Ravagnolo & Misztal (2000) presented a methodology for genetic evaluation of productive and reproductive traits under heat stress by combining test day records with public weather station records. The statistical model has two genetic effects, a regular genetic effect corresponding to performance under thermo-neutral conditions and a heat-stress effect corresponding to the rate of decline under thermal stress. Heat stress is measured using a temperature-humidity index (THI, NOAA 1976). In this context, a linear regression of phenotypic records on an environmental variable (THI) is fitted with the assumption that the phenotype is unaffected until a certain threshold level of THI, and above that level, the phenotypic performance declines linearly with increasing THI. Genetic variation is associated with the rate of decline under thermal stress conditions. Using this methodology, Ravagnolo & Misztal (2002) evaluated the effect of heat stress on non-return rate in US Holstein cows. Recently, using a similar approach, Ansari-Mahyari et al. (2019) investigated the genetic variability of conception rate and also non-return rate in Iranian dairy cows under heat stress conditions. Both studies indicated that there is a substantial genetic variation underlying conception rate under heat stress, and hence, dairy cow reproductive performance under thermal stress conditions could be genetically improved.
There is increasing knowledge about genes and pathways associated with conception in dairy cattle. For instance, Hoglund et al. (2015) identified genomic regions on BTA20 and BTA23 significantly associated with pregnancy rate in Danish Jersey cattle. Liu et al. (2017) identified individual genes on BTA8, BTA13 and BTA23, such as PRPF4B and PXDC1, as potential candidate genes affecting pregnancy establishment and maintenance in Chinese Holstein cattle. Recently, Kiser et al. (2019) reported a list of loci, positional candidate genes, including ARVCF, GJB4, GJB5 and JMY, and transcription factor binding sites associated with conception at first service and number of services per conception in US primiparous Holstein cows. On the other hand, there is no information regarding individual genes, biological pathways and molecular mechanisms affecting dairy cows' ability to successfully conceive under thermal stress. As such, the first objective of this study was to estimate genetic components for conception per insemination across lactations considering heat stress using random regressions as a function of THI values. The second objective of this study was to perform whole-genome scans and subsequent gene-set enrichment analyses in order to identify genomic regions, individual genes and functional gene-sets implicated in cow conception under heat stress. This relevant information could contribute to improving dairy cows' ability to conceive under heat stress using novel genomic strategies.

Phenotypic and genotypic data
Conception per insemination is an important fertility trait in dairy cattle, which is defined as the cow's ability to establish and maintain a successful pregnancy, provided the cow is inseminated at the time of ovulation (Averill et al. 2004). Phenotypic data consisted of 74 221 insemination records on 13 704 Holstein cows calved from 2006 to 2016 in two dairy herds in the State of Florida, USA. Service data includes 28 400 insemination records for first-lactation cows, 26 809 records for second-lactation cows and 19 012 records for third-lactation cows. A total of 6838 cows had insemination records in at least two consecutive lactations. All insemination records until 400 days in milk (DIM) were used in the analyses. The outcome of each insemination was recorded as 0 and 1, where 1 indicates that insemination achieved a successful conception and 0 indicates otherwise. The pedigree file included 28 620 animals based on a five-generation pedigree.
Genotype data for 60 671 SNP markers were available for a total of 4700 cows with insemination records and also 1592 sires in the pedigree. The SNP data were kindly provided by the Cooperative Dairy DNA Repository and the Council on Dairy Cattle Breeding. Those SNP markers that mapped to sex chromosomes, were monomorphic or had minor allele frequency less than 1% were removed from the genotype data. After quality control, a total of 58 046 SNPs were retained for whole-genome scans and subsequent gene-set analyses.

Weather information
Weather data were obtained from Florida Automated Weather Network for Alachua county (https://fawn.ifas.ufl.edu/) and hourly THI was calculated following Bohmanova et al. (2007) as THI ¼ ð1:8 temp þ 32Þ À ð0:55 À 0:0055 rhÞð1:8temp À26Þ, where 'temp' is the temperature in degrees Celsius and 'rh' is the relative humidity as a percentage. Mean daily THI of day of insemination was assigned to each insemination record as suggested by Ravagnolo & Misztal (2002).
A function of THI, denoted as f(THI), was created as a dummy variable in order to estimate decline in reproductive performance under heat stress conditions: where THI thr is equal to 68 and thus f(THI) is equal to max (0, THI À 68).

Statistical model
The following multitrait linear repeatability test-day model was used for all the analyses, considering the first three parities as different traits as in Aguilar et al. (2009): where y klmn is the outcome of the insemination (binary trait, 0 or 1), HTD kl is the fixed effect of herd-day k of insemination within parity l (l = 1, 2, 3), DIM m is the mth DIM class of insemination with classes defined every 20 days, a nl is the general random additive genetic effect (intercept) of animal n in parity l, pe nl is the general random permanent environmental effect (intercept) of cow n in parity l, f(THI) is a function of THI for herd-day k of insemination, v nl is the random regression additive genetic effect (slope) of conception per insemination per degree of THI above the threshold for the animal n in parity l (thermotolerance), q nl is the random regression permanent environmental effect (slope) of thermotolerance of the cow n in parity l and e klmn is the random residual effect.
nl be a vector of random additive genetic effects and p ¼ ½pe 0 nl; q 0 nl be a vector of random permanent environmental effects for lactations n = 1-3. The (co)variance structure was: where A is the numerator relationship matrix, Φ and w are 6 9 6 (co)variance matrices of random regression coefficients for additive and permanent environment effects respectively (three traits with two parameters, i.e. intercept and slope, per trait), R is a 3 9 3 diagonal matrix of residual variances corresponding to each trait and denotes the Kronecker product of matrices.

Genetic parameter estimation
Variance components were estimated using multitrait repeatability test-day models in a Bayesian framework using the software GIBBS2F90. Genotype data were not used for the estimation of variance components. Of a total of 500 000 samples, the first 100 000 samples were discarded as burn-in, and every 100th sample was retained to calculate features of the posterior distributions, such as posterior means and standard deviations. Convergence diagnostics of Markov chain Monte Carlo sampling output were performed by visual inspection of trace plots. Heritability (h 2 ) for conception per insemination at heat stress level f(i) was calculated as: where r 2 a is the variance of general additive genetic effects, r 2 v is the variance of thermotolerance additive genetic effects, r av is the additive genetic covariance between general and thermotolerance effects, r 2 pe is the variance of general environmental permanent effects, r 2 q is the variance of thermotolerance environmental permanent effects, r pq is the environmental permanent covariance between general and thermotolerance effects, f(i) is a function of THI and r 2 e is the residual variance. The genetic correlation between general and thermotolerance additive effects for conception per insemination was estimated as:

Gene mapping
The identification of genomic regions and individual genes affecting cow conception was performed using single-step genomic BLUP (ssGBLUP). In ssGBLUP, the statistical model is the same as the classic BLUP, but the inverse of the pedigree relationship matrix (A À1 ) is replaced with inverse of the realized relationship matrix (H À1 ) that combines both pedigree and genomic information (Aguilar et al. 2010). The combined pedigree genomic relationship matrix H À1 was calculated as follows: where G À1 is the inverse of the genomic relationship matrix and A À1 22 is the inverse of the pedigree relationship matrix of the animals with genotype data. In this study, G À1 has the dimension of 6292 9 6292 which includes 4700 cows with insemination records and 1592 sires in the pedigree. The A matrix has a dimension of 28 620 9 28 620, which was calculated based on a five-generation pedigree.
Candidate regions and individual genes associated with cow conception under thermo-neutral and thermo-stress conditions were identified based on the amount of genetic variance explained by 2.0 Mb windows of adjacent SNPs. Given the vector of genomic estimated breeding values, the SNP effects can be estimated as where b s is the vector of SNP marker effects, D is a diagonal matrix of weight of SNPs (in this study, D was set as an identity matrix), Z is a matrix relating genotypes of each SNP marker to observations and b u is the vector of genomic estimated breeding values (Wang et al. 2012). The percentage of genetic variance explained by a given 2.0 Mb window of adjacent SNPs was then calculated as: where u i is the genetic value of the ith genomic region under consideration, B is the total number of adjacent SNPs within 2.0 Mb region, z j is the genotype code of jth marker and s j is the marker effect of the jth SNP within the ith region. All of these calculations were performed using POSTGSF90 of the BLUPF90 family of programs (Aguilar et al. 2014).

Gene-set analysis
Gene-set enrichment analysis, also known as overrepresentation or pathway analysis, is a powerful tool that allows us to reveal biological pathways and molecular mechanisms underlying phenotypes of interest. Here, we conducted alternative gene-set analyses using different annotation databases, including GO, Medical Subject Headings, InterPro and Reactome, in order to obtain additional insights regarding biological processes that might explain cows' ability to conceive under thermal stress conditions. The gene-set analysis consisted basically of three steps (Han & Peñagaricano 2016): first, the assignment of SNP markers to annotated genes; then the assignment of genes to functional gene-sets or pathways; and finally the association between each functional term and the phenotype of interest. The UMD 3.1 bovine genome sequence assembly was used for SNP assignment using BIOCONDUCTOR R package BIOMART (Durinck et al. 2009). Here, SNPs were assigned to genes if they were located within the genomic sequence of the gene or at most 15 kb either upstream or downstream of the gene. The distance of 15 kb was chosen to capture regulatory regions that may lie close to but outside of the gene. An arbitrary threshold of 5% of the SNP effects distribution (in absolute value) was used to define relevant SNP markers, and potential thermotolerant genes were defined as those genes flagged by at least one relevant SNP. Finally, the identification of gene-sets or pathways significantly enriched with thermotolerant genes was performed using a Fisher's exact test, a test of proportions based on the cumulative hypergeometric distribution (Peñagaricano et al. 2012).

Genetic parameter estimation
Variance components for conception per insemination under thermo-neutral (intercept) and thermo-stress (slope) conditions were estimated using multitrait linear repeatability test day models (Table 1). Relevant genetic parameters include heritability and genetic correlations at THI = 78, i.e. moderate heat stress level, across the first three parities. Additive genetic variances under heat stress conditions increased between 80 and 200% across the first three parities, suggesting that heat stress has greater effects on conception rate in later lactations. Similarly, Biffani et al. (2016) reported that additive genetic variances for 56-day non-return rate under heat stress conditions increased from first to second parity in Italian Holstein cows. The increase in thermotolerance additive genetic variances across parities suggests that multiparous cows are more susceptible than primiparous cows to the negative effects of heat stress on reproductive performance.
Heritability estimates for conception per insemination at THI = 78, i.e. f(THI) = 10, were between 2 and 3% across the first three parities. Our findings are similar to those reported by Ravagnolo & Misztal (2002), who estimated a heritability equal to 1.4% for non-return rate 60 days after insemination at THI = 70 in primiparous US Holstein cows. Similarly, Hagiya et al. (2017) reported heritability estimates for conception rate between 1 and 3% in Holstein cattle in Japan under both mild and moderate heat stress conditions. Biffani et al. (2016) reported heritability estimates for 56-day non-return rate under heat stress of between 2.5 and 7.6% for the first two parities in Italian Holstein cows. Overall, all these studies suggest that there is a significant genetic variability underlying cow conception under heat stress, and therefore, genetic selection for improved fertility under thermal stress conditions is feasible.
Genetic correlations between general and thermotolerance additive genetic effects were always negative and ranged from À0.35 to À0.82 across the first three parities. These negative genetic correlations suggest that the continued selection for greater cow fertility under thermoneutral conditions will result in increasing even more the harmful effects of heat stress on reproductive performance. Our results are in agreement with the findings of several previous studies. For instance, for 56-day non-return rate in Italian Holstein cows, genetic correlations between general and thermotolerance stress effects were À0.31 for firstparity cows and À0.45 for second-parity cows (Biffani et al. 2016). Similarly, first-parity US Holstein cows showed a marked negative genetic correlation (À0.77) for 60-day non-return rate between regular and heat tolerance effects (Ravagnolo & Misztal 2002). The negative relationship between cows' ability to conceive under thermo-neutral vs. thermo-stress conditions emphasize the importance of selecting animals in the same environment under which they are going to reproduce. proportion of the additive genetic variance explained by 2.0 Mb SNP windows. Left plots report the genomic regions associated with cow conception under thermo-neutral conditions, whereas the right plots report the genomic regions affecting cow conception under heat stress.

Gene mapping
Genomic scans under thermo-neutral conditions clearly identified a peak on BTA23 at 3.60-5.60 Mb that explains more than 1.5% of additive genetic variance across the three parities. This genomic region harbors at least two putative genes, BMP5 and HCRTR2, that are directly implicated in early embryo development. Gene BMP5 encodes a member of the transforming growth factor-b superfamily and is highly expressed in placenta. Interestingly, previous research has shown that early addition of BMP5 to the embryo culture medium had a positive effect on the blastocyst rate, affecting the expression of several BMP target and pluripotency genes, suggesting that BMP5 plays an important role in the preimplantation development of bovine embryos (Garcia et al. 2015). Gene HCRTR2 encodes a G-protein coupled receptor involved in the regulation of feeding behavior and energy metabolism that appears to be directly implicated in uterine processes during early pregnancy (Smolinska et al. 2017). Another 2.0 Mb SNP window located on BTA12 at 13.48-13.66 Mb also explained a significant amount of general genetic variance across all three lactations. This region harbors at least two genes, LCP1 and TPT1, which are related to embryonic development. Gene LCP1 is involved in gastrulation and embryonic epidermal development (Baumgartner et al. 2019). Gene TPT1 is expressed in placenta and is implicated in calcium binding and homeostasis of trophoblast cells, which is important for growth and development of the fetus (Arcuri et al. 2005).
The ssGBLUP method also revealed several genomic regions and putative genes that partly explain the observed variation in cow conception under heat stress conditions (Fig. 1, right plots). Interestingly, most of the regions that explained a large amount of thermotolerance additive genetic variance are parity specific. For primiparous cows, the two most important regions are located on  and . The genomic region on BTA10 harbors at least two genes, EXD2 and SLC10A1, that are involved in cellular response to heat stress. Gene EXD2 is implicated in the cellular response to oxidative stress owing to heat stress; it protects the cells from oxidative stress by reducing reactive oxygen species (ROS) production. Indeed, the KO of EXD2 results in oxidative stress and cell death (Zlotorynski 2018). Gene SLC10A1 is involved in cholesterol biosynthesis and metabolism, and cholesterol plays an important role in cellular protection against heat stress through the activation of heat shock proteins production (Balogh et al. 2013). In addition, gene ADAM20, also located on BTA10 at 80.75-82.75 Mb, encodes a membrane-anchored protein involved in cell-cell and cell-matrix interactions, including fertilization and early embryonic development (Sha et al. 2018). Moreover, the genomic region on BTA21 at 21.31-23.31 Mb harbors at least two relevant genes, FSD2 and AP3B2. Gene FSD2 is an oxidative stress-related gene which acts as an ROS scavenger, protecting the cells from oxidative damage caused by heat stress (Myouga et al. 2008). Gene AP3B2 is implicated in the intracellular removal of heat-induced protein aggregates, maintaining protein homeostasis (Wilbe et al. 2015).
For second-lactation cows, two different regions on BTA11 (27.18-29.18 Mb) and BTA17 (56.50-58.50 Mb) explained large amounts of additive genetic variance for conception per insemination under heat stress. The genomic region on BTA11, which explained almost 1.5% of thermotolerance genetic variance, harbors gene EPAS1, which is an oxidative stress-related gene, directly implicated the hypoxic response, a common condition during thermal stress. The KO of EPAS1 in mice results in early embryonic mortality because of hypoxia (Xu et al. 2014). Notably, gene EPAS1 has already been associated with both fertility traits and rectal temperature/thermotolerance in dairy cattle HPD, Highest posterior density; r 2 a , general additive genetic variance; r 2 e , residual variance; 100 r 2 v , thermotolerance additive genetic variance; 10 r av , additive genetic covariance between general and thermotolerance effect; r G a;v , genetic correlation between general and thermotolerance effect; cor ht , thermotolerance additive genetic correlation; and cor gen , general additive genetic correlation.
© 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 511-520 (Dikmen et al. 2015;Ortega et al. 2016). Moreover, the genomic region on BTA17 at 56.50-58.50 Mb harbors at least two genes, namely TAOK3 and NOS1, that are directly implicated in the cellular response to heat stress. In fact, gene TAOK3 is involved in the repair of DNA owing to different insults that damage DNA such as thermal stress. The knockdown of TAOK3 inhibits the cellular response to DNA damage, resulting in cell death (Raman et al. 2007). Gene NOS1 is implicated in the maintenance of cellular redox homeostasis owing to its ability to neutralize harmful ROS produced during heat stress. Gene NOS1 enhances the antioxidant enzyme activities of glutathione reductase and ascorbate peroxidase, which have a protective role during heat stress (Shi et al. 2014).
The genomic region located on BTA1 at 139.09-141.00 Mb explained almost 1% of the additive genetic variance for conception per insemination under heat stress for third parity. Notably, this region harbors gene BRWD1, Figure 1 Whole-genome scans for conception per insemination for the first three parities (numbered vertically as lactation 1, lactation 2 and lactation 3). The left plots highlight genomic regions and candidate genes affecting cow conception under thermo-neutral conditions (general additive genetic effects), whereas the right plots highlight genomic regions and putative genes implicated in cow conception under heat stress conditions (thermotolerance additive genetic effects).

Figure 2
Functional terms and pathways significantly enriched with genes associated with dairy cow conception under thermal-stress conditions. Four gene annotation databases were analyzed: GO, Medical Subject Headings, InterPro and Reactome. The y-axis displays the names and the total number of genes of each gene-set. The black dots represent the significance of enrichment (Àlog 10 P-value, Fisher's exact test, top x-axis) and the bars represent the number of significant genes in each functional term (bottom x-axis).
which is involved in meiotic chromosomal stability in the oocytes, and hence plays a crucial role in ensuring a successful conception (Pattabiraman et al. 2015).
Overall, the whole-genome scans identified several regions implicated in dairy cow conception under heat stress conditions. Interestingly, most of these genomic regions harbor candidate genes that are directly involved in either the cellular response to heat stress and/or fertilization and early embryo development. Table S1 provides information about putative functional variants that are polymorphic in Holstein cows and are located either within or close to (AE5 kb) these candidate genes. This information was kindly provided by the 1000 Bull Genomes Project (http://www.1000bullgenomes.com) based on whole-genome sequence data from 824 Holstein bulls.

Over-representation analysis
Of the 58 046 SNP markers evaluated in the whole-genome association mapping, 27 488 SNPs were located either within or near annotated genes. This set of SNPs defined a total of 19 305 genes annotated in the bovine reference genome. A subset of 928 of these 19 305 genes were flagged by at least one relevant SNP (top 5% of the SNP effects distribution) in at least two parities, and hence these 928 genes were defined as relevant genes for cow conception under heat stress conditions. Figure 2 shows the most relevant functional terms and biological pathways enriched with thermotolerant genes. Note that the over-representation analysis interrogated four different databases, GO, Medical Subject Headings, InterPro and Reactome.
At least five different groups of gene-sets were identified; these functional terms are related to fertilization, development, heat shock proteins, cellular response to oxidative stress, and calcium ion homeostasis. Table S2 reports the full list of significant GO terms, including GO ID, GO name, total number of genes, number of relevant genes and Fisher's P-value.
Noticeably, some of most significant terms are directly involved in cow conception and early embryo development, such as fertilization (D005306), sperm-egg recognition (GO:0035036), single fertilization (GO:0007388), trophoblast giant cell differentiation (GO:0060707) and developmental biology (R-BTA-1266738). These significant reproductive terms had in common at least six thermotolerant genes, namely SPESP1, PLCZ1, HAND1, FYN, DOCK1 and PI. All of these genes are involved in the process of fertilization and subsequent early embryo development (Halet et al. 2008;Fujihara et al. 2010).
Of special interest, many gene-sets that showed an overrepresentation of thermotolerant genes are directly associated with cellular response to heat stress, such as protein refolding (GO:0042026), HSP70 heat-shock proteins (D018840), heat shock protein 70 family (IPR013126), HSP40/DnaJ peptide-binding (IPR008971), HSC70 heat-shock proteins (D050883) and chaperone DnaJ, c-terminal (IPR002939). These heat shock response terms contain at least four relevant genes, namely HSF1, DNAJB4, CAMK2D and HSPA1A, all of them directly involved in cellular response to heat stress by maintaining proper protein folding, preserving cytoskeletal integrity and removing misfolded proteins during thermal stress (Collier et al. 2008;Gray & Heller Brown 2014;Min et al. 2015).
There are also some relevant functional terms related to oxidative stress, including response to oxidative stress (GO:0006979), glutathione transferase (D005982), ionotropic glutamate receptor (IPR001320) and ionotropic activity of kainate receptors (R-BTA-451306). Oxidative stress occurs as a consequence of an imbalance between ROS production and the antioxidant defense, and heat stress is an environmental factor responsible for stimulating ROS production. These terms were significantly enriched with genes GRIK2, PRDX1 and PSEN1, all directly implicated in cell survival during heat stress.
Calcium ion homeostasis and metabolism was also identified as a biological process significantly enriched with genes implicated in cow conception under heat stress conditions. Indeed, terms such as calcium ion binding (GO:0005509), calcium (D002118), EGF-like calcium binding domain (IPR001881) and activation of Ca-permeable kainite receptor (R-BTA-451308) were among the most significant ones. These calcium-related gene-sets contain at least three relevant genes, MYL2, NELL2 and PLCZ1, that are associated with fertilization and embryonic development.

Conclusions
In this study, we performed a comprehensive genetic and genomic analysis in order to reveal the genetic and biological basis of reproductive performance of dairy cows under heat stress conditions. Heritability estimates for conception per insemination under heat stress ranged between 2 and 3%, suggesting that the ability of a dairy cow to conceive under heat stress is influenced by genetic factors, and hence it could be improved by genetic means. Notably, the genetic correlations between general and thermotolerance additive genetic effects were negative and ranged from À0.35 to À0.82. These findings reinforce the idea that there is a negative relationship between cows' ability to conceive under thermo-neutral vs. thermo-stress conditions, and hence the continued selection for greater fertility ignoring heat tolerance will result in even greater increases in susceptibility to heat stress. Therefore, the dairy industry should focus on selecting cows that conceive well regardless of the level of heat stress. The whole-genome scans identified at least six regions located on BTA1, BTA10, BTA11, BTA17, BTA21 and BTA23 that explained a large amount of thermotolerance additive genetic variances. Interestingly, most of these regions harbor genes directly implicated in either fertilization and early embryo development or cellular response to heat stress. The geneset analysis revealed at least five relevant processes, namely fertilization and development, protein misfolding and heat shock proteins, cellular response to oxidative stress and calcium ion homeostasis, as significantly enriched with genes associated with cow conception under heat stress. Overall, this study contributes to a better, deeper understanding of the effect that heat stress has on cow reproduction. Our findings may also contribute to the development of novel breeding strategies to improve cow reproductive performance under heat stress.