Identification of candidate genes on the basis of SNP by time-lagged heat stress interactions for milk production traits in German Holstein cattle

The aim of this study was to estimate genotype by time-lagged heat stress (HS) variance components as well as main and interaction SNP-marker effects for maternal HS during the last eight weeks of cow pregnancy, considering milk production traits recorded in the offspring generation. The HS indicator was the temperature humidity index (THI) for each week. A dummy variable with the code = 1 for the respective week for THI ≥ 60 indicated HS, otherwise, for no HS, the code = 0 was assigned. The dataset included test-day and lactation production traits from 14,188 genotyped first parity Holstein cows. After genotype quality control, 41,139 SNP markers remained for the genomic analyses. Genomic animal models without (model VC_nHS) and with in-utero HS effects (model VC_wHS) were applied to estimate variance components. Accordingly, for genome-wide associations, models GWA_nHS and GWA_wHS, respectively, were applied to estimate main and interaction SNP effects. Common genomic and residual variances for the same traits were very similar from models VC_nHS and VC_wHS. Genotype by HS interaction variances varied, depending on the week with in-utero HS. Among all traits, lactation milk yield with HS from week 5 displayed the largest proportion for interaction variances (0.07). For main effects from model GWA_wHS, 380 SNPs were suggestively associated with all production traits. For the SNP interaction effects from model GWA_wHS, we identified 31 suggestive SNPs, which were located in close distance to 62 potential candidate genes. The inferred candidate genes have various biological functions, including mechanisms of immune response, growth processes and disease resistance. Two biological processes excessively represented in the overrepresentation tests addressed lymphocyte and monocyte chemotaxis, ultimately affecting immune response. The modelling approach considering time-lagged genotype by HS interactions for production traits inferred physiological mechanisms being associated with health and immunity, enabling improvements in selection of robust animals.


Introduction
The tripled number of days with extreme temperatures in European countries from 1950 to 2018 [1] indicates the increasing importance of studying the effect of climate change on primary and functional traits in dairy cattle. Comparing to other domesticated livestock species, high productive Holstein Friesian cattle responded very sensitive to heat stress (HS), especially during the challenging early lactation period with temperatures beyond the thermoneutral zone lower than 5˚C or higher than 25˚C [2]. Selection on improved milk production simultaneously increased feed intake and body weight, implying a raised metabolic heat increment and a decline of the thermoneutral temperature range [3]. Furthermore, in addition to temperature, dairy cow trait alterations were influenced by humidity, suggesting consideration of a temperature-humidity index (THI) in genetic HS studies [4].
Response to direct HS, i.e., HS close to the trait recording date, was observed when THI exceeds a certain threshold. For example, in Holstein cows, protein yield decreased at THI 68 [5]. Genetically, random regression models were applied to infer genetic (co)variance components along a continuous THI scale [6]. Such modelling approach was also used to detect possible genotype x environment interactions (GxE), and THI dependent re-rankings of sires have been observed [7]. Genomically, genetic markers significantly associated with HS response and underlying candidate genes were detected via genome-wide association studies (GWAS). Four single nucleotide polymorphisms (SNPs) contributing to a milk yield decline under HS were identified by Hayes et al. [8] in Australian Holstein Friesian cattle. The four SNPs were located on Bos taurus autosomes (BTA) 8, 10, 25, and 29. Among them, the SNP BFGL-NGS-30169 on BTA29 was also identified in Jersey cattle, and the strongest annotated candidate gene for the variation of milk yield under HS was the fibroblast growth factor 4 [8]. Selection of cattle displaying the ability to regulate body temperature under HS will contribute to improved heat resistance genetically, implying only marginal detrimental effects on primary and functional traits. Thus, the genetic architecture of traits reflecting the thermoregulation ability under HS, such as rectal temperature and respiration rate, was investigated in Holstein and Gir x Holstein crossbred cattle [9,10]. Rectal temperature was used an indicator for HS [11], and respiration rate reflects the ability to maintain body temperature through evaporative cooling [10].
Most of the genetic studies focused on immediate responses to HS, i.e., considering THI at the measuring day, or the average THI shortly before the measuring day. Phenotypically, HS during late gestation in dams, termed as maternal or in-utero HS, significantly impaired growth, metabolism, immune functions and survival in offspring [12][13][14], since this period plays a crucial role in fetal growth [15]. Monteiro et al. [13] studied long-term effects of inutero HS on female fertility and performance traits recorded in offspring and identified most detrimental impact of HS from the last 6 weeks of gestation. Quantitative-genetically, Halli et al. [16] identified alterations of genetic co(variance) components for weight and growth traits of dual-purpose cattle due to HS impact on their dams during the late pregnancy period.
In this regard, the current study focusses on enhanced genomic modelling approaches to inferring time-lagged HS impact, and considers a large dataset of genotyped Holstein Friesian cows kept in large-scale contract herds. Specifically, we firstly aim on the detection of significant SNP for the direct (main) and the interaction effect in the context of time-lagged HS on milk production traits at the first official test-day after calving and on first lactation production records. The main effect represents the SNP effect for the respective trait expressed consistently, independent from time-lagged HS. The interaction component is the difference of the SNP effects between cows undergoing in-utero HS or not. Afterwards, on the basis of the identified SNP associations, we annotated potential candidate genes. Against this background, i.e.,

Genotype data
The cows were genotyped using the Illumina BovineSNP50 v2 BeadChip (3,775 animals), or the Illumina Bovine Eurogenomics 10K low-density chip (10,413 animals). The 10K SNP genotypes were imputed to the 50K SNP panel by project partner vit (Verden, Germany), as implemented in the process for national routine genetic evaluations [17]. SNP quality controls were performed using the preGSf90 program from the BLUPf90 package [18,19]. Filtering criteria for markers were as follows: consideration only of SNPs located on Bos taurus autosomes, minor allele frequency larger than 0.01, minimum animal and SNP call rate of 0.95 and no significant deviation from Hardy-Weinberg equilibrium (the difference between observed and expected heterozygous frequencies was smaller than 0.15). Furthermore, we deleted cows pairs

Meteorological data
Pairwise distances (in km) were calculated between weather stations and cow herds. The calculation based on coordinates for the respective longitude and latitude of each herd and weather station, and was performed using the GEOSPHERE package in R [21]. According to the minimal distances, we allocated 32 weather stations to 53 herds. The maximum distance between a herd and a weather station was 27.88 km, the minimum distance was 0.74 km and the average distance was 14.79 km. Hourly THI were calculated considering hourly temperature (T) and relative humidity (RH) as follows [22]: The daily THI was computed by averaging hourly THI over 24 hours. Afterwards, we calculated the average daily THI within the following weeks during late gestation of the respective dam, or, in other words, before the birth date of the genotyped cows: 0-7 days (WK1),  Table 2.

Statistical models
Variance components and genomic heritabilities. A genomic animal model was applied to estimate variance components for milk production and TSCS. In this regard, we used the AIREML algorithm as implemented in AIREMLf90 from the BLUPf90 package [19]. The statistical model in matrix notation was: where y = a vector of observations for the test-day or the lactation records; b = a vector of fixed effects including herd, calving year, calving month, age at first calving classes, lactation-calving-age classes for dams, classes for days in milk (for the test-day traits), lactation length classes

PLOS ONE
SNP by time-lagged heat stress interactions for milk production traits (for lactation traits) and a dummy variable indicating the presence of HS during 1 to 8 weeks before birth for the genotyped cows; g = a vector of random additive genetic effects following N(0, Gs 2 g ), where G = the genomic relationship matrix constructed according to VanRaden [20] and s 2 g = the genomic variance; and e = a vector of random residuals following N(0, Is 2 e ), where I = an identity matrix and s 2 e = residual variance. X and Z were incidence matrices for b and g, respectively.
Age at first calving was grouped into 6 classes (class 1: � 22 months, class 2: 23-24 months, class 3: 25-26 months, class 4: 27-28 months, class 5: 29-30 months, and class 6: � 31 months). Calving age and lactation (L) of the dams were combined to form 21 classes (class 1-6: the same as criteria for age at first calving in L1, class 7: Lactation length for lactation production traits was classified into 16 classes considering equal intervals between 275 to 305 days. The dummy variable indicating in-utero HS was 1 for the respective week for THI � 60, otherwise, a 0 was assigned. We performed separate runs for each week before calving. Alternatively, an interaction model (as introduced by Yao et al. [23]) considering interactions between genotype and in-utero HS, was defined: where g hs = a vector of genotype by HS interaction effects for cows with in-utero HS following N(0, G hs s 2 g hs ), where G hs = the genomic relationship matrix for the cows with in-utero HS and s 2 g hs = the variance of GxE; W = a design matrix allocating phenotypic records to g hs . The remaining effects are the same as described in model VC_wHS.
The genomic heritability (h 2 g ) for each trait was calculated as h 2 g ¼ s 2 g =ðs 2 g þ s 2 e Þ for estimates from model VC_nHS. For estimates from model VC_wHS, the heritability of the common genomic effects (h 2 c ) was h 2 c ¼ s 2 g =ðs 2 g þ s 2 g hs þ s 2 e Þ and the ratio of the variance for genotype by HS interaction effects (r hs ) was r hs ¼ s 2 g hs =ðs 2 g þ s 2 g hs þ s 2 e Þ. Genome wide associations. Subsequently, we estimated main and interaction effects for every SNP via generalized least squares (GLS) equations according to the algorithm as introduced by Yang et al. [24]. The models to estimate main and interaction SNP effects were: For solving the equations of model GWA_nHS, we considered the variance components estimated from model VC_nHS (algorithm is specified below). For solving the equations of model GWA_wHS, we considered the variance components estimated from model VC_wHS (algorithm is specified below).
Matrices and vectors of models for genome wide associations were as follows: x snpi = a vector of centered genotypes calculated as m snpi −2p snpi (m snpi = a vector of genotypes for marker i considering all genotyped animals; p snpi = a vector of allele frequency for marker i); b snpi = a regression coefficient for the i th SNP (the main SNP effect); x interi = a vector of centered genotypes for cows undergoing in-utero HS in the respective week before birth (consecutive runs for the 8 different weeks); and b interi = a regression coefficient for the i th SNP (the interaction effect) under in-utero HS. The remaining effects, vectors and matrices were identical with the effects as specified for models VC_nHS and VC_wHS.
In a self-written R program "GWAInter.R" (in S1 Appendix), we applied GLS to estimate the main and interaction effects for each marker, implying a loop with 41,139 repetitions for all SNP markers. The detailed procedure was: y (the last two fixed effects in the solution from GLS were main SNP and interaction effects for þ Iσ 2 e (the last two variances were the variances for main SNP and interaction effects for marker i, respectively); 2. The test statistic for the i th main SNP effect was: 3. The test statistic for the i th interaction effect was: The inflation factor (λ) was calculated based on the chi-squared statistics as: interi , and w 2 0:5;df ¼1 = 0.4549 (the statistic for a probability of 0.5 from a chi-squared distribution with 1 df). The P-value of the main and interaction effects for each SNP were determined by w 2 snpi and w 2 interi , respectively. Significantly associated SNPs were detected according to the Bonferroni correction calculated as P Bonf = 0.05/(number of SNPs) = 0.05/41,139 = 1.22 × 10 −6 (-log 10 (P Bonf )) = 5.92. Additionally, a less stringent correction was applied and defined as suggestive threshold with P sugg = 0.05/ (number of independent SNPs) = 0.05/3873 = 4.76 × 10 −6 (-log 10 (P sugg )) = 4.89. The number of independent SNPs was calculated based on restrictions for linkage disequilibrium (R 2 � 0.15) for consecutive genomic windows including 500 SNPs (calculated in PLINK 2.0 [25,26]).

Gene annotation
Suggestive SNPs according to P sugg were annotated to potential candidate genes as listed in the Ensembl genome database [27] for main and interaction SNP effects. In this regard, we used the BiomaRt R package [28,29]. Only genes located within a window of ±100 kb of suggestive SNP were considered as potential candidate genes. Afterwards, we submitted the identified potential candidate genes to gene ontology (GO) overrepresentation tests [30], using the GO web-tool [31,32]. The false discovery rate of P < 0.05 was considered to identify overrepresented GO terms for biological processes and reactome pathways. In addition, the windows with suggestive SNPs were mapped with the bovine QTL database through the online Data Analysis Tools [33] to elucidate the phenotypic contributions of the genomic segments. According to Hu et al. [33], the phenotypes were concisely classified into milk, meat and carcass, growth performance, reproduction, exterior and health trait categories.

Variance components
The genomic heritabilities from a model without consideration of in-utero HS (model VC_nHS) as well as common genomic, interaction and residual variances estimated from a model with in-utero HS effects (model VC_wHS) are plotted in Fig 1. Basically, common genomic and residual variances were quite constant for all studied weeks during late gestation. Moreover, both variance components from model VC_wHS were approximately the same when compared to the respective estimates from model VC_nHS. Hence, the genomic heritabilities from model VC_nHS and the common genomic heritabilities from model VC_wHS were 0.13 for TMY, 0.19 for TFP, 0.15 for TFY, 0.24 for TPP, 0.10 for TPY, and 0.08 for TSCS. Minor differences in heritability estimates from model VC_nHS and VC_wHS were observed for lactation production traits. In detail, the h 2 c from model VC_wHS ranged between 0.35 and 0.37 for LMY, between 0.70 and 0.72 for LFP, between 0.34 and 0.36 for LFY, between 0.69 and 0.70 for LPP, and between 0.26 and 0.28 for LPY. The h 2 g from model VC_nHS for the same traits were 0.37, 0.72, 0.36, 0.71, and 0.28, respectively. However, interaction variances varied to a large extent when HS during last eight weeks of pregnancy was considered in consecutive runs, leading to diverse in-utero HS ratios across the eight weeks. Among all traits, LMY displayed the largest r hs in a range from 0.00 (WK6) to 0.07 (WK5). The r hs estimates for TFP were consistently very close to zero in response to climatic impact from all eight last weeks of gestation. For all traits, the interaction variances explained small proportions of the total phenotypic variance. In contrast to small interaction variances reported in our study, environment-specific genomic variances for feed efficiency traits [23] were larger than the common genomic variances from the present studies. However, in this study [23], production and management environments different substantially, because datasets from three countries located in two different continents were merged. Herd specific effects explained very large proportions up to 30% of phenotypic variations in cow production traits in North America [34]. The maximum r hs of 0.07 in our study indicate that the effect of in-utero HS during late pregnancy is smaller than the herd effect due to geographical impact across latitudes. The time period from in-utero HS measurements until the trait recording date is longer for lactation records, but the r hs for lactation production traits were slightly larger than for the respective test-day traits. Nevertheless, the stronger impacts might be due to the accumulative in-utero HS effects during the whole lactation.

Inflation factors
The inflation factors for the main effects from model VC_nHS and VC_wHS ranged from 0.74 for LFP to 1.05 for TSCS (S1 Fig). An inflation factor of value 1 represents sufficient correction for the population structure in GWAS, and a value larger than 1.05 indicates overestimation with inflated false positive results [35]. Consideration of the interaction variance in model GWA_wHS had no impact on inflation factors for the main SNP effects (when compared to inflation factors from model GWA_nHS). However, a larger number of values exceeding 1.05 were observed for the SNP interaction effects from model GWA_nHS (Fig 2). Generally, with regard to model GWA_nHS, the inflation factors indicated more false positives for lactation compared to test-day traits. For test-day traits, less than 30% of all runs had a λ > 1.05, but 67.5% of all runs for lactation milk-production traits displayed such an increased inflation factor. Interestingly, when considering the interaction variances for in-utero HS in model GWA_wHS, inflated results decreased to 12.5% for test-day traits and to 5% for lactation traits. Due to the improvements from model GWA_wHS (i.e., the preferable inflation factors), only significant and suggestive SNPs for main and interaction effects from model GWA_wHS were considered in the ongoing gene annotations.

Significant and suggestive SNPs for main effects
For the main effects, suggestive SNPs from model GWA_nHS overlapped to a large extent with suggestive SNPs from model GWA_wHS. Therefore, only results from GWA_wHS are presented. A total of 380 suggestive SNPs and 237 significant SNPs (including repetitive SNPs for different traits) exceeded the respective thresholds ( Table 3). Number of suggestive and significant SNPs equaled to 171 and 98, respectively, when considering only unique SNPs. Detailed information for the suggestive SNPs is given in S1 Table. Among the eleven analyzed traits, LPP displayed the largest number of suggestive SNPs (91), while only three suggestive SNPs were detected for TSCS, indicating the strong polygenic genetic architecture of TSCS compared to milk production traits [36]. For both test-day and lactation records, results for fat and protein percentages indicated a likewise oligo-genic genetic architecture, as previously reported in genomic studies comparing percentage and yield traits [37]. Hence, the number of suggestive SNPs for milk composition traits including TFP, TPP, LFP, and LPP were always larger than for the respective yield traits TFY, TPY, LFY and LPY. It is well known that fat and protein percentages are controlled by some major genes with large effects, and many unknown genes with small effects. In this regard, the gene DGAT1 on BTA14 strongly contributed to milk fat content [38], and the bovine ABCG2 gene on BTA6 had a major effect on milk protein content [39]. Accordingly, BTA14 comprised the most suggestive SNPs for all traits, apart

PLOS ONE
from TPY and TSCS. SNPs on BTA6 were suggestively associated with six traits, including TMY, TFY, TPP, LMY, LPP and LPY. A larger number of SNPs were detected for lactation than for corresponding test-day milk production traits, reflecting the larger variation for accumulated lactation than for single test-day traits.
In general, for milk production traits, the identified suggestive SNPs for main effects, in particular associated SNPs on BTA3, 5, 6, 10, 14, and 20, were in agreement with the associated SNPs as reported in previous studies [40,41]. A quite large proportion (66.67%) of the suggestively associated SNPs was mapped by significant SNPs influencing milk production traits [40,41]. For TSCS, Meredith et al. [40,42] identified dozens of SNPs, but only one QTL region on BTA24 (spanning from 59.67 Mb to 59.76 Mb) was in close proximity to the SNP detected in our study. Meredith et al. [42] used daughter yield deviations for TSCS as phenotypes for sires in their studies. When considering original phenotypes in a cow population as dependent traits, the number of identified SNPs for TSCS decreased to nine [41].

Significant and suggestive SNPs for interaction effects
With regard to interaction effects, model GWA_nHS generated more suggestive SNPs (51) than model GWA_wHS. This was especially the case for lactation production traits, because the inflation factors for the interaction effects of lactation traits substantially decreased when switching from GWA_nHS to GWA_wHS (Fig 2). The suggestive SNPs from model GWA_nHS contained all suggestive SNPs from model GWA_wHS. Due to the improved inflation factors from GWA_wHS, only suggestive SNPs from GWA_wHS were considered in ongoing gene annotations. The number of suggestive and significant SNPs (31 and 4) for the interaction effects was smaller than for the main effects (Table 4), reflecting the variance proportions for the interaction and common genomic effects, respectively. In contrast to main effects, the number of suggestive and significant SNPs for the interaction effects (estimations from model GWA_wHS) was larger for test-day than for lactation traits. Three from the 31 suggestively associated SNP (model GWA_wHS) were associated with functional traits in Holstein cattle, including ARS-BFGL-NGS-88748 (at 61.48 Mb on BTA19) with body size [43], Table 3. Number of SNPs with suggestive main effects associated with milk production traits and somatic cell score. , had an effect on weight, height at withers, and weaning weight, respectively, in Hanwoo cattle [46]. Hap-map22875-BTA-155031 (at 112.90 Mb on BTA4) was significantly associated with tick resistance in Brazilian Braford and Hereford cattle [47] and ARS-BFGL-NGS-106479 (at 79.75 Mb on BTA11) with intramuscular fat percentage in Australian beef cattle [48]. Furthermore, the Table 4 Holstein haplotype of the gene APOB, detected as a causal mutation for cholesterol deficiency, embraces the SNP ARS-BFGL-NGS-106479 [49]. Interestingly, for the interaction effects, the eight SNPs as reported in the aforementioned studies are more relevant for functional than for primary traits in cattle, although our analyses mainly focused on milk production. Results for the suggestive SNP interactions and their associated traits indicate that robust cattle with superior functionality, e.g., favorable fertility and conformation, might be more tolerant to HS during late pregnancy than cattle with breeding focus on milk production. Mapping all suggestive SNPs for interaction effects with the bovine QTL database revealed that the two most important QTL regions were associated with TPP and TPY.

Potential candidate genes
For test-day milk production traits, a total of 373 unique potential candidate genes are located in regions spanning 100 kb upstream to downstream of suggestive main SNPs from model GWA_wHS (S1 Table). In addition to the DGAT1 gene, the annotated genes CPSF1, FOXH1, ARHGAP39, PPP1R16A, GRINA and MROH1 were frequently associated with milk production traits [41,50,51]. Because of the quite large genetic correlations (� 0.75; in S2 Table) between test-day and corresponding lactation traits, the percentage of overlapping candidate genes was larger than 78%. Interestingly, the genes GHR and ABCG2 were only detected as candidate genes for lactation, but not for test-day milk production traits. Only three candidate genes were located in close distance to suggestive SNPs for TSCS. However, no study reported a relationship between the three genes and TSCS. The gene ABCC9 harbors the most significant SNPs for production traits [41], the gene WDR7 is associated with sperm motility and concentration in Holstein-Friesian bulls [52], and the gene TXNL1 has an effect on moisture content in rainbow trout [53]. The response to oxygen-containing compound and the G-proteifan-coupled receptors signaling pathway were significantly detected based on all candidate genes for main effects. The first pathway represents biological processes inducing changes in cell or an organism activities (in terms of movement, secretion, enzyme production or gene expression), due to an oxygen-containing compound stimulus [54]. The second pathway recognizes a broad array of extracellular mediators including cationic amines, lipids, peptides, proteins and sensory agents [55]. Due to the limited number of suggestive SNPs for interaction effects from model GWA_wHS, a smaller number including 62 unique candidate genes in the vicinity of the SNP locations, were identified (± 100 Kb, S3 Table). In human, the members of ABCA6-like transporters, ABCA6, ABCA8, ABCA9 and ABCA10, are all integrated in cholesterol-related pathways, either directly or based on their dynamic regulation upon cholesterol application [56]. Chemokines, including CCL14 and CCL16, play important roles in regulating bovine endometrial functions during early pregnancy [57]. The genes CCL5, HEATR9, MMP28, GIMAP family and VEGFC activate immune responses [58][59][60][61][62]. Proteins encoded by the genes GAS2L2 and NRIP2 participate in growth process [63,64]. Mutations in the SCN1A and TTC21B genes cause several diseases [65,66]. Sanglard et al. [67] reported a strong correlation between ETNPPL variants and the top 20 differentially expressed genes in beef calves in response to energy restriction during late gestation. GPC4 from the same family of the identified GPC5 and the directly identified gene RASL10B were favorably associated with ETNPPL [67]. The genes LSAMP and GRM8 are involved in emotional and motivational functions [68,69]. The MDGA2 gene encodes novel proteins which regulate neuronal migration and neural development [70]. The GPC5 gene was identified by Naderi et al. [71], who focused on selection signatures in subpopulations of local dual-purpose black and white cattle from Germany. As a result from another selection signature analyses, the gene FOCAD was identified by Alshawi et al. [72] in Iraqi cattle. The gene C8H9orf85 was detected as a candidate gene for milk production traits in water buffalo [73]. The genes TULP3, IQSEC3 and FAM49A were identified as potential candidate genes for milk production traits in US Holstein cattle [41] and for milking speed in French Holstein cattle [74]. The genes RHNO1 and FOXM1 are linked to a variety of cancers and are defined as DNA damage repair regulators [75]. Demetriou et al. [76] explained the role of the FKBP4 gene in recurrent fetal losses in humans. Bourneuf et al. [77] identified de novo deleterious mutations in the FAM189A1 gene. The zinc-finger proteins constitute one of the most abundant groups of proteins in the mammalian genome and are involved in several cellular processes, differentiations of serval tissues, development of several diseases, as well as in tumorigenesis, cancer progression and metastasis formation [78].
From the overrepresentation test, two biological processes, including lymphocyte chemotaxis and monocyte chemotaxis, were excessively represented, referring to all genes from the bovine genome. Lymphocytes contribute to the development of immunity and to allergic inflammation of the lower respiratory tract [79]. Monocytes participate in both innate and adaptive immune responses, because of their phagocytic activity [80]. The identified annotations indicate the biological complexity of in-utero HS on productivity in offspring. Nevertheless, sensitivity to in-utero HS in offspring seems to be correlated with the health status of their dams.
In conclusion, the r hs in the range from 0 to 0.07 indicate quite small effects of HS during late pregnancy on test-day and lactation milk production traits. The accumulative effects due to in-utero HS may explain the generally larger interaction variances for lactation compared to test-day traits. Superiority in inflation factors for the SNP interaction effects from the model GWA_wHS over the model GWA_nHS suggests modelling of interaction variances when applying a GWAS with GxE interactions. Suggestive SNPs for the main effects and their corresponding annotated potential candidate genes reflect findings from previous GWAS focussing on milk production traits. For the SNP interaction effects, eight out of 31 suggestive SNPs are known to effect cattle functionality, indicating a possible favourable correlation between functional traits and heat tolerance. Furthermore, two biological processes as inferred on the basis of the identified 62 candidate genes located in close proximity to the SNPs with suggestive interaction effects, contribute to immune response mechanisms.
Supporting information S1 Table. Suggestive SNPs and potential candidate genes associated with main effects for test-day and lactation production traits. (XLSX) S2 Table. Genetic correlations between test-day and lactation production traits. (DOCX) S3 Table. Potential candidate genes associated with interaction effects for test-day and lactation production traits. a : TMY = first test-day milk yield; TFP = first test-day fat percentage; TFY = first test-day fat yield; TPP = first test-day protein percentage; TPY = first test-day protein yield; TSCS = first test-day somatic cell score; LMY = first lactation milk. (XLSX) S1 Appendix. An R script to perform genome wide associations considering interaction between each marker and a discrete environmental descriptor with two classes. (R) S1 Fig. Inflation factors for main effects from GWAS models with (GWA_wHS) and without genotype by heat stress interaction (GWA_nHS) during the last eight weeks of pregnancy. Dotted line = inflation factor of 1.05; TMY = first test-day milk yield; TFP = first testday fat percentage; TFY = first test-day fat yield; TPP = first test-day protein percentage; TPY = first test-day protein yield; TSCS = first test-day somatic cell score; LMY = first lactation milk yield; LFP = first lactation fat percentage; LFY = first lactation fat yield; LPP = first lactation protein percentage; LPY = first lactation protein yield. (TIF)