Genomic regions, candidate genes, and pleiotropic variants associated with physiological and anatomical indicators of heat stress response in lactating sows

Abstract

Background

Heat stress (HS) poses significant threats to the sustainability of livestock production. Genetically improving heat tolerance could enhance animal welfare and minimize production losses during HS events. Measuring phenotypic indicators of HS response and understanding their genetic background are crucial steps to optimize breeding schemes for improved climatic resilience. The identification of genomic regions and candidate genes influencing the traits of interest, including variants with pleiotropic effects, enables the refinement of genotyping panels used to perform genomic prediction of breeding values and contributes to unraveling the biological mechanisms influencing heat stress response. Therefore, the main objectives of this study were to identify genomic regions, candidate genes, and potential pleiotropic variants significantly associated with indicators of HS response in lactating sows using imputed whole-genome sequence (WGS) data. Phenotypic records for 18 traits and genomic information from 1,645 lactating sows were available for the study. The genotypes from the PorcineSNP50K panel containing 50,703 single nucleotide polymorphisms (SNPs) were imputed to WGS and after quality control, 1,622 animals and 7,065,922 SNPs were included in the analyses.

Results

A total of 1,388 unique SNPs located on sixteen chromosomes were found to be associated with 11 traits. Twenty gene ontology terms and 11 biological pathways were shown to be associated with variability in ear skin temperature, shoulder skin temperature, rump skin temperature, tail skin temperature, respiration rate, panting score, vaginal temperature automatically measured every 10 min, vaginal temperature measured at 0800 h, hair density score, body condition score, and ear area. Seven, five, six, two, seven, 15, and 14 genes with potential pleiotropic effects were identified for indicators of skin temperature, vaginal temperature, animal temperature, respiration rate, thermoregulatory traits, anatomical traits, and all traits, respectively.

Conclusions

Physiological and anatomical indicators of HS response in lactating sows are heritable but highly polygenic. The candidate genes found are associated with important gene ontology terms and biological pathways related to heat shock protein activities, immune response, and cellular oxidative stress. Many of the candidate genes with pleiotropic effects are involved in catalytic activities to reduce cell damage from oxidative stress and cellular mechanisms related to immune response.

Background

Global warming has become a major concern for the long-term sustainability of livestock production. The ten warmest years on record since 1880 have all occurred after 2010, and the nine years from 2014 to 2022 are the nine warmest years on record [1]. These higher temperatures have a direct impact on animal production systems due to increased incidence of heat stress (HS) [2]. The absence of functional sweat glands in pigs limits their ability to use evaporative heat loss through sweating, reducing their thermoregulation ability in hot environments [2, 3]. Furthermore, intensive genetic selection for greater productive and reproductive performance has been associated with an increased metabolic heat production, making the animals more sensitive to high environmental temperatures [4]. In this context, breeding for greater heat tolerance could improve the animals’ ability to cope with HS conditions.

Breeding pigs for improved resilience to HS is possible as there is genetic variability in their response to HS [5,6,7,8]. One alternative for genetically identifying animals that perform better under HS conditions is through the use of reaction norm approach that model the genetic merit of an individual over an environmental gradient [9]. In general, reaction norm models are based on variability in performance traits across a continuous environmental gradient. Estimated breeding values for heat tolerance (slope of the reaction norms) based on this approach tend to be unfavorably correlated with the level of production (e.g [7]). Therefore, identifying physiological, behavioral, and anatomical indicators of HS response that could be incorporated in selection indexes to breed for improved climatic resilience is promising in modern pig breeding programs [10]. Breeding for indicators of HS response could capture additional biological mechanisms involved in HS response and enable improvements in climatic resilience without compromising productivity [11]. Johnson et al. [10] described several traits as HS response indicators, such as skin and automatically-recorded vaginal temperatures, and respiration rate, and Freitas et al. [11] showed that these traits are heritable, i.e., they can be improved through genetic selection. Understanding the genetic background of these novel traits is a key step before including them in a breeding scheme. This includes the detection of quantitative trait loci (QTL) and the identification of candidate genes and biological pathways controlling the phenotypic expression of these HS indicator traits.

Genome-wide association studies (GWAS) enable the identification of genomic markers statistically associated with the traits of interest. GWAS is based on the principle of identifying genomic markers in linkage disequilibrium (LD) with the QTL affecting the trait of interest [12]. Various studies have reported genomic regions associated with heat tolerance based on reaction norm models for reproductive traits [6] and productive traits such as hot carcass weight [8, 13]. Kim et al. [14] did not identify significant markers associated with respiration rate, rectal temperature, and skin temperature in crossbred gilts, however, they did report some genes known to be involved in physiological adaptation to general stressors. However, there is still a lack of studies aiming to identify candidate genes associated with physiological indicators of HS, such as vaginal temperature, skin temperature, and respiration rate in pigs using many phenotypic records and genomic markers.

Heat stress indicator traits may share metabolic pathways. A single locus can be involved in the expression of more than one trait, which is a phenomenon known as pleiotropy [15]. Pleiotropy occurs when a single causal variant or multiple variants within the same gene (or region) are associated with different phenotypes. Furthermore, pleiotropy also occurs when one phenotype causes a second phenotype, and a genetic variant is directly associated with the initial phenotype [16]. Pleiotropy is one of the causes of genetic correlation between traits. In this context, Freitas et al. [11] reported positive and moderate to high genetic correlations among physiological indicators of HS response such as skin temperature, vaginal temperature, panting score, and respiration rate. Identifying genes with pleiotropic effects and the genetic correlation between traits allows breeders to understand the biological results of the correlated response in multiple-trait selection. In this context, the primary objectives of this study were to identify genomic regions associated with physiological and anatomical indicators of HS response in lactating sows using imputed whole-genome sequence (WGS) data and identify pleiotropic variants. In addition, candidate genes, gene ontology (GO) terms, and metabolic pathways were identified in the significant genomic regions associated with the HS responses of lactating sows.

Results

Overview

We performed GWAS for 18 indicators of HS response and anatomical traits in crossbred lactating sows (Large White x Landrace) and found candidate genes and pleiotropic variants affecting these traits. The thermoregulatory indicator traits evaluated were ear skin temperature (T_ES), shoulder skin temperature (T_SS), rump skin temperature (T_RS), tail skin temperature (T_TS), respiration rate (RR), panting score (PS), vaginal temperature automatically measured every 10 min (T_Vall), the average of the six records per hour corresponding to 0800, 1200, 1600, and 2000 h during four days (T_V4days), and single records corresponding to 0800, 1200, 1600, and 2000 h measured at the first day of collection (T_V8h, T_V12h, T_V16h, T_V20h, respectively). The anatomical traits evaluated were hair density (HD), body size (BS), body condition score (BCS) using a sow caliper (BCS_cal), visual BCS (BCS_vis), ear length (EL), and ear area (EA). The variance components and genetic parameters for the traits evaluated in this study have been previously reported by Freitas et al. [11]. A description of all the studied traits and their corresponding heritability estimates, as previously reported by Freitas et al. [11], is presented in Table 1. The significant SNPs were identified based on their respective p-values, with adjustment for multiple testing using the Bonferroni correction at a significance level of α = 0.05 for genome-wide significance. To take into account the dependence among tests due to linkage disequilibrium, α was divided by the number of independent chromosome segments ($Me$).

Table 1 Description of the evaluated traits and heritability estimates reported by Freitas et al. [11]

Genome-wide association and post-GWAS indicators of heat stress response in lactating sows

The medium-density genotypes were imputed to WGS using the Swine Imputation (SWIM) Server [17] and the imputed genotypes were used to perform GWAS for all traits included in the study. A total of 1,388 SNPs across 16 chromosomes were associated with T_ES, T_SS, T_RS, T_TS, RR, PS, T_Vall, T_V8h, HD, BCS_cal, and EA (Table S1). T_SS and T_RS are influenced by many significant genomic regions located on many chromosomes (Fig. 1.b-c), indicating that these traits are highly polygenic. For T_ES, significant peaks located on chromosomes 5, 7, 12, and 16 were identified (Fig. 1.a), while for T_TS, we observed significant peaks located on chromosomes 4, 6, and 15 (Fig. 1.d). Two significant genomic regions located on chromosomes 7 and 8 were identified for RR (Fig. 1.e) and one significant region on chromosome 9 was found to be associated with PS (Fig. 1.f).

Fig. 1

Manhattan plots for skin temperatures, respiration rate, and panting score using whole-genome sequence data. -log10(p-values) for ear skin temperature (T_ES; a), shoulder skin temperature (T_SS; b), rump skin temperature (T_RS; c), tail skin temperature (T_TS; d), respiration rate (RR; e), and panting score (PS; f). Genome-wide significance level is shown in a red line

The Manhattan plots for vaginal temperature traits are shown in Fig. 2. One significant marker was observed on chromosome 1 (Fig. 2.a) for T_Vall, and one significant peak (228 SNPs) on chromosome 2 for T_V8h (Fig. 2.c). No significant SNPs were observed for T_V4days, T_V12h, T_V16h, and T_V20h after Bonferroni correction for multiple tests (Fig. 2). Figure 3 shows the results of GWAS for anatomical traits. Two significant SNPs located on chromosomes 2 and 7 and a significant peak (13 SNPs) on chromosome 15 were identified for HD (Fig. 3.a). Four significant SNPs located on chromosome 13 were found to be associated with BCS_cal (Fig. 3.c), while three significant SNPs on chromosome 7 were identified to be associated with EA (Fig. 3.e). No significant SNPs were observed for BS, BCS_vis, and EL after multiple testing correction (Fig. 3). The Q-Q plots for the GWAS analyses for all evaluated traits are presented in Figures S1-S3 (Additional File 1). The lambda values (genomic inflation factor) ranged from 1.00 to 1.32.

Fig. 2

Manhattan plots of GWAS for vaginal temperatures using whole-genome sequence data. -log10(p-values) for all measures (every 10 min) of vaginal temperatures for four days (TV_all; a), four-time measures of vaginal temperatures for four days (T_V4days; b), vaginal temperature measured on the first day at 8:00 (T_V8h; c), at 12:00 (T_V12h; d), at 16:00 (T_V16h; e), and at 20:00 (T_V20h; f). Genome-wide significance level is shown in a red line

Fig. 3

Manhattan plots of GWAS for anatomical traits using whole-genome sequence data. -log10(p-values) for hair density (HD; a), body size (BS; b), body condition score using a sow caliper (BCS_cal; c) and visual (BCS_vis; d), ear area (EA; e), and ear length (EL; f). Genome-wide significance level shown in a red line

Genomic windows of 100 kb on each side of significant genomic markers were used to identify candidate genes associated with each trait, as shown in Table S2. A total of 26, 151, 102, 35, 63, eight, two, 93, four, four, and 11 genes were identified to be associated with T_ES, T_SS, T_RS, T_TS, RR, PS, T_Vall, T_V8h, HD, BCS_cal, and EA, respectively. The functional enrichment analyses enabled the identification of 20 GO terms (five biological processes, four cellular components, and 11 molecular functions), and 11 metabolic pathways (Table 2). The complete list of all GO and metabolic pathways identified, including those considered to be statistically significant and suggestive, are shown in Table S3.

Table 2 Significant gene ontology terms and pathways for heat stress related traits

Pleiotropic effects

As traits are correlated and pleiotropic effects are expected, a multiple-trait analysis, as proposed by Bolormaa et al. [18] was also performed. First, the traits were decorrelated using Cholesky transformation [19, 20] and grouped into seven categories: skin surface temperature traits (T_ES, T_SS, T_RS, and T_TS), vaginal temperatures (T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, and T_V20h), indicators of animal temperature (i.e., all skin surface and vaginal temperature traits), respiration rate traits (RR and PS), thermoregulatory indicators (i.e., all skin surface temperatures, vaginal temperatures, and respiration traits), anatomical traits (BCS_cal, BCS_vis, HD, BS, EA, and EL), and a combination of all traits.

We identified one significant SNP with potential pleotropic effect for the skin temperature group of traits, six significant variants for the vaginal temperatures group of traits, six significant markers for the group of traits related to animal temperature, one significant marker for respiration rate traits, four significant markers for the thermoregulatory traits, and 21 significant markers for anatomical traits (Fig. 4). The significant marker observed for the skin temperatures group was also observed for the group containing all traits. The same significant region on chromosome 7 was observed for vaginal temperature and animal temperature groups, which was also observed for the thermoregulatory traits group.

Fig. 4

Manhattan plots for multiple-trait analysis for different categories of traits. (a) The category skin temperatures included the traits ear skin temperature (T_ES), shoulder skin temperature (T_SS), rump skin temperature (T_RS), tail skin temperature (T_TS). (b) The category vaginal temperatures included all measures (every 10 min) of vaginal temperatures for four days (TV_all), four-time measures of vaginal temperatures for four days (T_V4days), vaginal temperature measured on the first day at 8:00 (T_V8h), at 12:00 (T_V12h), at 16:00 (T_V16h), and at 20:00 (T_V20h). (c) The category animal temperature included the traits T_ES, T_SS, T_RS, T_TS, T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, and T_20h. (d) The category respiration traits included respiration rate (RR) and panting score (PS). (e) The category thermoregulatory traits included the traits T_ES, T_SS, T_RS, T_TS, T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, T_20h, RR, and PS. (f) The category anatomical traits included the traits hair density (HD), body size (BS), body condition score using a sow caliper (BCS_cal) and visual (BCS_vis), ear area (EA), and ear length (EL). (g) The category all traits included the traits (T_ES, T_SS, T_RS, T_TS, T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, T_20h, RR, PS, HD, BS, BCS_cal, BCS_vis, EA, and EL). Genome-wide significance level shown in a red line

Seven, five, six, two, seven, 15, and 14 genes were identified for the trait groups of vaginal temperature, animal temperatures, respiration traits, thermoregulatory traits, anatomical traits, and all traits, respectively (Table 3). Five common genes (ENPP5, ENPP4, RCAN2, U6, and CLIC5) were annotated for the vaginal temperatures, animal body temperatures, and thermoregulatory traits groups. Only one significant (p = 0.0269) GO term of biological process (GO:0007605 - Sensory perception of sound) for the animal temperatures group, involving the annotated genes ADGRV1 and CLIC5 was identified.

Table 3 Candidate genes with pleiotropic effects for physiological and anatomical indicators of heat stress response

Discussion

Genome-wide association for heat stress indicators

The traits evaluated in this study include physiological and anatomical indicators of heat stress response in lactating sows. The physiological indicators are response to mechanisms used by the animal to dissipate heat, such as increasing body temperature and respiration rate [21, 22]. The anatomical indicators are anatomical characteristics that affect the ability of the animal to dissipate heat, such as the animal’s surface area, hair density, and body mass [23]. The heritability estimates for the traits evaluated in this study were previously reported by Freitas et al. [11]. Skin temperatures (i.e., T_ES, T_SS, T_RS, and T_TS), RR, and PS exhibited low heritability estimates, ranging from 0.04 to 0.06; vaginal temperatures (i.e., T_Vall, T_V4days, T_V8h T_V12h, T_V16h, and T_V20h) were moderately heritable, with estimates ranging from 0.15 to 0.29 [11]; and, anatomical traits (i.e., BCS_cal, BCS_vis, HD, BS, EA, and EL) showed moderate to high heritability estimates, ranging from 0.25 to 0.40 [11].

The statistical power of GWAS is impacted by various factors, including trait heritability [24, 25]. A higher heritability increases statistical power, as heritability is a key determinant of effect size [24]. Consequently, as heritability increases, the effect size of each causal variant can also increase, thereby increasing the probability of identifying causal variants [24, 25]. Additionally, the performance of GWAS can be improved by using higher density genotyping. Several studies have shown that performing GWAS using whole-genome sequence (WGS) data can improve QTL detection [26,27,28] because it contains most causal mutations in the genome and can enable the identification of variants that are highly linked and closer to the causal mutations as compared to genotyping platforms with less SNPs. However, sequencing a large number of animals is still expensive. A more cost-effective approach is the use of genotype imputation from low-density to WGS data [28, 29].

In this study, the genomic inflation factor across the 18 evaluated traits ranged from 1.00 to 1.32, which suggested that the results might be slightly inflated for some traits. Although the inflation factor may indicate small population structure effects, Yang et al. [30] showed that significant inflation of test statistics is expected under polygenic inheritance even when there is no population stratification. Additionally, the genomic inflation factor reflects the trait heritability and the number of causal variants [30]. The range of observed values are within those observed in the literature for other complex traits. Multiple markers located on sixteen chromosomes, each with relatively small effect sizes (Table S1), were associated with the physiological and anatomical indicators of HS response in lactating sows, suggesting that these traits are highly polygenic. For some of the evaluated traits, no significant markers were identified, which could be due to the study sample size or the polygenic nature of the study traits. Polygenic traits are influenced by numerous loci with small effects on trait expression, reducing the likelihood of identifying significant markers in association studies [31]. Furthermore, identifying markers associated with such traits requires larger sample sizes and more phenotypic records to enhance the power of the association analyses [31]. The amount of data used in this study might not be sufficient to identify associated markers with some of these traits (e.g., T_V4days, T_V12h).

Some genomic regions were commonly found to be significant for the skin surface temperature traits, which was already expected since these traits are highly genetically correlated (from 0.78 to 0.91) [11]. Additionally, Johnson et al. [10] reported moderate to high phenotypic correlations for these same traits (from 0.56 to 0.76). As these traits are highly genetically correlated, genetically selecting for skin temperature measured at any of the locations evaluated in this study should result in similar genetic progress on all of them due to indirect genetic responses.

A significant genomic region located on chromosome 2 was associated with vaginal temperature measured at 800 h. The genes harboring this region included RASA1 (Ras GTPase-Activating Protein 1), a regulator of blood vessel and lymphatic vessel growth in adult mice and humans. RASA1 induced-deficiency results in a lymphatic vessel disorder characterized by hyperplasia and leakage [32], which can affect the ability to vasodilate and dissipate excess body heat through the skin. Another gene found in the significant genomic region for T_V8h was LYSMD3 (LysM Domain Containing 3). This gene plays a crucial role in the innate immune response by recognizing chitin and fungal spores, mediating the production of cytokines, and promoting early inflammatory responses in human lung epithelial cells [33]. For T_Vall, the candidate gene AKAP7 (A-kinase Anchoring Protein 7), which encodes a protein kinase A-binding scaffolding molecule, was identified. This gene is crucial in cellular regulatory pathways during viral infections and influences neuroprotection by enhancing mitochondrial networks and inhibiting apoptosis under glutamate-induced oxidative stress, being reported as a cytoprotective factor, potentially preventing tissue damage during viral infections [34,35,36]. Due to AKAP7’s cytoprotective functions, it might also contribute to preventing tissue damage in HS conditions. Investigating the same population as in this study, Wen et al. [37] used random regression models to analyze all automatically-recorded vaginal temperature records as longitudinal traits. Wen et al. [37] identified two important genomic regions on chromosomes 10 and 16, containing genes associated with immunity, protein transport, and energy metabolism. However, in this study, we did not find these same genomic regions and genes, which might be due to different modeling approaches and the small effect size of the identified genomic regions.

It is well-known that heat shock proteins (HSP) play an important role in cell survival under stress conditions, especially HS [38,39,40]. HSP protects cells from damage by binding to unfolded or denatured proteins, preventing their aggregation and promoting their refolding [41, 42]. In addition, HSP helps to restore protein structure and function, ensuring that essential processes can continue despite the stress [41, 42]. HSP also acts in the degradation of proteins that cannot be refolded properly and can interact with other proteins and modulate their activity or stability, influencing various cellular processes such as cell cycle progression and apoptosis [41]. HSP activities have been reported to minimize HS-induced cellular death and improve thermotolerance in livestock [43] and to regulate the immune response, including stimulatory and regulatory functions [44]. In this study, we found genomic regions associated with T_SS, T_RS, and T_TS overlapping with HSBP1L1 (Heat Shock Factor Binding Protein 1 Like 1) and regions associated with EA overlapping with the HSP90AB1 (Heat Shock Protein 90 Alpha Family Class B Member 1) gene. HSBP1L1 is presumed to suppress the heat shock factor transcription under stress, suggesting that HSBP1L1 may regulate cellular responses to stress and cellular heat acclimation [45]. HSP90AB1 encodes a member of the heat shock protein 90 family, a crucial protein for signal transduction, protein folding, and degradation, playing a significant role in maintaining protein stability and promoting cell survival, particularly under cellular stress [46]. Genetic variations in HSP90AB1 have been linked to stress-induced mortality, heat tolerance, and productivity in animals [46].

One of the main responses to HS and mechanisms used for pigs to maintain their body temperature is increased RR, allowing them to dissipate heat through evaporation [4, 47, 48]. In the genomic windows around the two significant regions associated with RR we identified 63 genes, including HSPA1L (Heat Shock Protein Family A -Hsp70 - Member 1 Like), which encodes a 70 kDa heat shock protein crucial for cellular and molecular responses to HS in livestock. HSPA1L ensures proper protein folding, unfolding, and refolding, having protective functions alongside other HSPs during HS events [49]. In the functional analyses, which includes all identified candidate genes for RR, a biological process related to protein refolding and cellular response to unfolded protein was revealed based on the GO analyses. Molecular functions of misfolded protein binding and amino acid binding were also identified, and these functions are linked to HSP activity.

It is already known that HS affects the immune response in pigs [50,51,52]. The genes found to overlap with the significant regions observed for T_ES include BTBD17 (BTB domain containing 17), PTPRE (Protein Tyrosine Phosphatase Receptor Type E), and ARNT2 (Aryl Hydrocarbon Receptor Nuclear Translocator 2). BTBD17 has been reported to be involved in the negative regulation of viral genome replication and viral response [53, 54]. PTPRE has been shown to repress M2 macrophage activation and may have a critical role in regulating inflammatory responses and local homeostasis, particularly in the context of macrophage activation and function [55]. ARNT2 regulates immune-related genes, phagocytic cell differentiation, and lymphocyte activity [56]. A polymorphism in ARNT2 has been linked to reduced macrophage phagocytic activity, impacting immunity, platelet activation, and phagocytosis in pigs [56]. Additionally, we also found that the NR1D1 (Nuclear Receptor Subfamily 1 Group D Member 1) gene harbors genomic markers associated with T_SS and T_RS. NR1D1 acts as a transcriptional repressor and is required for establishing and maintaining body temperature rhythm in a manner adaptable to environmental demands, allowing mammals to adapt to environmental temperature changes [57].

Previous studies have reported that some anatomical characteristics, such as the animal’s surface area, hair density, and body mass, can affect the ability of the animal to dissipate heat [23]. In this study we identified that regions associated with EA overlapping with the HSP90AB1 (Heat Shock Protein 90 Alpha Family Class B Member 1) gene that encodes a member of the heat shock protein 90 family, as previously discussed. Additionally, we found that regions associated with HD overlap with the gene DOCK10 (dedicator of cytokinesis 10). DOCK10 is part of the dedicator of cytokinesis (DOCK) family proteins, which are expressed at varying levels in lymphocytes and have diverse effects on immune functions [58]. DOCK10 influences the B cells development, which play a crucial role in the adaptive immune response by producing antibodies [58].

In the functional genomic analyses, the molecular functions of intracellular calcium activated chloride channel activity and chloride channel activity were found to be associated with T_ES. The chloride anion plays a key role in regulating cellular homeostasis, affecting diverse cellular functions, including reactive oxygen species levels [59]. Oxidative stress stimulates calcium activated chloride channels [60,61,62], and promotes cellular and mitochondrial oxidative damage [63]. In addition, the molecular function of DNA replication origin binding was associate with T_SS and T_RS and the biological process of centrosome cycle was associated with T_TS. These two GO terms are related to cell division and are affected by HS. HS can induce partial DNA re-replication and centrosome amplification in early S-phase cells [64]. Although HS causes centrosome degradation, HSP70 protects centrosome proteins from denaturation [65]. In the functional genomic analysis for RR, we identified two pathways: the pathway of Staphylococcus aureus infection and the pathway of complement and coagulation cascades. The latter is a proteolytic cascade in blood plasma and serves as a mediator of innate immunity, a nonspecific defense mechanism against pathogens [66,67,68].

Candidate genes with pleiotropic effects for heat stress-related traits in lactating sows

Pleiotropy is a well-known effect that occurs when a locus influences more than one trait [16, 69]. In cases of pleiotropy, if a gene has effects on more than one trait, the genetic variation associated with one trait may be correlated with the genetic variation of another trait, resulting in a genetic correlation between the two traits. Freitas et al. [11] reported high genetic correlations among skin temperature traits (> 0.78) and among vaginal temperature traits (> 0.75). Similarly, high genetic correlations were observed for RR with PS (0.87), with BCS measures (0.92), and with ear length and ear area (0.95). Additionally, low to moderately positive genetic correlations were observed between skin temperatures and vaginal temperatures (ranging from 0.25 to 0.76). Low genetic correlations were estimated between vaginal temperatures and BCS (from − 0.01 to 0.06), while all other trait combinations were lowly genetically correlated [11].

Identifying candidate genes with pleiotropic effects is important to understand how a single gene can affect multiple traits and the common biological pathways involved in the expression of these traits. In GWAS studies, it is possible to identify the same regions or genes affecting different traits, which is termed cross-phenotype associations, however these associations might not be pleiotropy [16]. Solovieff et al. [16] defined three categories of pleiotropy: biological pleiotropy, mediated pleiotropy, and spurious pleiotropy. Biological pleiotropy refers to a genetic variant that directly influences more than one trait [16, 70]. Mediated pleiotropy is when a variant directly affects a first phenotype that affects the second phenotype [16, 70]. Spurious pleiotropy can occur due to some sources of bias, identifying genes that are not truly pleiotropic [16, 70].

In this study, to identify the pleiotropy among the traits, we grouped them as skin surface temperature traits (T_ES, T_SS, T_RS, and T_TS), vaginal temperatures (T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, and T_V20h), animal body temperatures (i.e., all skin surface and vaginal temperature traits), respiration traits (RR and PS), thermoregulatory indicators (i.e., all skin surface temperatures, vaginal temperatures, and respiration traits), anatomical traits (BCS_cal, BCS_vis, HD, BS, EA, and EL), and a combination of all traits.

A common genomic region located on chromosome 7 was associated with vaginal temperature, animal body temperatures, and thermoregulatory traits groups. This region comprised the genes RCAN2 (Regulator of Calcineurin 2), ENPP5 (Ectonucleotide Pyrophosphatase/Phosphodiesterase Family Member 5), ENPP4 (Ectonucleotide Pyrophosphatase/Phosphodiesterase 4), CLIC5 (Chloride Intracellular Channel 5), and U6 (U6 small nuclear RNA). ENPP5 and ENPP4 belong to the ectonucleotide pyrophosphatase/phosphodiesterase (ENPP) family of genes encoding a group of proteins that have enzymatic activity involved in the hydrolysis of nucleotides [71]. Based on their sequence homology, ENPP4 and ENPP5 have been described as putative ENPP ectoenzymes [72]. The ectoenzymes are found on leukocytes and endothelial cells and play a crucial role in regulating leukocyte trafficking [73]. RCAN2 is a gene that encodes a protein involved in the regulation of the calcineurin signaling pathway. This protein is involved in many physiological processes by binding to the catalytic domain of calcineurin A and inhibiting calcineurin-mediated nuclear translocation of the transcription factor NFATC1 (nuclear factor of activated T cells 1) [74, 75]. Inactivation of the RCAN2 gene in mice reduced age- and diet-induced obesity by causing a reduction in feed intake [76]. CLIC5 encodes a member of the chloride intracellular channel (CLIC) family of chloride ion channels that regulates chloride transport and interacts with the cortical actin cytoskeleton in polarized epithelial cells [77]. U6 is a small nuclear RNA (snRNA) involved in catalytic steps of pre-mRNA splicing [78]. This common genomic region identified for vaginal temperatures, animal body temperatures, and thermoregulatory traits groups might be relevant for thermoregulatory indicators since the genes annotated in this region are involved in immune response and in catalytic activities to reduce cell damage caused by oxidative stress.

The candidate genes MCTP1 (Multiple C2 And Transmembrane Domain Containing 1), and ELL2 (Elongation Factor for RNA Polymerase II 2), both located on chromosome 2, were candidate genes for animal temperature traits that also presented pleiotropic effects. MCTP1 regulates endocytic recycling in central nervous system neurons and synapses, and possibly cellular vesicle retrieval and oxidative stress [79], which might influence HS response since it also causes cellular oxidative stress. ELL2 has been reported to contribute to prostate homeostasis and potentially acts as a tumor suppressor in the development and progression of prostate cancer cells [80, 81]. Additionally, it regulates plasma cell gene expression, including immunoglobulin production in plasma cells [82]. Therefore, it is possible that the ELL2 gene may influence immune response when the sows are under HS conditions.

We also annotated the protein-coding genes NTN1 (Netrin-1) and STX8 (Syntaxin 8) located on chromosome 12 for the thermoregulatory traits group. NTN1 belongs to a family of laminin-related secreted proteins and contributes to the regulation of inflammation by influencing cell migration [83, 84]. NTN1 is involved in immune response modulation, regulation of inflammation, and leukocyte infiltration [83]. STX8 has been identified as a key player in the efficient sorting and trafficking of cytotoxic molecules to functional lytic granules in cytotoxic T lymphocytes in humans [85]. Both NTN1 and STX8 contribute significantly to the immune response, particularly in combating pathogens through the actions of macrophages and cytotoxic T cells.

For the RR group of traits, the gene AKAP7 located on chromosome 1 was annotated. AKAP7 is a gene encoding a protein kinase A-binding scaffolding molecule associated with a cellular regulatory pathway during viral infections [34, 35]. In addition, this gene promotes neuroprotection by enhancing mitochondrial networks and blocking apoptosis in glutamate-induced oxidative stress [36]. According to Gusho et al. [34], AKAP7 may serve as a cytoprotective factor and potentially prevent tissue damage during viral infections. Several candidate genes with pleiotropic effects for physiological indicators of HS response in lactating sows were associated with immune response, supporting the already reported effects of HS on immune responses [50,51,52].

The gene TMEM161B (Transmembrane Protein 161B), located on chromosome 2, was found to be a candidate gene with pleiotropic effects for skin temperature traits. In zebrafish and mice, TMEM161B is linked with cardiac rhythm, preventing arrhythmic calcium oscillations, and ensuring proper blood flow for the developing embryo [86]. Moreover, studies indicate that increasing vasodilatation elevates skin blood flow, increasing skin temperature [87], which facilitates heat dissipation. In humans, TMEM161B has been reported to regulate cerebral cortical gyration, Sonic Hedgehog signaling, and ciliary structure in the developing central nervous system [88]. This gene might be important under HS conditions, where the sympathetic nervous system coordinates adjustments in cardiac rhythm and skin blood flow as response to thermal challenges [89, 90]. The rate of skin blood flow is primarily determined by body temperatures, and above an internal body temperature threshold for vasodilation, the increase in skin blood flow is proportional to the increase in internal temperature managing the impact of increased skin blood flow and heat transfer requirements on the circulatory system [90, 91].

For the anatomical traits group, the candidate genes with pleiotropic effects found included ZMYND11 (Zinc finger MYND-type containing 11) and DIP2C (Disco Interacting Protein 2 homolog C) located on chromosome 10, and GDNF (Glial Cell Derived Neurotrophic Factor) and CDH18 (Cadherin 18) on chromosome 16. ZMYND11, a multidomain protein, modulates RNA polymerase activity and regulates splicing. This gene was differentially regulated in response to thermal stress in sockeye and pink salmon [92] and in chicken pituitary under heat stress [93]. ZMYND11 has been associated with neurodevelopmental disorders, global developmental delay, seizures, and hypotonia, and it is a potential causative gene in complex neurodevelopmental phenotypes [94,95,96]. While the biological functions of DIP2C are not entirely clear, it may regulate neuronal differentiation and early embryonic nervous development [97], as well as overall brain development and function [98]. GDNF is essential for maintaining neuronal morphological and neurochemical phenotypes, and protecting dopaminergic neurons from toxic damage [99]. It has several positive effects on the nervous system, including viability, proliferative activity, and migratory ability of cells [100]. CDH18 has been associated with neuropsychiatric disorders in humans [101] and tumor suppressors [102, 103]. The candidate genes with pleiotropic effect for anatomical traits influence nervous system functions. The response to heat stress involves complex interactions between the brain and various neurotransmitter systems [104]. Furthermore, prolonged exposure to HS conditions can negatively affect central nervous system, leading to cognitive and neurological sequelae [105].

Implications and limitations

Understanding the genetic background of physiological and anatomical traits related to HS response is valuable for breeding programs aimed to improve HS resilience in pigs. In this study we identified many genomic regions associated with indicators of HS response in lactating sows, but the QTL presented small effect sizes indicating that the traits are highly polygenic. The candidate genes for indicators of HS response evaluated are associated with catalytic activities to reduce cell damage from oxidative stress and cellular mechanisms related to immune response. Furthermore, gene ontology terms and pathways involved in heat shock protein activities, immune response, and oxidative cellular stress were found, which inform biological mechanisms related to heat tolerance and resilience in pigs. Additionally, in this study we identified variants with potential pleiotropic effects and influence in multiple indicators of HS response in lactating sows. This information allows breeders to design programs that could maximize genetic gain across multiple traits for more comprehensive improvements in animal populations. The utilization of imputed whole-genome sequencing (WGS) enabled the identification of significant variants that were absent on the medium-density SNP array. Additionally, pleiotropic variants can be filtered and included in the SNP array [106], which would benefit breeding programs focusing on selecting animals with greater climatic resilience. However, further investigations are required to assess pleiotropy between HS response indicators and traits related to production and reproduction. This evaluation is crucial to determine the importance of these variants for other relevant traits.

The use of crossbred animals represents a limitation of this study, as different genomic regions may be associated with the same indicators of HS response in purebred populations. This discrepancy arises because allelic frequencies vary across populations and breeds. Additionally, the dataset is not large, and the animals are from a single farm. Therefore, caution should be taken when extrapolating findings from crossbred animals to purebred populations, as the genetic background of the populations can significantly influence the observed associations between genomic regions and HS response indicators.

In future studies, the incorporation of economically important traits associated with reproductive and productive performance would allow the understanding of pleiotropy between these traits and those related to HS response. This expanded knowledge would provide valuable insights for creating selection indexes, enabling the identification of animals that are not only more resilient to HS but also capable of maintaining high productivity.

Conclusions

In this study we identified numerous genomic regions significantly associated with indicators of HS response in lactating sows. The significant variants have small effects, indicating the highly polygenic nature of the traits evaluated. The candidate genes for physiological and anatomical indicators of HS response are associated with gene ontology terms and pathways involved in heat shock proteins activities, immune response, and oxidative cellular stress. The indicators of HS response exhibited pleiotropic effects, suggesting that multiple genomic regions influence various physiological and anatomical indicators of HS response. Many candidate genes with potential pleiotropic effects are associated with catalytic activities to reduce cell damage from oxidative stress and cellular mechanisms related to immune response were identified, providing valuable insights for breeding programs aimed at improving HS resilience in pigs. The findings have significant implications for advancing breeding programs aimed at improving heat tolerance in pigs.

Materials and methods

Animals and genotypes

The data used in this study was previously described by Johnson et al. [10] and Freitas et al. [11]. In summary, the phenotypes were collected from 1,645 multiparous lactating sows (Large White x Landrace cross) at a commercial sow farm located in Maple Hill, North Carolina, USA. A total of 1,639 sows were genotyped using the PorcineSNP50K (50,703 SNPs) Bead Chip (Illumina, San Diego, CA, USA). The quality control applied to the genotypes consisted of removing animals with a call rate lower than 0.80 and SNPs with a call rate lower than 0.90, and SNPs located on non-autosomal chromosomes. After quality control, 36,258 SNPs for 1,622 animals remained for further analyses. The genotypes were imputed to WGS containing 34,615,361 SNPs using the Swine Imputation (SWIM) Server [17]. The reference population consisted of 2,259 animals from various breeds, including 615 Landrace and 543 Large White, and 25 Landrace x Large White crossbred animals [17].

After genotype imputation, additional quality control was performed, removing variants with impute scores lower than 0.85, animals and SNPs with missing rates greater than 0.10, minor allele frequency (MAF) lower than 0.05, and extreme p-values for Hardy-Weinberg equilibrium tests lower than 1 × 10^− 15. Data from 7,065,922 SNPs for 1,622 animals remained for subsequent analyses. LD-based pruning was performed considering a pairwise correlation coefficient threshold of 0.90 between markers, and a window size of 10 markers, with a step size of five markers, retaining 1,710,312 SNPs. Some advantages of performing LD pruning are to reduce redundancy and confounding factors that may interfere with statistical analyses and reduce the number of tests performed improving the power of the analyses [107]. Genotyping data was pre-processed using the Plink software [108].

Phenotypic records

The description of the phenotypes evaluated in this study are presented in Table 1. In summary, the thermoregulatory indicators of HS response phenotypes collected were respiration rate (RR), ear skin temperature (T_ES), shoulder skin temperature (T_SS), rump skin temperature (T_RS) tail skin temperature (T_TS), vaginal temperatures, and panting score (PS). RR and skin surface temperatures were collected at 0800, 1200, 1600, and 2000 h daily during a period of four consecutive days. Vaginal temperature was automatically recorded in 10 min intervals over four days using calibrated thermochron temperature recorders. This data was used to compute six traits in this study: whole data measured at every 10 min (T_Vall), the average of the six records per hour corresponding to 0800, 1200, 1600, and 2000 h during four days (T_V4days), and single record corresponding to 0800, 1200, 1600, and 2000 h measured at the first day of collection (T_V8h, T_V12h, T_V16h, and T_V20h, respectively). Panting score was collected during four consecutive days at 15:30 h and scored from 0 to 3 (0: if it was with closed mouth and normal breathing; 1: closed mouth and rapid breathing; 2: open mouth and rapid breathing; 3: open mouth and rapid breathing with obvious salivation).

We also collected phenotypes for sow anatomical traits that might influence on the ability of the animal to dissipate metabolic heat. The traits measured include hair density (HD), body condition score (BCS), body size (BS), ear length (EL), and ear area (EA). HD was measured as a subjective visual score ranging from 0 to 2, being the score 0 = hairless or limited hair cover, 1 = normal or moderate hair cover, and 2 = sow with greater than normal hair cover. BCS was evaluated using a sow caliper (BCS_cal) [109] and it was also evaluated as a visual BCS (BCS_vis) considering five categories: 1 = emaciated, 2 = thin, and 3 = ideal, 4 = fat, and 5 = overly fat. The animals were also score into three categories according to their BS as small, medium, or large. The descriptive statistics for the traits used in this study are presented in Table 4. More details about the data collection procedures are presented in Johnson et al. [10].

Table 4 Descriptive statistics for continuous heat stress indicators in Landrace x Large White crossbred pig population

The phenotypes were pre-corrected for fixed effects to be used as input for the GWAS analyses. The fixed effects for each trait were described by Freitas et al. [11]. In summary, the fixed effects fitted for skin surface temperature (T_ES, T_SS, T_RS, and T_TS) and RR were trait recorder; concatenation of week, day, and time of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For PS, the fixed effects were trait recorder; concatenation of week and day of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For T_Vall, the fixed effects were the concatenation of week and day of measurement; parity; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For T_V4days, the fixed effects were concatenation of week, day, and time of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For T_V8h, T_V12h, T_V16h, and T_V20h, the fixed effects were the concatenation of week and day of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For HD, the fixed effects fitted were trait recorder and parity. For BS, the fixed effects fitted were trait recorder, week of measurement, and parity. For BCS_cal and BCS_vis, the fixed effect fitted were trait recorder, week of measurement, parity, days in lactation, and concatenation of barn type and room. For EA and EL, the fixed effects fitted were trait recorder and picture quality. For traits that presented repeated measurements (i.e., the thermoregulatory traits: skin surface temperatures, vaginal temperatures, RR, and PS), we computed the average adjusted phenotype to be used as response variables for GWAS [110].

Genome-wide association studies

The GWAS was performed considering the WGS dataset using the leave-one-chromosome-out approach on mixed linear model association analysis (MLMA-LOCO) implemented in the GCTA software [111]. The model can be described as:

$${\bf{y}} = {\bf{1}}\mu + {\bf{X}}{\rm{b}} + {\bf{Zu}} + \varepsilon ,$$

,

Where $\mathbf{y}$ is the vector of pre-corrected phenotypes for each analyzed trait; ${\bf{1}}$ is a vector of ones; ${\mu }$ is the overall mean; $\text{b}$ is the fixed effect of each tested SNP; $\mathbf{X}$ is a vector containing the genotype for SNP; $\mathbf{u}$ is the random vector of polygenic effect with $\mathbf{u} \sim \text{N}(0,\mathbf{G}{{\sigma }}_{u}^{2})$, where $\mathbf{G}$ is the genomic-based relationship matrix (GRM) as previously described [111], ${{\sigma }}_{u}^{2}$ is the additive genetic variance; $\mathbf{Z}$ is the incidence matrix for $\mathbf{u}$; and $\varvec{\epsilon }$ is a random vector of residual effects with $\varvec{\epsilon } \sim \text{N}(0,\mathbf{I}{{\sigma }}_{\epsilon }^{2})$, where $\mathbf{I}$ is an identity matrix and ${{\sigma }}_{\epsilon }^{2}$ is the residual variance.

The MLMA-LOCO approach is based on fitting a different genomic relationship matrix (GRM) to account for population stratification for the association tests of each chromosome by excluding the markers located on the chromosome where the tested SNP are located when creating the GRM for that chromosome [112]. The inclusion of the marker being tested in the GRM can reduce the statistical test power due to the double fitting of the marker in the model: both as a fixed effect for association and as a random effect as part of the GRM [112,113,114]. This phenomenon was termed “proximal contamination” by Listgarten et al. [113], who demonstrated that a mixed linear model excluding the marker from GRM is the mathematically correct approach and provided an efficient algorithm for mixed linear model analysis.

Adjustment for multiple testing

To check if there was potential population stratification due to the population structure, the genomic inflation factor (λ) [115] was evaluated, along with its 95% confidence interval. Due to the multiple testing nature of the GWAS method employed in this study, there is a cumulative probability of false-positive associated SNPs [116, 117]. In this sense, a Bonferroni correction at α = 0.05 for the genome-wise significance level [118] was applied, by dividing α by the number of independent chromosome segments ($Me$) to account for the dependence among tests due to LD. The $Me$ was calculated as proposed by Goddard et al. [119]:

$$Me= \frac{2NeLk}{\text{l}\text{o}\text{g}\left(NeL\right)},$$

where $Ne$ is the effective population size, which was computed as 85 according to populational LD [120]; $L$ is the average length of chromosome expressed in Morgan; and $k$ is the number of chromosomes. A SNP was considered as statistically significant if its p-value was smaller than 0.05/$Me$, which was equal to 1.21 × 10^− 6.

Pleiotropic effects

The multiple-trait analysis proposed by Bolormaa et al. [18] allows to evaluate the pleiotropy in groups of more than two traits. This approach is based on combining the results from the single-trait association analyses to calculate a multi-trait statistic and it has been demonstrated that the multi-trait approach proposed by them has a lower FDR than a single trait analysis (at the same significance test p-value) [18, 121]. First, the phenotypes were Cholesky-decorrelated because the Bolormaa et al. [18] approach was originally proposed to perform meta-analysis, i.e., independent studies [19, 20]. Then, the traits were grouped in seven different categories: skin surface temperature traits (i.e., T_ES, T_SS, T_RS, and T_TS), vaginal temperatures (i.e., T_Vall, T_V4days, T_V8h, T_V12h, T_V16h, and T_V20h), animal temperatures (all skin surface and vaginal temperature traits), respiration traits (i.e., RR and PS), thermoregulatory indicators (i.e., all skin surface and vaginal temperature traits, RR and PS), anatomical traits (i.e., BCS_cal, BCS_vis, HD, BS, EA, and EL), and all traits to evaluate the presence of pleiotropic effects in each trait group. The null hypothesis tested states that the SNP does not affect any of the $n$ traits in the group based on the ${\chi }^{2}$ with $n$ degrees of freedom. The multi-trait statistics for each SNP was calculated as presented by [18]:

$${\rm{Multi - trait}}\,{\chi ^2} = t_i^\prime {V^{ - 1}}{t_i},$$

where ${\varvec{t}}_{i}$ is a $n$ x 1 vector of signed t-values for the ${i}^{th}$ SNP across the $n$ trait, ${\varvec{t}}_{i}^{{\prime }}$ is the transpose of ${\varvec{t}}_{i}$, and ${\varvec{V}}^{-1}$ is the inverse of the $n$ x $n$ correlation matrix over all estimated SNP effects. As this method presents a distribution associated to the ${\chi }^{2}$ statistics, each marker presents a statistic and a p-value to infer if the marker is significant or not for pleiotropic effects for the group of traits evaluated. A Bonferroni correction was also applied in the same way described for the GWAS, then a variant was considered with significant pleiotropic effect if its p-value was smaller than 1.21 × 10^− 6.

Gene annotation and functional analyses

Gene annotation was performed within a region of ± 100 kb around the significant SNPs using the GALLO package [122]. The annotated data for Sus scrofa from the Ensembl database (www.ensembl.org/Sus_scrofa/Info/Index) were used in the analyses. The web-based tool Database for Annotation, Visualization, and Integrated Discovery (DAVID) [123, 124] was used for conducting Gene Ontology (GO) and KEGG pathway enrichment (P < 0.05) analyses to identify biological processes, molecular functions, cellular components, and biological pathways associated with the positional candidate genes identified. DAVID is a web-based tool for gene-enrichment analysis that utilizes a knowledge base aggregating information from various heterogeneous and widely distributed public databases [124].