Genome-wide associations for heat stress response suggest potential candidate genes underlying milk fatty acid composition in dairy cattle

Contents of milk fatty acids (FA) display remarkable alterations along climatic gradients. Detecting candidate genes underlying such alterations might be beneficial for the exploration of climate sensitivity in dairy cattle. Consequently, we aimed on the definition of FA heat stress indicators, considering FA breeding values in response to temperature-humidity index (THI) alterations. Indicators were used in GWAS, in ongoing gene annotations and for the estimation of chromosome-wide variance components. The phenotypic data set consisted of 39,600 test-day milk FA records from 5,757 first-lactation Holstein dairy cows kept in 16 large-scale German cooperator herds. The FA traits were C18:0, polyunsaturated fatty acids (PUFA), saturated fatty acids (SFA), and unsaturated fatty acids (UFA). After genotype quality control, 40,523 SNP markers from 3,266 cows and 930 sires were considered. Meteorological data from the weather station in closest herd distance were used for the calculation of maximum hourly daily THI, which were allocated to 10 different THI classes. The same FA from 3 stages of lactation were considered as different, but genetically correlated traits. Consequently, a 3-trait reaction norm model was used to estimate genetic parameters and breeding values for FA along THI classes, considering either pedigree ( A ) or genomic ( G ) relationship matrices. De-regressed proofs and genomic estimated breeding values at the intermediate THI class 5 and at the extreme THI class 10 were used as pseudophenotypes in ongoing genomic analyses for thermoneutral (TNC) and heat stress conditions (HSC), respectively. The differences in de-regressed proofs and in genomic estimated breeding values from both THI classes were pseudophenotypes for heat stress response (HSR). Genetic correlations be-tween the same FA under TNC and HSC were smallest in the first lactation stage and ranged from 0.20 for PUFA to 0.87 for SFA when modeling with the A matrix, and from 0.35 for UFA to 0.86 for SFA when modeling with the G matrix. In the first lactation stage, larger additive genetic variances under HSC compared with TNC indicate climate sensitivity for C18:0, PUFA, and UFA. Climate sensitivity was also reflected by pronounced chromosome-wide genetic variances for HSR of PUFA and UFA in the first stage of lactation. For all FA under TNC, HSC, and HSR, quite large genetic variance proportions were explained by BTA14. In GWAS, 30 SNP (within or close to 38 potential candidate genes) overlapped for HSR of the different FA. One unique potential candidate gene ( AMFR ) was detected for HSR of PUFA, 15 for HSR of SFA ( ADGRB1 , DENND3 , DUSP16 , EFR3A , EMP1 , ENSBTAG00000003838 , EPS8 , MGP , PIK3C2G , STYK1 , TMEM71 , GSG1 , SMARCE1 , CCDC57 , and FASN ) and 3 for HSR of UFA ( ENSBTAG00000048091 , PAEP , and EPPK1 ). The identified unique genes play key roles in milk FA synthesis and are associated with disease resistance in dairy cattle. The results suggest consideration of FA in combination with climatic responses when inferring genetic mechanisms of heat stress in dairy cows.


INTRODUCTION
Heat stress as a serious economic issue in dairy cattle farming is associated with a decline in reproductive performances, milk production, and milk quality (Aguilar et al., 2010;Nardone et al., 2010;Negri et al., 2021).Intensive long-term selection in dairy cattle for increased milk yield impairs the maintenance of homeothermy under heat stress, due to the tremendous metabolic heat load of high-yielding cows (Purwanto et al., 1990).Thus, especially high yielding cows are susceptible to the challenging effects of thermal stress (West, 2003;Garner et al., 2016).Heat stress hampers DMI, inducing a state of a negative energy balance (West, 2003).Subsequently, negative energy balance stimulates the mobilization of body fat reserves, implying an increase Genome-wide associations for heat stress response suggest potential candidate genes underlying milk fatty acid composition in dairy cattle M. Bohlouli, 1 K. Halli, 1 T. Yin, 1 N. Gengler, 2 and S. König 1 * of the plasma nonesterified fatty acid concentrations (Bell, 1995;Bielak et al., 2016).In lactating cows, such changes in metabolic pathways result in increasing levels of UFA and, conversely, lower concentrations of SFA in milk (Soyeurt et al., 2008;Gross et al., 2011).Plasma nonesterified fatty acid concentrations in early lactation were significantly associated with milk fat synthesis (Pullen et al., 1989;Adewuyi et al., 2005).Hence, alterations in milk fatty acid (FA) profiles during early lactation are stronger under heat stress conditions (HSC).
The availability of SNP markers on a large scale through the rapidly increasing number of genotyped female cattle enables the detection of potential candidate genes affecting dairy cow performances.Several GWAS focused on the identification of genomic regions and biological mechanisms underlying milk FA traits.For instance, Li et al. (2014) identified 20 novel candidate genes significantly associated with at least one of the short-to long-chain (C4 to C18) SFA and UFA in Chinese Holstein dairy cattle.In the context of a multipopulation GWAS considering a sample of Chinese, Danish, and Dutch Holstein, Gebreyesus et al. (2019) identified potential candidate genes on different chromosomes associated with FA (i.e., DGAT1 on BTA14, ACLY, STAT5A, PRKAA1, and GH on BTA19, and ELOVL3 and ACLS5 on BTA26).
For traits controlled by genes with major effects, the respective chromosome may explain larger proportions of the genetic variance (Visscher et al., 2007;Pimentel et al., 2011).In the context of heat stress, pronounced chromosome-wide genetic variances reflect environmental sensitivity, indicating genomic regions and candidate genes involved in dairy cow trait responses under HSC.Potential candidate genes for heat stress were annotated based on HSR for milk production traits (Macciotta et al., 2017;Sigdel et al., 2019) and rectal temperature (Dikmen et al., 2013;Otto et al., 2019).However, considering the environmental sensitivity of FA in milk for genomic analyses and gene annotations is a novel approach.Consequently, the aim of the present study was to define pseudophenotypes reflecting thermotolerance, heat stress tolerance, and HSR based on breeding values for FA under different climatic conditions.The pseudophenotypes were used as dependent variables in ongoing chromosome-wide genetic parameter estimations, GWAS, and gene annotation approaches.

Cow Traits
The phenotypic data consisted of 39,600 test-day milk FA records (measured in g/100 g of milk yield) for C18:0, PUFA, SFA, and UFA from 5,757 first-lactation Holstein dairy cows.The cows were kept in 16 largescale cooperator herds, located in the region of former East Germany in the federal states Berlin-Brandenburg and Mecklenburg-West Pomerania.The 16 herds reflect the cow genotype and phenotype data basis used for the implementation of cow reference sets in national genomic evaluations in German Holsteins (Yin and König, 2018).The milk samples were collected for 3 yr from 2014 to 2016.Milk FA contents from official test-days were determined by mid-infrared spectra from Foss MilkoScan FT6000 spectrometers at the milk recording center in Güstrow, Germany, using the preinstalled calibration equations.Cows had at least 4 test-day records for the same FA between 5 and 305 DIM.Age at first calving ranged from 20 to 40 mo.The pedigree was traced back for 4 generations, including 27,396 animals (2,291 sires and 19,350 dams).

Genotypes
Genotyping was conducted using either the Illumina Bovine 50K SNP BeadChip V2 (419 animals) or the Illumina Bovine Eurogenomics 10K low-density chip (3,777 animals).Animals with low-density genotypes were imputed to 50K (incorporated into the routine procedure for official national genetic evaluations in Germany and considering 50K genotypes from more than 44,000 bulls) as implemented by project partner VIT Verden (Segelke et al., 2012).After imputation, 45,613 SNP were available from 4,196 genotyped animals.The quality controls of SNP were performed using the software package PLINK (Purcell et al., 2007).The SNP located on the X and Y chromosomes, SNP with minor allele frequency lower than 0.05, and SNP significantly deviating from the Hardy-Weinberg equilibrium (P < 0.001), were discarded.All genotyped animals and SNP had call rates larger than 95%.Finally, 40,523 SNP from 3,266 cows (2,977 cows with phenotypic records) and 930 sires were available for the genomic analyses.

THI
Longitude and latitude for weather stations and farms were fed into the Geosphere package in R (Hijmans et al., 2016) to identify the nearest weather station for each farm.In this regard, 13 different weather stations were allocated to the 16 different farms.The minimum distance between a farm and the nearest weather station was 6.5 km and the maximum distance was 26.8 km (average distance: 16.6).The hourly THI was calculated using meteorological data as follows (NRC, 1971): where T is the temperature and RH is the relative humidity recorded in hourly intervals.The maximum hourly THI from 4 d (the test-day and the previous 3 d) were averaged (as suggested by Bohmanova et al., 2008) and merged with the respective test-day for the milk FA traits.
The averaged maximum THI ranged from 31.62 to 75.4.Small proportions of phenotypic records (<1%) were available at both extreme ends of the THI scale.Therefore, to avoid possible artifacts of reaction norm models on (co)variance component estimates at extreme THI, we created 10 distinct classes for the averaged maximum THI.Each THI class consisted of about 3,960 test-day records per FA trait.Cows had phenotypic records in at least 3 different THI classes.Of all cows, 35% and 29% had repeated records in 5 and 6 different THI classes, respectively.

Estimation of Genetic Parameters
The same traits from different stages of lactation (3 stages were defined: 5-105 DIM, 106-205 DIM, and 206-305 DIM) were considered as different, but genetically correlated traits.Consequently, a 3-trait reaction norm model for the same FA in different lactation stages was fitted to estimate genetic parameters and breeding values along THI.In matrix notation, model 1 was defined as follows: where y i is the vector of observations at the ith stage of lactation (i = 1, 2, or 3) considered separately for each FA trait (C18:0, PUFA, SFA, or UFA); b i is the vector of fixed effects including herd-test-day, milking frequency (2 or 3 times per day), 10-d intervals for DIM within lactation (30 classes), calving season and fixed regressions on age at first calving; a i and p i are the vectors of random regression coefficients for additive genetic and permanent environmental effects on THI classes, respectively; X i is an incidence matrix for fixed effects; Z i and W i is the covariable matrices for additive genetic and permanent environmental, respectively, and e i is the random residual effect.The rows of the Z i and W i matrices contain the quadratic Legendre polynomials.Random effects were assumed to be normally distributed with zero means.The covariance structure of random effects was defined as follows: , where a, p, and e are additive genetic, permanent environmental, and residual effects, respectively; W a and W p are 9 × 9 covariance matrices for the random regression coefficients of additive genetic and permanent environmental effects, respectively; R is a 3 × 3 covariance matrix of random residual effects; K is the pedigree (A) or genomic (G) relationship matrix among individuals; I is an identity matrix for the permanent environmental and residual effects, and ⊗ denotes the Kronecker product.The G matrix was constructed as proposed by Yang et al. (2010), that is, using the "option which G3." Genetic parameters were estimated using the AIREMLF90 program (Misztal et al., 2002), considering a convergence criterion of 10 −12 .The negative inverse of the average information matrix provides estimates for the variances of estimated parameter.Approximate standard errors were computed based on the estimated parameters and their variances according to Fischer et al. (2004).

Accuracy of Breeding Values and Calculation of DRP
For each FA at lactation stage 1, 2, and 3, accuracies of EBV r i t ( ) from the A matrix and of genomic esti- mated breeding values from the G matrix (GEBV) for the ith animal in the tth THI class were calculated as , where PEV i t was the prediction error variance for the ith animal in the tth THI class, and σ a t side of the mixed model equations.According to Tier and Meyer (2004), PEV i t was the diagonal of QT i Q′, where Q was a 10 × 3 matrix with the values of 3 coefficients of the second-order Legendre polynomials for THI classes 1 to 10, and T i was a 3 × 3 matrix of random regression coefficients of prediction error (co)variances for the ith animal.The de-regressed proof (DRP) for the ith animal in the tth THI class DRP i t ( ) was calculated for the pedi- gree-based EBV according to Garrick et al. (2009) as , where EBV i t was the EBV for the ith animal in the tth THI class.

Pseudophenotypes
In a previous study, we identified remarkable genotype by environment interactions for C18:0, PUFA, and UFA between THI class 10 and THI class 5 (Bohlouli et al., 2021).Consequently, in the present study, the THI class 5 (THI range from 50.5 to 52.5) was chosen as an environmental descriptor for thermoneutral conditions (TNC) and the THI class 10 (THI range from 68.4 to 75.4) for HSC.The descriptive statistics for the defined climate conditions are given in Table 1.The DRP and GEBV of the genotyped cows and sires for the same FA in THI class 5 and THI class 10 were used as pseu-dophenotypes in ongoing genomic analyses for TNC and HSC, respectively.In addition, the differences in GEBV and in DRP from both THI classes of these animals were considered as pseudophenotypes for HSR.To perform reliable GWAS, we only considered EBV and GEBV with accuracies larger than 0.60 for the creation of pseudophenotypes.Accordingly, the number of records for DRP was smaller than for GEBV (Table 2).

Genome-Wide Associations
The GWAS was performed using the linear mixed model association method as implemented in the GCTA software (Yang et al., 2011a).The statistical model 2 was where y is the vector of pseudophenotypes (DRP or GEBV) for C18:0, PUFA, SFA, or UFA under specific environmental conditions (TNC, HSC, or HSR) at the given lactation stage (stage 1, 2, or 3); 1 is the vector of ones; µ is the overall mean effect; x i is the vector of genotypes coded as 0, 1, or 2; s i is the effect of the ith SNP to be tested for association; ) is the vector of random polygenic effects, in which K is the genetic relationship matrix and σ u 2 is the polygenic vari- ance estimated from the null model (y = 1µ + Zu + e); Z is an incidence matrix relating phenotypes to the corresponding random polygenic effects, and e is the residual random effect.Matrix K was constructed based on the A matrix for DRP and on the G matrix for GEBV.
Manhattan plots for the −log 10 P-values of the tested SNP were created using the ggplot2 package in R (Wickham, 2016).The genomic inflation factor was defined as the median of the observed chi-squared test statistics divided by the expected median chi-squared distribution (0.4549).Chi-squared test statistics were calculated from P-values assuming one degree of freedom.A Bonferroni correction was applied to account for multiple testing in the GWAS.Accordingly, the Bonferroni-corrected significance threshold was defined as 0.05/N, where N is the total number of SNP.

Genetic Parameters for Individual Chromosomes
Genomic restricted maximum likelihood as implemented in GCTA was used to estimate variance components for each of the 29 autosomes.All chromosomes were considered simultaneously in joint analyses via model 3 as follows: where y is a vector of pseudophenotypes (DRP or GEBV) for C18:0, PUFA, SFA, or UFA under specific environmental conditions (TNC, HSC, or HSR) at the given lactation stage (stage 1, 2, or 3); µ is the overall mean effect; u i is a vector of additive genetic effects attributed to the ith chromosome with variance of , where G i is the genomic relationship matrix constructed from the SNP located on the ith chromosome, and σ u i 2 is the additive genetic variance explained by the SNP located on the ith chromosome.

Annotation of Potential Candidate Genes
Significant SNP (Bonferroni-corrected level) from model 2 were used to annotate potential candidate genes.For the annotation of potential candidate genes, windows of 400 kb consisting of 200 kb upstream and 200 kb downstream of the significant SNP were defined.The SNP were mapped to corresponding genes from the Bos taurus ARS-UCD1.2annotation release 96 assembly from the Ensembl database (http: / / www .ensembl.org/biomart/ martview), using the R package biomaRt (Durinck et al., 2009).The VennDiagram package in R (Chen, 2018) was used to illustrate overlapping and unique significant SNP and potential candidate genes across traits, across stages of lactation, and across environmental conditions (TNC, HSC, and HSR).

FA Responses by THI and Lactation Stage
Phenotypically, the PUFA content was constant across lactation stages.Contents for C18:0 and UFA were larger in the first than in the third stage of lactation (Table 1).Conversely, the highest SFA content was observed in later lactation stages.The C18:0 content was larger under TNC than under HSC (0.35 ± 0.10 vs. 0.31 ± 0.09 g/100 g of milk).A decline under HSC was also observed for SFA, especially in the early lactation stage.Accordingly, Soyeurt et al. (2008) reported the lowest SFA content during summer, whereas the content of UFA was larger during spring and summer compared with the winter season.Heat stress conditions cause a negative energy balance, with ongoing effects on milk yield and FA contents (Penasa et al., 2015).In the stage of a negative energy balance, the mobilization of adipose FA increases, and contributes to an increase of UFA and C18 contents in milk (Bastin et al., 2013).Thus, especially in the early lactation stage, heat stress intensifies alterations in milk FA contents.

Genetic Parameters for FA
The estimated variance components and heritabilities based on the A and G matrices are presented in Table 3.For the FA in the first stage of lactation under TNC, the heritabilities ranged from 0.11 for PUFA to 0.21 for SFA with the A matrix, and from 0.11 for PUFA to 0.30 for SFA with the G matrix.Among all traits in the first lactation stage under HSC, the largest heritability of 0.32 was estimated for SFA when considering genomic relationships.Generally, heritabilities for all FA were larger in the later than in the early lactation stage.In the early lactation stage, genetic variances were larger for C18:0, PUFA, and UFA under HSC compared with the respective estimates under TNC.Oppositely, in the late lactation stage, genetic variances for all FA were larger under TNC than under HSC.Independent from lactation stage effects, Hammami et al. (2015) and Bohlouli et al. (2021) estimated largest genetic variances for C18:0 and UFA at high THI, whereas estimates were smallest in the thermoneutral zone.Generally, heritabilities for FA are in line with previously reported estimates in other cattle breeds located in middle Europe (Soyeurt et al., 2008;Bastin et al., 2013).Bastin et al. (2013) reported quite large heritabilities of 0.43 for C16:0 and 0.46 for SFA, and moderate heritabilities of 0.26 for MUFA, 0.31 for PUFA, and 0.27 for UFA.Significant genetic effects on variations for SFA were identified by Palmquist et al. (1993) and Dewhurst et al. (2007), explaining the moderate to large heritabilities.Conversely, UFA contents were more strongly affected by diet compositions and environmental components (Krag et al., 2013;Hammami et al., 2015), explaining the smaller heritabilities for PUFA and UFA.With regard to the largest genetic variances for C18:0, PUFA, and UFA in the first lactation stage under HSC, selection on these traits in harsh heat stress environments might contribute to improvements of thermotolerance.
Using the A matrix, genetic correlations between the same traits in the first stage of lactation under TNC and HSC ranged from 0.20 for PUFA to 0.87 for SFA (Table 3).For all FA, the largest genetic correlations (>0.84) between the same traits under TNC and HSC were estimated in the second stage of lactation.In the first stage of lactation, genetic correlations between same traits under TNC and HSR were larger with the G matrix than with the A matrix.For the remaining stages, genetic correlations from the G matrix were similar to corresponding values from the A matrix.These results reflect estimates by Bohlouli et al. (2021), who used a combined pedigree-genomic relationship matrix.Contrary to the A matrix, unrelated individuals are also connected genetically through the G matrix.In addition, only a subset of phenotyped animals from pedigree-based analyses were considered when constructing the G matrix.Such differences between G and A matrices might explain the minor differences in genetic parameters estimates.Furthermore, approximated standard errors of genetic parameters were smaller from the G matrix, as previously reported (Forni et al., 2011;Aldridge et al., 2020).With regard to the first stage of lactation, the large additive genetic variances under HSC and the low genetic correlations between yields under TNC and HSC indicate climate sensitivity for C18:0, PUFA, and UFA.Hence, genetic merits for animals for these traits in the first lactation stage can differ under TNC and HSC.As indicated above, the bioenergetic status during heat stress is comparable with metabolic stress, which cows experience during early lactation (Moore et al., 2005;Hammami et al., 2015).

Chromosome-Wide Genetic Variances
Using both GEBV (Figure 1) and DRP (Figure 2), larger chromosome-wide genetic variances were estimated for FA under TNC and HSC when compared with the estimates for the respective traits for HSR.
In the first stage of lactation, chromosome-wide genetic variances were larger for PUFA and UFA under HSC than under TNC.For SFA under TNC and HSC, genetic variances were smallest in the early lactation stage.In the later lactation stages, chromosome-wide genetic variances were very similar for the same FA under TNC and HSC.Substantial chromosome-wide genetic variances were estimated for HSR of PUFA and UFA in the first lactation stage (Figure 1), reflecting climate sensitivity.For HSR of C18:0 and SFA, the largest genetic variances were estimated in the third stage of lactation.
The relationship between the amount of genetic variance attributed to a chromosome and the chromosome length was weak, supporting results for milk production traits (Jensen et al., 2012).Larger genetic variances for some chromosomes than expected (according to their chromosomal length) indicate that genes with major effects on FA profiles are unequally distributed across the genome.The very large genetic variance explained by BTA14 for all FA was the main reason for the weak relationship between chromosomal variance and chromosome length.It is well known that BTA14 carries several major genes with large effects on milk yield and fat percentage in dairy cattle (Grisart et al., 2002;Pimentel et al., 2011).A major part of the genetic variance for BTA14 is due to the DGAT1 gene, which has a major effect on milk production traits (Grisart et al., 2002) and milk FA profiles (Gebreyesus et al., 2019).Grisart et al. (2002) reported that a nonconservative K232A polymorphism in DGAT1 explains 18%, 15%, and 8% of the variances in daughter yield deviations for milk yield, fat yield, and protein yield, respectively.
For all FA in the first stage of lactation under different environmental conditions (TNC, HSC, and HSR), more than 5% of the total genetic variance was assigned to BTA3.The amount of genetic variance explained by BTA5 (more than 5% of the total genetic variance) was notable for the FA in different lactation stages under TNC and HSC.However, for HSR of FA in the early lactation stage, BTA5 explained less than 5% of the total genetic variation.In a study by Pimentel et al. (2011), BTA5 contributed to a large proportion of genetic variation for fat percentage and SCC in a pre-selected group of Holstein dairy cattle.For UFA under different environmental conditions, large proportions of genetic variation were due to the effects on BTA11.For instance, BTA11 explained 5.83%, 4.80%, and 4.93% of the total genetic variances for HSR of UFA in the first, second, and third stages of lactation, respectively.As explained later, 2 potential candidate genes (ENSBTAG00000048091 and PAEP) on BTA11 are associated with HSR of UFA.With regard to HSR, the contribution of each chromosome to the trait heri- ( ) variances with corresponding approximate SE (in parentheses) for C18:0, PUFA, SFA, UFA, and genetic correlations (r g ) for the same traits between thermoneutral (TNC) and heat stress (HSC) conditions

Genome-Wide Associations
The −log 10 P-values (Manhattan plots) for SNP associations for C18:0, PUFA, SFA, and UFA using GEBV are presented in Supplemental Figure S1 (https: / / jlupub .ub.uni-giessen .de// handle/ jlupub/ 559), and in Figures 3, 4, 5, and 6 for DRP.For all FA, the plots display the most pronounced peak on BTA14.In comparison to the association analyses based on GEBV, the number of significant SNP was larger when using DRP as pseudophenotypes.The method for calculating DRP (i.e., using the EBV reliability in the denominator) implies larger variations for DRP, which may contribute to a larger number of significant SNP.Apart from a few exceptions, all significant SNP (according to the Bonferroni correction) for GEBV were also sig-nificant for DRP.One exception addresses HSR of UFA in the first stage of lactation, where 5 significant SNP were detected on BTA14 for GEBV, but not for DRP.For both pseudophenotypes DRP and GEBV of PUFA and UFA in the first lactation stage, a larger number of significant SNP was detected under HSC than under TNC.For DRP and GEBV in the second and the third lactation stage, an almost identical number of significant SNP was detected for the same FA under TNC and HSC.With regard to the whole lactation period and HSR, the number of significant SNP for DRP varied from 31 (SNP located within or in close distance to 39 potential candidate genes) for PUFA to 121 (SNP located within or close to 65 potential candidate genes) for SFA.For GEBV, the number of significant SNP varied from 36 (SNP located within or close to 40 potential candidate genes) for PUFA to 51 (SNP located within or close to 48 potential candidate genes) for SFA.For both pseudophenotypes DRP and GEBV, 52% of significant SNP for all FA are located within genes.
Inflation factors for GWAS based on DRP in different lactation stages ranged from 1.08 to 1.36 for C18:0, from 1.03 to 1.24 for PUFA, from 1.02 to 1.32 for SFA, and from 1.05 to 1.28 for UFA, indicating only small to moderate inflated P-values.For GEBV, inflation factors ranged from 0.79 to 1.02.In a GWAS for cattle weight traits (Yin and König, 2019), inflation factors were quite large (range from 1.33 to 1.96) when using DRP based on EBV from the A matrix.In the same study, inflation factors were slightly lower than 1.0 when using DRP derived from GEBV, supporting the results from the present GWAS for FA and results by Nayeri et al. (2016) for milk production traits.In GWAS, inflation factors depend on the population structure, the sample size, the polygenic nature of the trait, the heritability and the linkage disequilibrium pattern (Yang et al., 2011b).For both pseudophenotypes of the same traits in different lactation stages, inflation factors for HSR were close to 1.0, probably due to the smaller additive genetic variances for HSR compared with TNC and HSC (Figure 1).Furthermore, pseudophenotypes as used in the GWAS under TNC and HSC had quite large accuracies, which might explain the deviation of the test statistics from its expected values (Jardim et al., 2018).
Generally, the Manhattan plots were quite similar for both pseudophenotypes DRP and GEBV.Minor differences address the GBLUP methodology, where connections through the G matrix among unrelated animals may cause fluctuations in genetic evaluations (Misztal et al., 2020).Furthermore, when computing DRP, heterogeneous variances are taken into account due to differences in EBV reliabilities (Garrick et al., 2009).

Potential Candidate Genes for Heat Stress Response Across FA Traits
For HSR in the first lactation stage, all significant SNP for C18:0 were also significantly associated with PUFA (Figure 7).In the early lactation stage, 27 potential candidate genes overlapped between PUFA and C18:0.No significant SNP was detected for HSR of SFA and UFA in the first lactation stage.For the second lactation stage, 52 significant SNP (located within or close to 45 potential candidate genes) overlapped for C18:0, SFA, and UFA.For the third lactation stage, all of the 25 potential genes for C18:0 were detected as potential candidate genes for SFA and UFA.For HSR in all lactation stages, 30 SNP (located within or close to 38 potential candidate genes) overlapped between the different FA.Results were very similar when using GEBV as pseudophenotypes (Supplemental Figure S2, https: / / jlupub .ub.uni-giessen .de// handle/ jlupub/ 559).
Most of the potential candidate genes detected for the FA support results from previous GWAS for milk production (Nayeri et al., 2016;Atashi et al., 2020;Shabalina, 2021) and FA (Cruz et al., 2019;Freitas et al., 2020) traits in Holstein dairy cattle.In the present study, 3 SNP (rs109421300, rs109146371, and rs109350371) on BTA14 showed the highest associations with all FA.The most significant SNP (rs109421300) is located within an intron of DGAT1.Two potential candidate genes (PPP1R16A and FOXH1) are located in close distance to the second most significant SNP (rs109146371).The third most significant SNP (rs109350371) is located in the PLEC gene.The mentioned 3 significant SNP display the largest effects on the FA in all lactation stages under different environmental conditions, explaining the large genetic variances attributed to BTA14.
The genes DGAT1, PPP1R16A, FOXH1, and PLEC play key roles in the lactation process (Cases et al., 1998;Jiang et al., 2019;Shabalina, 2021).Accordingly, associations of these genes with FA were reported in previous studies (Bovenhuis et al., 2016;Palombo et al., 2018;Sanchez et al., 2018;Cruz et al., 2019).DGAT1 is a major gene influencing inflammatory response and lipid metabolism and encodes an enzyme that catalyzes the final step of triacylglycerol synthesis (Cases et al., 1998), and significantly affected energy metabolism in early lactation in high yielding cows (Lešková et al., 2013).
Some genes were detected as potential candidates for HSR of specific FA.The autocrine motility factor receptor (AMFR) gene on BTA18 was associated with PUFA only.The expression of the AMFR gene was found to be significantly associated with fat percentage and protein percentage (Cole et al., 2011).Twelve potential candidate genes (ADGRB1, DENND3, DUSP16, EFR3A, EMP1, ENSBTAG00000003838, EPS8, MGP, PIK3C2G, STYK1, TMEM71, and GSG1) on BTA5 and 3 potential candidate genes (SMARCE1, CCDC57, and FASN) on BTA19 only were associated with HSR of SFA.These genes were reported as potential candidate genes for milk production traits (Nayeri et al., 2019;Oliveira et al., 2019;Shabalina, 2021).Additionally, PIK3C2G was associated with C10:0 contents (Li et al., 2014).Our results also agree with the significant associations for FA as outlined by Bouwman et al. (2014), Knutsen et al. (2018), andBhuiyan et al. (2018), who identified CCDC57 and FASN as candidate genes for C14:0.A significant SNP (rs41921177) is positioned in the CCDC57 gene, and located in close distance to the FASN gene.The CCDC57 gene is expressed in mammary glands of lactating cows (Medrano et al., 2010) and FASN gene plays a key role in fat synthesis.The FASN gene encodes a multifunctional enzyme that catalyzes the de novo biosynthesis of long-chain SFA (Semenkovich, 1997).The SNP rs110436636, located in close distance to 2 potential candidate genes (ENS-BTAG00000048091 and PAEP) on BTA11, was only significantly associated with HSR of UFA.The effect of PAEP on milk production traits was reported by Korkuć et al. (2021) and by Shabalina (2021).Knutsen et al. (2018) andSanchez et al. (2018) reported sig-nificant associations between PAEP and FA traits.The PAEP gene encodes the milk protein β-LG, which is the major whey protein in milk (Knutsen et al., 2018) and an intracellular transporter of FA (Le Maux et al., 2014;Knutsen et al., 2018).ENSBTAG00000048091 is a novel gene and is a potential duplicate of PAEP (Xiang et al., 2020).

Potential Candidate Genes for Heat Stress Response Across Lactation Stages
Most of the potential candidate genes were detected consistently in all lactation stages, under TNC as well as under HSC for both pseudophenotypes DRP and GEBV.Hence, the same genes are affecting physiological processes of milk FA synthesis along the entire lactation period under both climatic conditions.However, with regard to HSR of FA at different lactation stages, we detected different sets of potential candidate genes (Figure 8).
For HSR of C18:0 in the second lactation stage, 53 significant SNP (located within or close to 45 potential candidate genes) were detected.A subset of 20 and 18 SNP was also significant in the first and the third lactation stages, respectively.Thirty SNP on BTA14 and 1 SNP (rs109072828) located in the gene AMFR on BTA18 were significantly associated with HSR of PUFA in the first and the third lactation stages, respectively (Figure 8).From the 121 significant SNP for HSR of SFA, 77 SNP overlapped between the second and the third lactation stages.No significant SNP was detected for HSR of SFA in the early lactation period.For HSR of UFA, 71 significant SNP (located within or close to 53 potential candidate genes) were detected in the third lactation stage.subset of 64 SNP (located within or close to 49 potential candidate genes) was significant in the second lactation stage.Results were very similar when basing the association analyses on GEBV (Supplemental Figure S2).
The sign of significant SNP effects (β) was the same for each trait under both TNC and HSC, and no changes in the sign were observed across lactation stages (not shown).For HSR of FA, the signs of SNP effect for estimates in the first lactation stage partly differed from the signs from the remaining stages.For example, the effect of the SNP rs109421300 located in the DGAT1 gene on C18:0 in the first lactation stage (β = 4.0e-03) was opposite to those estimated for the second (β = −9.6e-04)and third (β = −5.6e-04)stages.With regard to HSR, the significantly associated SNP located within or in close distance to the genes PPP1R16A, FOXH1, CYHR1, ZNF623, ZNF696, GLI4, KCNK9, TSNARE1, and TRAPPC9 display differing signs for the estimates across lactation.The results indicate that the genes related to FA are differentially expressed in response to heat stress in different stages of lactation as reflected by the low genetic correlations between the same traits from the early and later lactation periods (Bohlouli et al., 2021).

Potential Candidate Genes Across TNC, HSC, and HSR
The number of significant SNP and genes associated with the same FA in the entire lactation overlapping across TNC, HSC, and HSR are shown in Figure 9 (for DRP) and in Supplemental Figure S2 (for GEBV).For the entire lactation, all 53 significant SNP (located within or close to 45 potential candidate genes) associ- ated with HSR of C18:0 were also significant under TNC and HSC (Figure 9).For PUFA, the potential candidate gene NPFFR2 (neuropeptide FF receptor 2) and the group specific component genes on BTA6 were detected under TNC as well as under HSC.These 2 genes are located within a QTL, which is associated with milk production and clinical mastitis in dairy cattle (Olsen et al., 2016).
The KCNK9 (potassium channel, subfamily K, member 9) gene was inferred as a potential candidate gene for all FA across lactation under TNC as well as under HSC.Regarding HSR, KCNK9 was associated with PUFA only in the first lactation stage, and with C18:0, SFA, and UFA in the later lactation stages.The KCNK9 gene is a candidate gene for milk yield and fat percentage (Jiang et al., 2014) and milk FA (Freitas et al., 2020) in Holstein dairy cattle.The KCNK9 gene encodes a protein that modulates aldosterone secretion.Aldosterone plays a key role in regulating sodium homeostasis in the distal nephron.Ravanelli et al. (2021) showed that heat stress causes an increase in aldosterone secretion.The genes TSNARE1 (T-SNARE domain-containing protein 1) and TRAPPC9 (traffick-ing protein particle complex 9) were associated with all FA in different lactation stages under both TNC and HSC.TSNARE1 and TRAPPC9 are described as candidate genes for milk production (Jiang et al., 2014(Jiang et al., , 2019;;Shabalina, 2021) and FA (Freitas et al., 2020) in Holstein dairy cattle.TSNARE1 plays a role in intracellular protein transport and synaptic vesicle exocytosis (Luo et al., 2021).Luo et al. (2021)  .Previous studies demonstrated that heat stress reduces the antioxidant body capacity, with detrimental effects on immune functions (Zhang et al., 2014;Dahl et al., 2020).Most of the identified potential candidate genes for FA across climatic conditions (TRAPPC9, TSNARE1, LY6D, LYNX1, LYPD2, ARC, GC, and NPFFR2) were related with mastitis susceptibility in dairy cattle (Tiezzi et al., 2015;Wang et al., 2015;Olsen et al., 2016).Consequently, our finding supports reports on similar genetic mechanisms underlying HSR and mastitis resistances (Dahl et al., 2020;Rakib et al., 2020).

CONCLUSIONS
In the first stage of lactation, larger additive genetic variances were estimated for C18:0, PUFA, and UFA under HSC compared with TNC, reflecting the environmental sensitivity during early lactation in high yielding cows.On the chromosome-wide level, larger genetic variances were estimated for HSR of PUFA and UFA in early lactation compared with the later lactation stages.The significant SNP markers on BTA14, which are located within or close to the potential candidate genes DGAT1, PPP1R16A, FOXH1, and PLEC, had the largest effects on the FA across lactation under TNC, HSC, and HSR.With regard to GWAS for HSR, we identified 38 potential candidate genes associated with all FA traits.The results from the present study also  sequence data for the exploration of physiological functions and on expressions of the identified potential candidate genes.
their approximate SE are multiplied by 10 −3 .tability varied across lactation stages, indicating different genetic mechanisms for FA trait responses under heat stress in the course of lactation.The partitioning of genetic variances onto individual chromosomes contributes to the detection of the chromosomes harboring potential candidate genes.

Figure 3 .
Figure 3. Manhattan plot of −log 10 P-values of the tested SNP markers for C18:0 FA in different lactation stages under thermoneutral conditions (TNC), under heat stress conditions (HSC), and for heat stress response (HSR).The red dashed line is the significance threshold for the Bonferroni correction.De-regressed proofs were used as pseudophenotypes.

Figure 4 .
Figure 4. Manhattan plot of −log 10 P-values of the tested SNP markers for PUFA in different lactation stages under thermoneutral conditions (TNC), under heat stress conditions (HSC), and for heat stress response (HSR).The red dashed line is the significance threshold for the Bonferroni correction.De-regressed proofs were used as pseudophenotypes.

Figure 5 .
Figure 5. Manhattan plot of −log 10 P-values of the tested SNP markers for SFA in different lactation stages under thermoneutral conditions (TNC), under heat stress conditions (HSC), and for heat stress response (HSR).The red dashed line is the significance threshold for the Bonferroni correction.De-regressed proofs were used as pseudophenotypes.
evaluated physiological indicators of heat stress in Holstein dairy cattle and reported a strong association of TSNARE1 with rectal temperature.The TRAPPC9 gene plays a key role in the regulation of cell differentiation and might be relevant for the control of immune system mechanisms [S.Mi, M. Song, Y. Dong, Z. Zhang, L. Fan, L. Shi (China Agricultural University), X. Wang (Northwest Agriculture and Forestry University), K. Shi, and Y. Yu (China Agricultural University), unpublished data]

Figure 6 .
Figure 6.Manhattan plot of −log 10 P-values of the tested SNP markers for UFA in different lactation stages under thermoneutral conditions (TNC), under heat stress conditions (HSC), and for heat stress response (HSR).The red dashed line is the significance threshold for the Bonferroni correction.De-regressed proofs were used as pseudophenotypes.

Figure 7 .
Figure 7. Venn diagrams for the number of SNP markers and genes associated with heat stress response overlapping across C18:0, PUFA, SFA, and UFA.De-regressed proofs were used as pseudophenotypes.

Figure 8 .
Figure 8. Venn diagrams for the number of SNP markers and genes associated with heat stress response for C18:0, PUFA, SFA, and UFA overlapping across the first (S1), second (S2), and third (S3) stages of lactation.De-regressed proofs were used as pseudophenotypes.
Bohlouli et al.: HEAT STRESS RESPONSE OF FATTY ACIDS

Table 1 .
Bohlouli et al.: HEAT STRESS RESPONSE OF FATTY ACIDS Number of records (N) and means (SD in parentheses) for C18:0, PUFA, SFA, and UFA at different stages of lactation under different climate conditions

Table 2 .
Number of genotyped animals, sires, cows, and SNP markers used in genome-wide association analyses based on de-regressed proofs (DRP) and genomic estimated breeding values (GEBV) for C18:0, PUFA, SFA, and UFA

Table 3 .
Bohlouli et al.: HEAT STRESS RESPONSE OF FATTY ACIDS Bohlouli et al.: HEAT STRESS RESPONSE OF FATTY ACIDS Heritabilities