Longitudinal genomic analyses of automatically-recorded vaginal temperature in lactating sows under heat stress conditions based on random regression models

Abstract

Background

Automatic and continuous recording of vaginal temperature (T_V) using wearable sensors causes minimal disruptions to animal behavior and can generate data that enable the evaluation of temporal body temperature variation under heat stress (HS) conditions. However, the genetic basis of T_V in lactating sows from a longitudinal perspective is still unknown. The objectives of this study were to define statistical models and estimate genetic parameters for T_V in lactating sows using random regression models, and identify genomic regions and candidate genes associated with HS indicators derived from automatically-recorded T_V.

Results

Heritability estimates for T_V ranged from 0.14 to 0.20 over time (throughout the day and measurement period) and from 0.09 to 0.18 along environmental gradients (EG, − 3.5 to 2.2, which correspond to dew point values from 14.87 to 28.19 ˚C). Repeatability estimates of T_V over time and along EG ranged from 0.57 to 0.66 and from 0.54 to 0.77, respectively. T_V measured from 12h00 to 16h00 had moderately high estimates of heritability (0.20) and repeatability (0.64), indicating that this period might be the most suitable for recording T_V for genetic selection purposes. Significant genotype-by-environment interactions (GxE) were observed and the moderately high estimates of genetic correlations between pairs of extreme EG indicate potential re-ranking of selection candidates across EG. Two important genomic regions on chromosomes 10 (59.370–59.998 Mb) and16 (21.548–21.966 Mb) were identified. These regions harbor the genes CDC123, CAMK1d, SEC61A2, and NUDT5 that are associated with immunity, protein transport, and energy metabolism. Across the four time-periods, respectively 12, 13, 16, and 10 associated genomic regions across 14 chromosomes were identified for T_V. For the three EG classes, respectively 18, 15, and 14 associated genomic windows were identified for T_V, respectively. Each time-period and EG class had uniquely enriched genes with identified specific biological functions, including regulation of the nervous system, metabolism and hormone production.

Conclusions

T_V is a heritable trait with substantial additive genetic variation and represents a promising indicator trait to select pigs for improved heat tolerance. Moderate GxE for T_V exist, indicating potential re-ranking of selection candidates across EG. T_V is a highly polygenic trait regulated by a complex interplay of physiological, cellular and behavioral mechanisms.

Background

Heat stress (HS), which occurs when an individual cannot dissipate body heat adequately to maintain thermal equilibrium under hot conditions, is a common problem in livestock production [1]. It compromises animal welfare and causes significant economic losses in animal production, reproductive performance, and animal health [2,3,4]. Furthermore, intensive genetic selection for a limited number of performance traits, such as litter size and milk yield (through litter weight at weaning), has contributed to increased metabolic heat production in modern sows and may have led to their greater sensitivity to environmental conditions [5]. Global warming further aggravates these effects and highlights the need for selecting more climatic resilient animals. Over the past decades, numerous studies have investigated the underlying mechanisms of the response to HS in various species [6,7,8,9,10,11]. Dozens of genes and quantitative trait loci (QTL) associated with response to HS have been identified [12, 13].

Measurements of core body temperature have been widely used and are considered reliable indicators for HS in livestock [14]. Unlike other core body temperature indicators (e.g., rectal temperature), automatically-measured vaginal temperature (T_V) reduces disruptions in animal behavior, captures diurnal changes in body temperature, and allows continuous data to be recorded [15, 16]. A sensitive HS indicator that is not affected by disturbances from animal behaviors can play an important role in breeding programs, especially for lactating animals that are more heat-sensitive due to greater metabolic heat production compared to non-lactating sows [17]. However, little is known about the genetics of continuously-monitored T_V (i.e., longitudinal data) as a HS indicator in lactating sows.

Random regression models (RRM) have been used in breeding to investigate longitudinal data, as they can provide more precise results over time and across environmental gradients (EG) than other methods [18, 19]. Legendre orthogonal polynomials (LEG) are common functions that are used to model the fixed and random trajectories for longitudinal traits. However, B-spline functions (BS) can provide a similar fit to the data without having implausible estimates at the extremes of the trajectory curves that often result when fitting high-order polynomials [20]. However, fitting too many BS parameters can result in convergence problems and high computational requirements.

Gene expression data from broiler small intestine tissue and caenorhabditis elegans has been shown to exhibit a temporal dynamic profile during HS conditions [21, 22], since some genes have a greater role in specific regulatory functions and stages of HS response or recovery. This suggests that the effects of HS-related genes may vary with time or with environmental conditions. In addition, Oliveira et al. [23] showed that, when longitudinal traits are analyzed in dairy cattle (e.g., milk production and somatic cell score), differential sets of candidate genes are identified for different lactation stages.

Intra-vaginal sensors [24] enable automatic recording of Tv, which enables evaluation of the genetics of response to HS in lactating sows from a longitudinal perspective, and to identify genomic regions, candidate genes, biological processes, and metabolic pathways associated with this trait. Such analyses could contribute to understanding potential changes in the biological mechanisms that underlie functions involved in coping with HS and to optimize genetic progress for climatic resilience. Worldwide breeding programs are conducted on farms with considerably different environmental conditions and management choices [25], and can provide data to investigate the presence of significant genotype-by-environment interactions (GxE) for heat tolerance in lactating sows. Therefore, the main objectives of this study were to: (1) define statistical models and estimate genetic parameters for automatically-recorded Tv in lactating sows at different time points throughout the day and along EG; (2) assess potential GxE interactions for heat tolerance of lactating sows; and (3) identify genomic regions and biological processes that contribute to response to HS in maternal-line lactating sows from a longitudinal perspective.

Methods

Animals and datasets

All datasets used in this study were collected on a commercial farm in Maple Hill, NC, USA (34.70738° N, 77.73653° W) as previously described by Johnson et al. [24]. Sows were provided ad libitum access to feed and water. The T_V of 1381 sows from the studied population, which included 1645 lactating sows (parities 2 to 7; Landrace × Large White cross), was automatically measured every 10 min for 5 consecutive days between days 8 to 20 of lactation, as previously described [24]. In total, 932,681 T_V records were obtained. First parity sows were not included in this study because the farm where the experiment was conducted had only second and later parity sows (a common practice in some regions in North America). The ambient temperature (T_a) and relative humidity (RH) of each room were recorded every five min [24], from which dew points were calculated using the Magnus–Tetens equation. Dew point is a thermal index that takes both temperature and humidity into account, and was chosen as environmental indicator due to its ability to accurately represent heat stress conditions, as previously demonstrated [30]. Outliers with environmental and T_V records that deviated by more than 3.5 standard deviation (SD) from the mean were discarded. T_V records below 37 °C were also deleted.

In total, 1639 sows (including all sows with phenotypic records) were genotyped using the PorcineSNP50K Bead Chip (50,703 SNPs; Illumina, San Diego, CA, USA). For quality control (QC) of the genotype data, single nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) higher than 0.01, or that had an extreme difference between observed and expected heterozygous frequencies less than 0.15, and SNPs and samples with a call rate greater than 0.90 were kept for further analyses. The QC was implemented using the BLUPF90 + family software [26] and 49,547 SNPs and 1639 animals remained for further analyses. A schematic representation of all analyses performed in the study is shown in Additional file 1: Fig. S1.

Statistical models to estimate genetic parameters across time

Twelve single-trait RRM models were evaluated, with random regressions fitted as time-derived covariates using LEG (linear, quadratic, and cubic) and BS (linear, quadratic, and cubic with 5, 6, or 7 knots). The general formulation of the RRM used for the estimation of variance components and genome-wide association study (GWAS) is:

$${y}_{ijklmn}={Par}_{k}+{Loc}_{l}+{Date}_{m}+\sum_{r=1}^{R}{\gamma }_{ir}{\phi }_{r}\left({t}_{ij}\right)$$

$$+\sum_{r=1}^{R}{\alpha }_{ir}{\phi }_{r}\left({t}_{ij}\right)+\sum_{r=1}^{R}{\rho }_{ir}{\phi }_{r}\left({t}_{ij}\right)+{e}_{ijklmn},$$

where ${y}_{ijklmn}$ is the $n$th T_V observation, recorded at the $j$th time point of animal $i$ from the $k$th parity, at location $l$, and date of measurement $m$; ${Par}_{k}$ is the fixed effect of the $k$th parity of animal $i$; ${Loc}_{l}$ is the fixed effect of the $l$th location (e.g., barn and room within barn) of animal $i$; ${Date}_{m}$ is the fixed effect of the $m$th measurement date of animal $i$; ${\gamma }_{ir}$ is the $r$th fixed regression coefficient specific for the time trajectory; ${\alpha }_{ir}$ and ${\rho }_{ir}$ are the $r$th random regression coefficients that describe the trajectory of the additive genetic effects and permanent environmental effects of animal $i$, which were assumed to follow ${\varvec{\upalpha}}\sim {\text{N}}\left(0,\mathbf{G}{\upsigma }_{{\text{a}}}^{2}\right)$ and ${\varvec{\uprho}}\sim {\text{N}}(0,\mathbf{I}{\upsigma }_{{\text{pe}}}^{2})$, respectively, where $\mathbf{G}$ is the genomic relationship matrix built according to VanRaden’s Method 1 [27]; ${\phi }_{r}({t}_{ij})$ is the covariate of the regression function with time using LEG or BS; ${t}_{ij}$ is the $j$th time point at which T_V was measured throughout the day (24 h) for animal $i$; and $R$ is the order of the LEG or BS for the fixed regression effects, genetic additive random effects, and permanent environmental effects. The same polynomial order was fitted for the additive genetic and permanent environmental effects, as previously suggested [28]. The random residual effect (${e}_{ijklmnr}$) was assumed to be distributed homogeneous or heterogeneous throughout the day. For the latter, residual variance (RV) was modelled considering six periods of 4 h (00h00–04h00, 04h00–08h00, 08h00–12h00, 12h00–16h00, 16h00–20h00, 20h00–24h00)].

The LEG were obtained as proposed by Kirkpatrick et al. [29] and the first three polynomial orders were evaluated. B-spline functions [20] with equally-spaced knots (n = 5, 6, or 7) and degrees (linear—L, quadratic—Q, cubic—C) for the additive genetic and permanent environmental effects were evaluated. The RRM models using BS are referred to as ‘BSX K_a K_pe’, where X = L, Q, or C, that is, the degree of each polynomial segment, and K_a and K_pe are the numbers of random regression coefficients for the additive genetic and permanent environmental effects, respectively. The RRM models using LEG are referred to as ‘LEG X’, where X is the Legendre orthogonal polynomial order for the additive genetic and permanent environmental effects. Thus, model BSQ88 consists of a quadratic B-spline function with eight random regression coefficients for the additive genetic and permanent environmental effects, and model LEG4 fits additive genetic and permanent environmental regressions based on fourth-order Legendre orthogonal polynomials. The (co)variances components and genetic parameters for T_V over 24 h were estimated based on the best model (lowest Bayesian information criterion—BIC and highest accuracy of genomic predictions) using the AIREML algorithm implemented in the BLUPF90 + program [26].

Statistical models to estimate genetic parameters across environmental gradients

A reaction norm model (RNM) was implemented using the BLUPF90 + program [26] to evaluate GxE interactions, with EG derived by standardizing the dew point as: standardized dew point = (actual dew point – mean dew point)/standard deviation of dew point. The mean and standard deviation of dew point was 23.052 °C and 2.335 °C, respectively. The statistical model used can be described as:

$${y}_{ijklmno}=\alpha 1+{Par}_{k}+{Loc}_{l}+{Hour}_{m}+\beta {\varphi }_{1k}$$

$$+{a}_{0i}1+{a}_{1i}{\varphi }_{1k}+{p}_{0i}1+{p}_{1i}{\varphi }_{1k}+{e}_{ijklmno},$$

where ${y}_{ijklmno}$ is the $o$th observation for T_V for animal $i$; $\alpha$ is the intercept; ${Par}_{k}$ is the fixed effect of the $k$th parity of animal $i$; ${Loc}_{l}$ is the fixed effect of the $l$th location (e.g., barn and room within barn) of animal $i$; ${Hour}_{m}$ is the fixed effect of the $m$th measurement hour for animal $i$. $\beta$ is the fixed regression coefficient on EG, ${\varphi }_{1k}$ is the EG vector value $k$, ${a}_{0l}$ and ${a}_{1l}$ are the random regression coefficients for the intercept and slope of the additive genetic effect of individual $i$, ${p}_{0i}$ and ${p}_{1i}$ are the random regression coefficients for the intercept and slope of the random permanent environmental effect of individual $i$, and ${{\text{e}}}_{{\text{ijklmno}}}$ the residual effect, which was assumed distributed homogeneous or heterogeneous throughout the HS period, as described above. The following assumptions were made for the additive genetic effects: $\left[\begin{array}{c}{{\text{a}}}_{0}\\ {{\text{a}}}_{1}\end{array}\right]\sim {\text{N}}\left(0, \mathbf{G}\otimes {\mathbf{K}}_{\mathbf{a}\mathbf{b}}\right)$ and $\left[\begin{array}{c}{{\text{p}}}_{0}\\ {{\text{p}}}_{1}\end{array}\right]\sim {\text{N}}\left(0, \mathbf{I}\otimes {\mathbf{K}}_{\mathbf{c}\mathbf{d}}\right)$, and ${\mathbf{K}}_{\mathbf{a}\mathbf{b}}$ and ${\mathbf{K}}_{\mathbf{c}\mathbf{d}}$ are both 2 × 2 (co)variance matrices for the intercept and slope effects, with ${\mathbf{K}}_{\mathbf{a}\mathbf{b}}=\left[\begin{array}{cc}{\upsigma }_{0}^{2}& {\upsigma }_{01}\\ {\upsigma }_{10}& {\upsigma }_{1}^{2}\end{array}\right]$, where ${\upsigma }_{0}^{2}$ is the additive genetic variance for the intercept term, ${\upsigma }_{1}^{2}$ is the additive genetic variance for the slope term, ${\upsigma }_{10}$ and ${\upsigma }_{01}$ are the covariances between the two aforementioned effects, and ${\mathbf{K}}_{\mathbf{c}\mathbf{d}}=\left[\begin{array}{cc}{\upsigma }_{{\text{p}}0}^{2}& {\upsigma }_{{\text{p}}0{\text{p}}1}\\ {\upsigma }_{{\text{p}}1{\text{p}}0}& {\upsigma }_{{\text{p}}1}^{2}\end{array}\right]$, where ${\upsigma }_{{\text{p}}0}^{2}$ is the permanent environmental variance for the intercept term, ${\upsigma }_{{\text{p}}1}^{2}$ is the permanent environmental variance for the slope term, ${\upsigma }_{{\text{p}}1{\text{p}}0}$ and ${\upsigma }_{{\text{p}}0{\text{p}}1}$ are the covariances between the two aforementioned effects.

Genetic and environmental (co)variance matrices for each time point were defined and calculated as in Oliveira et al. [19]. A one-tailed t-test was performed to determine if there was significant GxE interactions, indicated by the variance of the RNM slope significantly differing from zero. After obtaining (co)variance components, estimates of heritability and repeatability for T_V over time were calculated based on the following equations:

$${h}_{i}^{2}=\frac{{\widehat{\sigma }}_{{a}_{i}}^{2}}{{\widehat{\sigma }}_{{a}_{i}}^{2}+{\widehat{\sigma }}_{{pe}_{i}}^{2}+{\widehat{\sigma }}_{{e}_{i}}^{2}},$$

$${r}_{i}=\frac{{\widehat{\sigma }}_{{a}_{i}}^{2}+{\widehat{\sigma }}_{{pe}_{i}}^{2}}{{\widehat{\sigma }}_{{a}_{i}}^{2}+{\widehat{\sigma }}_{{pe}_{i}}^{2}+{\widehat{\sigma }}_{{e}_{i}}^{2}},$$

where ${h}_{i}^{2}$, ${r}_{i}$, ${\widehat{\sigma }}_{{a}_{i}}^{2}$, ${\widehat{\sigma }}_{{pe}_{i}}^{2}$, ${\widehat{\sigma }}_{{e}_{i}}^{2}$ are the estimates of heritability, repeatability, additive genetic variance (VAG), permanent environmental variance (VPE), and RV for the $i$th time point. The estimate of heritability for EG_i was calculated as follows [31], ${{\text{h}}}_{{\text{i}}}^{2}=\frac{{\widehat{\upsigma }}_{{{\text{u}}}_{{\text{i}}}}^{2}}{\sum {\widehat{\upsigma }}_{{{\text{n}}}_{{\text{i}}}}^{2}+{\widehat{\upsigma }}_{{\text{e}}}^{2}}$, where ${\widehat{\upsigma }}_{{{\text{u}}}_{{\text{i}}}}^{2}$ is the estimate of the additive genetic variance, which was computed as ${\widehat{\upsigma }}_{{{\text{u}}}_{{\text{i}}}}^{2}={\widehat{\upsigma }}_{{{\text{a}}}_{0}}^{2}+2{\widehat{\upsigma }}_{{{\text{a}}}_{0}{{\text{a}}}_{1}}{\widehat{\uptheta }}_{{\text{i}}}+{\widehat{\upsigma }}_{{{\text{a}}}_{1}}^{2}{({\widehat{\uptheta }}_{{\text{i}}})}^{2}$, and the denominator is the estimate of the phenotypic variance, with $\sum {\widehat{\upsigma }}_{{{\text{n}}}_{{\text{i}}}}^{2}=\sum {\widehat{\upsigma }}_{{{\text{n}}}_{0}}^{2}+2{\widehat{\upsigma }}_{{{\text{n}}}_{0}{{\text{n}}}_{1}}{\widehat{\uptheta }}_{{\text{i}}}+{\widehat{\upsigma }}_{{{\text{n}}}_{1}}^{2}{({\widehat{\uptheta }}_{{\text{i}}})}^{2}$, where ${\text{n}}$ refers to the random effects fitted for T_V. For the model fitted with heterogenous residual variances, the component ${\widehat{\upsigma }}_{{\text{e}}}^{2}$ was calculated as ${\widehat{\upsigma }}_{{{\text{e}}}_{{\text{i}}}}^{2}={\text{exp}}({{\text{d}}}_{0}+{{\text{d}}}_{1}{\widehat{\uptheta }}_{{\text{i}}})$. The estimate of the genetic correlation between EG_i and ${{\text{i}}}^{\mathrm{^{\prime}}}({{\text{r}}}_{{{\text{ii}}}^{\mathrm{^{\prime}}}})$ was calculated as: ${{\text{r}}}_{{{\text{ii}}}^{\mathrm{^{\prime}}}}=\frac{{\widehat{\upsigma }}_{{{\text{u}}}_{{{\text{ii}}}^{\mathrm{^{\prime}}}}}}{\sqrt{{\widehat{\upsigma }}_{{{\text{u}}}_{{\text{i}}}}^{2}{\widehat{\upsigma }}_{{{\text{u}}}_{{{\text{i}}}^{\mathrm{^{\prime}}}}}^{2}}}$, where ${\widehat{\upsigma }}_{{{\text{u}}}_{{{\text{ii}}}^{\mathrm{^{\prime}}}}}$ is the estimate of the covariance of additive genetic effects between EG_i and ${{\text{i}}}^{\mathrm{^{\prime}}}$, which was computed as ${\widehat{\upsigma }}_{{{\text{u}}}_{{{\text{ii}}}^{\mathrm{^{\prime}}}}}={\widehat{\upsigma }}_{{{\text{a}}}_{0}}^{2}+{\widehat{\upsigma }}_{{{\text{a}}}_{0}{{\text{a}}}_{1}}{\widehat{\uptheta }}_{{\text{i}}}+{\widehat{\upsigma }}_{{{\text{a}}}_{0}{{\text{a}}}_{1}}{\widehat{\uptheta }}_{{{\text{i}}}^{\mathrm{^{\prime}}}}+{\widehat{\upsigma }}_{{{\text{a}}}_{1}}^{2}{\widehat{\uptheta }}_{{\text{i}}}{\widehat{\uptheta }}_{{{\text{i}}}^{\mathrm{^{\prime}}}}$.

Comparison of models

The BIC [32] was used to choose the best models, with models with a lower BIC value providing a better fit of the data. The BIC for a model with k estimated parameters and n observations was calculated as:

$${\text{BIC}}=-2{\text{logL}}+{\text{ln}}\left({\text{n}}\right){\text{k}},$$

where $-2{\text{logL}}$ is the restricted maximum log likelihood value of the model. In addition, fivefold cross-validation was used to compare prediction accuracies of the RRM using LEG and BS functions (LEG4 and BSQ88) [33]. For this purpose, the population was randomly allocated to five groups. For a given fold, one group was in turn assigned with missing phenotypic values and used as a validation dataset, and the other four groups were used as the training dataset. Genomic estimated breeding values (GEBV) for the full dataset were calculated for all animals considering the phenotypes from all sows in all groups. The Pearson correlations between GEBV for each time point for sows from the validation dataset and the corresponding GEBV for sows from the full dataset were calculated. The obtained average prediction accuracies across time points were used to select the best model for further analyses.

Genome-wide association studies

Genome-wide association studies were implemented with a sliding window of 10 consecutive SNPs to map genomic regions that are associated with T_V, using the postGSf90 software [26]. The same models used to estimate the variance components were implemented for the GWAS. Four HS time-periods throughout the day were used to investigate the genomic regions that are associated with T_V at each of these time-periods: (a) from 23h00h to 06h30, which represents the period when T_V starts to decrease together with T_a and with relatively more thermally comfortable conditions; (b) from 06h30 to 09h30, which represents the period when T_V is maintained at thermoneutral levels because the environment is not too hot for the sows; (c) from 09h30 to 18h30, which represents the period when T_V increases with increasing T_a; and (d) from 18h30 to 23h00,which represents the period T_V is maintained at a relatively high level because the environment is usually above thermoneutral conditions. Furthermore, three classes of EG were defined to investigate the genomic regions that are associated with T_V under different EG conditions. The following three EG classes were defined based on standardized Dew Point units: no to mild HS: − 3.5 to 1.5; moderate HS: − 1.5 to 0.5; and severe HS: 0.5 to 2.5. The GEBV for Tv for an animal was obtained by summing its GEBV at different time points or for different EG classes. For instance, the GEBV for animal ${\text{i}}$ for the four HS time-periods, ${{\text{GEBV}}1}_{{\text{i}}}$, ${{\text{GEBV}}2}_{{\text{i}}}$, ${{\text{GEBV}}3}_{{\text{i}}}$, ${{\text{GEBV}}4}_{{\text{i}}},$ was computed by summing the GEBV at each time point (10-min interval) from 23h00 to 06h30, 06h30 to 09h30, 09h30 to 18h30, 18h30 to 23h00, respectively, as: ${{\text{GEBV}}1}_{{\text{i}}}={{\text{GEBV}}}_{{\text{i}}2305}+{{\text{GEBV}}}_{{\text{i}}2315}+\dots +{{\text{GEBV}}}_{{\text{i}}0635},$

$${{\text{GEBV}}2}_{{\text{i}}}={{\text{GEBV}}}_{{\text{i}}0645}+{{\text{GEBV}}}_{{\text{i}}0655}+\dots +{{\text{GEBV}}}_{{\text{i}}0935},$$

$${{\text{GEBV}}3}_{{\text{i}}}={{\text{GEBV}}}_{{\text{i}}0945}+{{\text{GEBV}}}_{{\text{i}}0955}+\dots +{{\text{GEBV}}}_{{\text{i}}1835},$$

$${{\text{GEBV}}4}_{{\text{i}}}={{\text{GEBV}}}_{{\text{i}}1845}+{{\text{GEBV}}}_{{\text{i}}1855}+\dots +{{\text{GEBV}}}_{{\text{i}}2255}.$$

The SNP effects for each time-period were calculated for T_V using the postGSf90 program [26] as:

$${\widehat {\bf{u}}_{\bf{m}}} = {\bf{D}}{{\bf{Z}}^\prime }{[{\bf{ZD}}{{\bf{Z}}^\prime }]^{ - 1}}{\bf{GEB}}{{\bf{V}}_m},$$

where ${\widehat{\mathbf{u}}}_{\mathbf{m}}$ is the vector of the estimated SNP effects for $m$th HS time-period; $\mathbf{D}$ is a diagonal matrix of the weights for the variances of SNP effects ($\mathbf{D}=\mathbf{I}$ for GBLUP), $\mathbf{Z}$ is a matrix relating the genotype (gene content) at each SNP to the individual, and ${\mathbf{G}\mathbf{E}\mathbf{B}\mathbf{V}}_{{\varvec{m}}}$ is the vector of the estimated GEBV for $m$th HS time-period, which includes the GEBV for all animals. The genetic variance explained by SNP $k$ (${\widehat{\sigma }}_{u,k}^{2}$) was estimated as:

$${\widehat{\sigma }}_{u,k}^{2}={\widehat{u}}_{k}^{2}2{p}_{k}\left(1-{p}_{k}\right),$$

where ${p}_{k}$ is the observed allele frequency for the first allele of SNP $k$, and ${\widehat{u}}_{k}^{2}$ is the square of the estimated effect of SNP $k$. The percentage of the total additive genetic variance explained by the $i$th 10-SNP moving window was calculated as:

$$\frac{{{\text{Var}}(\widehat {{a_i}})}}{{\widehat {\sigma _a^2}}} \times 100\% = \frac{{{\text{Var}}\left( {\sum _{{\text{j}} = 1}^{10}{{\bf{z}}_{\text{j}}}{{\widehat {\text{u}}}_{\text{j}}}} \right)}}{{\widehat {\sigma _a^2}}} \times 100\% ,$$

where $\widehat{{a}_{i}}$ is the genetic variance explained by the ${\text{i}}$th window that consists of 10 consecutive SNPs; $\widehat{{\sigma }_{a}^{2}}$ is the total additive genetic variance; ${\mathbf{z}}_{\mathbf{j}}$ is the gene content vector of the ${\text{j}}$th SNP for all individuals, and ${\widehat{{\text{u}}}}_{{\text{j}}}$ is marker effect of the ${\text{j}}$th SNP within the ${\text{i}}$th region. The same method was used to obtain the genomic window variance for the EG classes. Genomic windows that explained more than 0.25% of the genetic variance were considered as relevant genomic regions [34]. To avoid double-counting of SNP effects, only non-overlapping genomic windows were kept for further analyses, as suggested by Fragomeni et al. [35].

Functional genomic analyses

The genes within the relevant genomic regions were annotated based on the latest pig reference genome Sscrofa 11.1 assembly (http://useast.ensembl.org/Sus_scrofa/Info/Index). Gene Ontology (GO) [36] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [37] enrichment analyses for candidate genes were carried out via the “clusterProfiler” package of R with a cutoff P-value < 0.05 and a false discovery rate (FDR) < 0.20 [38].

Results

The number of records and the average Tv for each time point and each EG are presented in Fig. 1. Each time point included more than 6000 records. The T_V had a circadian rhythm and remained relatively constant from 06h30 to 09h00 and from 18h30 to 23h00. The number of records was lower for the extreme EG. The average T_V showed a distinct pattern along the EG scale, which is based on standardized dew point values. Initially, T_V increased significantly at low dew point values. This increase slowed down as the EG value approached 1.7. Beyond this point, the T_V’s increase became sharp once again (Fig. 1b).

Fig. 1

Number of records and average vaginal temperature throughout the day and along environmental gradients. a Number of records (green bars) and average vaginal temperature (T_V) (red dots) for time points (every 10 min) throughout the day and b number of records (green bars) and average T_V (red dots) along environmental gradients based on dew point values

Statistical model comparisons

Table 1 summarizes the results of the model comparisons based on the logarithm of the likelihood (Log L) and BIC. The number of parameters ranged from 21 to 111 across the models evaluated. The largest number of knots that converged when fitting quadratic BS was 7. For the BS models, BIC decreased as the number of knots and order of the function increased. The BSQ88 model outperformed (lowest BIC) all other models (Table 1). The patterns of the accuracy of genomic prediction over time are shown in Additional file 2: Fig. S2 and the accuracy for T_V were 0.72 and 0.13 with models LEG4 and BSQ88, respectively (see Additional file 2: Fig. S2). These results indicate that BSQ88 might have overfitted the data. Thus, LEG4 was used for further analyses.

Table 1 Number of model parameters (N), knot position, log likelihood (Log L), Bayesian information criterion (BIC) for random regression models fitted to longitudinal Tv data using B-splines functions or Legendre orthogonal polynomials

Variance components and genetic parameters

Figure 2 shows estimates of variance components, heritability, and repeatability for T_V throughout the day and along the EG when fitting heterogeneous and homogeneous RV for lactating sows under HS conditions. The model with heterogeneous RV had lower BIC values than those considering homogeneous RV, which implied better goodness of fit. Estimates of VAG, VPE, heritability, and repeatability for T_V obtained by the models that fit heterogeneous versus homogeneous RV had similar trends and values. Estimates of heritability (from 0.14 to 0.20) and repeatability (from 0.57 to 0.66) for T_V fluctuated over time, with slightly different trends. Heritability estimates for T_V increased from 0.09 along the EG scale, with the highest estimate (0.18) at a standardized dew point of 1.2 (25℃), after which it stabilized. Repeatability estimates of T_V ranged from 0.54 to 0.77 along the EG.

Fig. 2

(a left column) Additive genetic variance, permanent environmental variance, residual variance, heritability, and repeatability along the time trajectory throughout the day and (b right column) continuous environmental gradient for vaginal temperature of lactating sows based on random regression model analyses. Red dot points (●) represent models considering heterogeneous residual variance and green triangles (▲) represent models considering homogeneous residual variance

The trajectories of the estimates of genetic correlations for T_V over time and along the EG are shown in Fig. 3a and b, respectively. Over time, the genetic correlation estimates ranged from 0.77 to 1.00, and were lowest for Tv between 09h00–12h00 and 01h00–04h00, and between 04h00–08h00 and 18h00–24h00. For the EG, genetic correlation estimates were positive and moderate to high, ranging from 0.45 to 1.00 (mean correlation: 0.75 ± 0.31) with the lowest estimate observed between the most divergent EG values. The variance of the slope for Tv (0.017) in the RNM model was significantly different from zero using a one-tailed t-test (P < 0.05). The estimate of the genetic correlation between the intercept and the slope for T_V was 0.33.

Fig. 3

Estimates of genetic correlations (Cor) for vaginal temperature in lactating sows under heat stress between a pairs of time points (each 10 min) throughout the day (range: 0.76–1) and b pairs of environmental gradients (standardized dew point) (range: 0.44–1)

Association analyses

As mentioned above, LEG4 was used for the association analyses. In total, 12, 13, 16, and 10 relevant 10-SNP windows across more than 10 Sus scrofa chromosomes (SSC) were identified for T_V across the four time-periods 1, 2, 3, and 4, respectively (Table 2). The total additive genetic variance explained by all relevant genomic windows for the four HS time-periods was 16.9%. Based on the candidate genes identified, four and one KEGG pathways were significantly enriched for the HS time-periods 3 and 4, respectively (see Additional file 3: Tables S1 and S2). For the GO analyses, 12, 21, and 13 terms were significantly enriched for HS time-periods 1, 3 and 4, respectively (see Additional file 3: Table S3). The results from the KEGG and GO analyses indicate that development (embryonic skeletal system development), neural system (glutamatergic synapse), and cardiac disease related pathways (hypertrophic cardiomyopathy, dilated cardiomyopathy) were the most enriched pathways based on the candidate genes identified.

Table 2 Relevant genomic regions identified for the four heat stress (HS) time-periods

Eighteen, 15, and 14 10-SNP windows were identified for T_V for the three EG classes (mild HS: − 3.5 to 1.5; moderate HS: − 1.5 to 0.5; and severe HS: 0.5 to 2.5, respectively) (Table 3). The total additive genetic variance explained by all the relevant genomic windows for the three EG classes was 16.3%. No significant KEGG pathway and, respectively, 30, 23, and 58 GO terms (see Additional file 3: Table S4) were significantly enriched based on the candidate genes (see Additional file 3: Table S5) identified for the three EG classes. Enriched GO terms were mainly associated with development (e.g., skeletal system development), energy metabolism (e.g., fatty acid biosynthesis processes), and hormone regulation (e.g., regulation of secretion).

Table 3 Relevant genomic regions identified for three environmental gradient (EG) classes

Two important genomic regions were identified to be associated with Tv across all HS time-periods and all EG classes, on SSC10 (59.370–59.998 Mb) and SSC16 (21.548–21.966 Mb). Candidate genes within these regions are involved in immunity (CDC123 and CAMK1d), protein transport (SEC61A2), and energy metabolism (NUDT5) functions. The overlapping and unique candidate genes regulating T_V across the four HS time-periods and the three EG classes are presented in a Venn diagram in Fig. 4. HS time-stage 3 had the most uniquely enriched genes (18), while HS time-stage 2 had the least number of uniquely enriched genes (7). One, 16, 21, and 7 significant KEGG pathways were detected based on the uniquely enriched genes for the four respective HS time-periods (Table 4). Furthermore, EG classes 1, 2, and 3 had 20, 1, and 9 uniquely enriched genes, respectively. Twenty-nine genes were shared by all three EG classes. Three, five, and three significant KEGG pathways were detected based on the uniquely enriched genes for EG classes 1, 2, and 3, respectively (Table 5).

Fig. 4

Venn diagram of the count of candidate genes identified based on a heat stress (HS) time-periods and b environmental gradient class. The four potential HS time-periods throughout the day were defined as (1) from 23h00 to 06h30, representing when vaginal temperature (T_V) starts to decrease with decreased ambient temperature (T_a) and with relatively more comfortable environmental conditions; (2) from 06h30 to 09h30, representing T_V maintained at a relatively low level when environment is not too hot for the sows; (3) from 09h30 to 18h30, representing the time in which T_V starts to increase with increasing T_a; (4) from 18h30 to 23h00, representing T_V maintained at a relatively high level when T_a is usually above thermoneutral conditions. Each two standardized Dew point units was considered as an environmental gradient class (no to mild HS: − 3.5 to 1.5; moderate HS: − 1.5 to 0.5; severe HS: 0.5–2.5)

Table 4 Significantly enriched (P < 0.05) KEGG pathways using uniquely enriched positional candidate genes for Tv for the four heat stress (HS) time-periods

Table 5 Significantly enriched (P < 0.05) KEGG pathways using uniquely enriched candidate genes identified for the three environmental gradient classes

Discussion

Climate change threatens worldwide livestock production and causes substantial economic and animal welfare issues. Understanding the genetic basis of an animal’s HS response, including identification of novel indicator traits such as T_V, is paramount for designing effective strategies for breeding more resilient and productive pigs. In this study, genetic parameters and GWAS analyses were conducted for T_V that was automatically measured every 10 min in Landrace × Large White lactating sows under HS conditions. The datasets of automatically-measured T_V and within-barn environmental measurements provided an opportunity to explore the genetic basis of longitudinal variability in Tv during lactation.

Comparison of models

The order of the BS had a greater impact on estimates of variance components than the number of knots. The genetic parameter estimates followed similar trends for all cubic BS models. Lower BIC values for models with a larger number of parameters have been reported previously (e.g., [20, 39, 40]). Although LEG models had larger BIC values than BS models, the computing times of the LEG models were much shorter. The lower accuracy of genomic prediction obtained with the RRM model using BS indicated that BS models may overfit the data, and based on the cross-validation analyses the LEG models were considered to be better.

Estimates of genetic parameters and variance components

To our knowledge, this study presents the first estimates of heritability and variance components for T_V over time and along EG under HS conditions in lactating sows. The use of longitudinal automatically-measured T_V is expected to yield more accurate estimates. Estimates of variance components along the EG scale were generally larger than those across time periods, while the heritability estimates across periods were higher. The moderate level of the average heritability estimate for T_V obtained with the RRM with LEG (0.18) indicates that T_V can be included in selection indexes as a heat tolerance indicator. Another measure of internal temperature, rectal temperature has been reported to have similar heritability estimates in lactating Holstein cows (0.17 [41] and 0.15 to 0.31 [42]). In those studies, the temperature was collected only once a day, which did not allow changes in heritability estimates and genetic variances throughout the day or with climatic conditions to be evaluated. Nevertheless, further studies are needed to evaluate the genetic correlations between T_V and other economically important traits (e.g., weaning weight) in pigs. If T_V can only be measured during a short time-window, the time-period from 1200 to 1600 h is recommended because it had the highest average heritability (0.20) and repeatability (0.64).

Trends in heritability estimates across time and across EG were similar to trends in average Tv (Fig. 1). A substantial increase in average T_V was observed after the 1.7 value for EG (dew point: 27 Celsius degree), which suggests that animals may experience severe HS above this dew point and may not be able to maintain body temperature homeostasis, thus needing cooling devices to maintain animal welfare. Trends in heritability estimates along the EG indicated that genetics has an increasing role in controlling T_V as the dew point increases. Thus, greater selection response to T_V is expected in a hotter environment due to the greater genetic variation, likely because more physiological and behavioral processes are involved in body temperature regulation under severe HS conditions [43]. Previous studies reported that heritability estimates for rectal temperature under HS conditions ranged from 0.11 to 0.36 in many species, which is in line with our results [41, 44, 45]. Also, the high VPE and repeatability estimates for Tv at the extremes of the EG scale indicate that relatively constant and accurate records could be obtained under extreme or mild HS conditions (Fig. 2b).

The variance for the slope for T_V from the RNM models was significantly different from zero, which indicates that GxE interactions exist when lactating sows are kept under HS conditions. The estimate of the genetic correlation being lowest (0.47) between extreme EG also supports the existence of GxE interactions and implies that animals will be ranked differently along the EG scale (Fig. 3b). Heat resilient individuals that can maintain performance levels under HS contribute to a greater profitability under HS conditions [46]. Cooling devices (i.e., sprinkler systems and evaporative cooling system) and changes in animal behavior (i.e., reduced feed intake, reduced milk production, and laying down for longer times) can also mitigate the negative effects of HS and improve the animals’ actual performance by reducing air temperature [2] or metabolic heat production [47]. Future research on the genetics of longitudinal Tv should consider larger datasets with greater environmental variability and different populations.

GWAS

Understanding the genetic mechanisms that underlie changes in T_V of lactating sows under HS is important to alleviate adverse effects from HS and to select heat stress resilient animals. This is the first study to report GWAS results for continuously-recorded T_V in lactating sows under HS conditions, mainly because obtaining relevant data is difficult. In this study, we performed a GWAS to investigate the dynamic regulatory functions of HS-linked genomic regions, genes, and QTL. Automatically-measured body temperature data and environmental records capture more variation and allow a more accurate detection of genomic regions as compared to traditional reaction norm models based on performance traits and data from public weather stations.

Multiple genomic regions on different chromosomes but with relatively small contributions to the total additive genetic variation of T_V were identified across the different measurement periods, indicating that T_V is a highly polygenic trait (Tables 2 and 3). Previous reports have indicated that rectal temperature, which is strongly genetically correlated with T_V in dairy cattle, sows, and sheep (r ≥ 0.8), is also a highly polygenic trait [15, 48,49,50]. Our results showed that some effects of the identified genomic regions were consistent throughout the entire HS period, while most genes or genomic regions associated with T_V regulation are involved in different HS time-periods and EG classes. This suggests that genetic regulation of T_V is complex and influenced by both common and unique genomic regions. Regulation of T_V can also be affected by farm management (i.e., cooling device), environmental conditions (i.e., good ventilation), and animal physiological state [16, 43].

Six genes were identified for T_V under HS conditions: CAMK1D, RNU6ATAC39P, CDC123, NUDT5, SEC61A2, and RANBP3L, many of which have previously been reported to be related to HS. The genes CAMK1D and CDC123 are mainly related to immunity. Heat stress can result in intestinal and systematic inflammation and thus in suppressing the innate immune function and increasing the animals’ susceptibility to diseases [51]. Many studies have also demonstrated that HS has detrimental effects on immunity and negatively affects health in cattle [52, 53], pigs [54, 55], sheep [56, 57], goats [58, 59], and poultry [60, 61]. Another study reported that SEC61A2 was expressed during HS [62].

GWAS enrichment over time

The results for KEGG pathway enrichment were similar for the first two HS time-periods, and the enriched pathways were mainly related to the nervous system (i.e., glutamatergic synapse) and platelet activation. The pathways that were enriched for HS time-period 3 were mainly related to signaling (i.e., oxytocin signaling pathway, mTOR signaling pathway, and hippo signaling pathway), cardiac disease (i.e., hypertrophic cardiomyopathy, dilated cardiomyopathy), and hormone synthesis and regulation (i.e., aldosterone synthesis and secretion). Previous studies revealed that the cardiovascular system plays a key role in human thermoregulation because heat dissipation is accompanied by cardiovascular adjustments, including increased cardiac output through a combination of elevations in heart rate and cardiac contractility and elevations in sympathetic activity for reduced blood flow and blood volume [63,64,65]. The cardiovascular system can be negatively affected by HS [66]. Occurrences of cardiovascular and respiratory mortality have been shown to increase significantly under HS conditions, especially in tropical areas [67, 68]. Under HS conditions, animals can experience severe dehydration, which would activate the renin–angiotensin–aldosterone system to maintain fluid and electrolyte balance [69]. Previous studies in cattle have shown that aldosterone concentration declines substantially under HS [70]. In our study, the pathways that were enriched for the HS time-period 4 were associated with various functions, including signaling (i.e., calcium signaling pathway and oxytocin signaling pathway), protein export, and platelet activation.

Many GO terms that were enriched for the four HS time-periods were related to development (i.e., embryonic skeletal system morphogenesis, embryonic skeletal system development), metabolic process (i.e., glucocorticoid metabolic process, carbohydrate metabolic process), nervous system (i.e., dendrite development, regulation of neuron projection development, and dendritic spine morphogenesis), and immunity (i.e., lymphocyte mediated immunity, natural killer cell mediated immunity, and leukocyte mediated immunity). This suggests that a sow’s response to HS is a complex process that involves a series of nervous, physiological, cellular, and molecular processes that can result in greater stress resilience, which is in line with previous studies [71, 72].

To identify which biological pathways play key and unique roles in response to HS and to better understand their functions, KEGG enrichment analyses of uniquely enriched genes were performed for each HS time-period. Although many similarities were found between the four HS time-periods, specific biologic functions were identified for each period. For HS time-period 1, only one significant pathway was identified, i.e. mitophagy-animal, which is known to increase resistance to diverse stressors and improve longevity [73]. Several studies have reported that mitophagy is impaired upon exposure to HS, oxidative stress, or other stressors [74,75,76].

Several pathways related to synaptic transmission were identified for HS time-period 2, including glutamatergic synapse, serotonergic synapse, dopaminergic synapse, GABAergic synapse, cholinergic synapse, and retrograde endocannabinoid signaling. Previous studies have shown that thermal stress is associated with suppressed overall synaptic transmission (especially GABAergic, glutamatergic synapse) and receptor loss [77]. In addition, GABAergic and glutamatergic synapse have been shown to be involved in the synthesis of heat shock proteins (HSP), which contributes to the repair of stress-induced synaptic protein damage and facilitates neuroprotective mechanisms [78, 79]. Cheruiyot et al. [80] showed that the glutamatergic synapse pathway is highly enriched for heat tolerance in Holstein cows [80]. Another pathway, circadian entrainment was identified for HS time-period 2, which is relevant since there is increasing evidence that thermoregulation is controlled by circadian rhythm and that animals with well-entrained circadian rhythms may be better able to cope with HS by optimizing their physiological responses to the stressor [81, 82]. In addition, Li et al. [83] used KEGG and GO enrichment analyses to demonstrate an association between circadian rhythm and heat stress response in Hu sheep.

Disease and development related pathways were mostly enriched for HS time-period 3, with particular emphasis on pathways related to cancer and cardiovascular disease, such as breast cancer and gastric cancer. Heat stress has been shown to activate HSP, which can promote cell proliferation and survival and may be involved in the development of cancer [84, 85]. In addition, as previously mentioned, HS increases the prevalence of cardiovascular diseases, potentially through the pro-inflammatory and pro-atherosclerotic activities of hormones, such as aldosterone [86]. Moreover, during HS time-period 3, sows were observed to suffer from the accumulated HS and produced more heat, leading to the activation of multiple pathways, including thermogenesis, hippo signaling pathway, mTOR signaling pathway, and Wnt signaling pathway. Machado et al. [87] showed that the thermogenesis pathway is related to stress-induced hyperthermia. Thermogenesis can also increase the capacity of endotherms to keep thermal homeostasis under cold stress conditions [88]. Hippo, mTOR, and Wnt signaling pathways are involved in various cell processes, such as cell proliferation, differentiation, and cell death. These pathways also regulate tissue homeostasis in multiple species [89, 90].

The main pathways identified for HS time-period 4 were related to death (nucleotide metabolism, pyrimidine metabolism, and purine metabolism) and hormone functions (e.g., renin secretion). During this HS time-period, sows maintained high body temperatures for several hours and suffered from HS for a longer period. Das et al. [91] showed that increased levels of pyrimidine and purine can activate a tolerance mechanism to protect nucleic acids and ultimately protein synthesis from HS. Stasolla et al. [92] found that alterations in pyrimidine metabolism could be an early signal of apoptosis and could be caused by an increase in endogenous nitric oxide (NO). HS increases the level of reactive oxygen species (ROS) and leads to dysfunction and destruction of cell membranes and finally causes cell death [93]. Both the pyrimidine and purine metabolism pathways have been reported to be enriched in response to HS in buffalo [94], beef cattle [95], and dairy goats [96]. It is not surprising that renin secretion was enriched for HS time-period 4, since animals may undergo severe dehydration during this stage, which can activate the renin–angiotensin–aldosterone system and cause hormonal responses to retain water and maintain mineral homeostasis [97]. The spliceosome pathway was also identified for this HS time-period, which is consistent with the report of Hu et al. [10] who showed that alternative splicing is an important transcriptional mechanism in heat stressed dairy cattle.

GWAS enrichment along the EG scale

Numerous HS-related pathways were identified based on the uniquely enriched genes for each EG class. Three pathways (other types of O-glycan biosynthesis, Fc gamma R-mediated phagocytosis, and signaling pathways regulating pluripotency of stem cells) were enriched exclusively for EG class 1. These are related to glycan metabolism, immunity, and stem cell. Heat stress induces the synthesis of HSP, which could affect stem cells, including their self-renewal, differentiation, sensitivity to environmental stress, and aging [98].

Five significantly enriched pathways were identified for genes that were uniquely enriched for EG class 2, and their functions were mainly associated with energy (i.e., carbohydrate digestion and absorption, mineral absorption, and bile secretion) and cell proliferation (e.g., oocyte meiosis, cell cycle). Carbohydrate is one of the most important energy sources and it is clear that HS has complex effects on carbohydrate metabolism through increased demands for glucose and changes in insulin signaling and sensitivity [99]. Bile functions are important for fat digestion and metabolism and contribute to growth and the intestine development via the hormones and pheromones that are excreted in bile. Supplementation with bile acids has been shown to reduce the harmful effects of HS in broiler chickens but this has not been studied in cattle or pigs [100]. Heat stress response markedly alters postabsorptive carbohydrate, lipid, and protein metabolism independently of reduced feed intake through coordinated changes in fuel supply and utilization by multiple tissues [99].

Nervous system related pathways (neuroactive ligand-receptor interaction and long-term depression) and nicotinate and nicotinamide metabolism were enriched for genes that were uniquely enriched for EG class 3. The nervous system plays a vital role in controlling and regulating body temperature by continuously monitoring body temperature and initiating appropriate responses, such as vasodilatation, sweating, and reduced metabolic rate, to maintain homeostasis [43]. In ducks, a neuroactive ligand-receptor interaction that is involved in maintaining energy homeostasis during HS has been identified in a previous study [101]. Studies in humans have demonstrated that HS induces mental issues such as chronic depression or chronic anxiety disorders [102, 103]. Nicotinamide provides anti-inflammatory and neuroprotective effects by increasing oxidative phosphorylation, buffering and preventing metabolic stress, and increasing mitochondrial size and motility [104, 105].

Taken together, the results suggest that the genes that have a dynamic function are important in regulating HS across various time-periods or EG classes and in the related pathways for HS in mammals. These findings are of interest for future research towards understanding and managing HS for coping with rising global temperatures. Greater knowledge of the biological basis of response to HS can help farmers and breeders design improved breeding programs and make selection decisions. However, heritability estimates for T_V are not available for the whole lactation period because of the data collection scheme implemented in our study. Such results could be useful and contribute to investigating how lactation affects body temperature under HS conditions. Moreover, more potential novel indicators of heat tolerance using body temperature data should be explored. As global warming intensifies, there is an increasing need for breeding more climatically-resilient livestock. Such animals should not only be minimally affected by potential disruptions due to climate changes and farm management but also demonstrate better health, welfare, and improved production efficiency. Thus, the next step is to derive novel indicators of climatic resilience in lactating sows based on variability in Tv and investigate the genetic relationship of these novel traits with production indicators and other routinely recorded traits in swine breeding programs.

Conclusions and implications

Automatically-recorded vaginal temperature is a promising indicator of HS response for breeding programs due to its moderate heritability and repeatability estimates. The RRM with Legendre orthogonal polynomials (i.e., LEG4) was identified as the best model due to its low BIC value and computational requirements. GxE interactions for T_V were identified, indicating that HS could result in re-ranking of breeding values under different climatic challenges or timing of HS. T_V is a highly polygenic trait and is controlled by multiple genomic regions of small effects. Twelve, 13, 16, and 10 relevant 10-SNP windows located on more than 10 Sus scrofa chromosomes (SSC) were estimated to explain more than 0.25% of the genetic variance for T_V across the four respective HS time-periods. For the three respective EG classes, 18, 15, and 14 10-SNP windows were identified for T_V. Two genomic regions, on SSC10 (59.370–59.998 Mb) and SSC16 (21.548–21.966 Mb), were identified across all HS time-periods and EG classes. Most genomic regions related to T_V have dynamic regulatory functions over time and under varying climatic conditions, and specific biological functions were identified for each HS time-period and EG class, such as those related to immunity, metabolism, and hormone. This study provides important insights to conduct genomic selection for improved heat tolerance in pigs, especially during lactation, based on automatically-measured T_V.