All datasets were provided by the Smithfield Premium Genetics company (Rose Hill, North Carolina, US). Phenotypic nucleus-herd records were obtained from January 2004 to December 2019 for Large White animals distributed among 33 farms located across North America. The farm latitudes range from 34° to 42° north, and the longitudes range from -77° to -113°. Different management strategies have been used across farms, such as equipment to alleviate heat stress. However, these practices changed over time, and this information was not recorded. Moreover, the systematic effects included in the statistical models (e.g., contemporary group–CG) are assumed to account for these effects. The traits analyzed were related to body composition [ultrasound backfat thickness (BFT; mm) and ultrasound muscle depth (MDP; mm)], growth [weaning weight of the piglets (WW; kg) and off-test weight (OTW; kg; approximately measured at around 5.5 months of age)], and reproduction [interval between farrow (IBF; days), total number of piglets born (TNB), number of piglets born alive (NBA), number of piglets born dead (NBD), number of piglets weaned (WN), and weaning to estrus interval (IWE; days)]. Contemporary groups (CG) were defined by the concatenation of farrowing year, farrowing season, and farrowing farm for the reproductive traits and of birth year, birth season, and birth farm for the growth and body composition traits, respectively. The phenotypic datasets were edited independently for each trait by removing records deviating ±3.5 SD from the mean. In addition, CG with less than 10 records were removed from further analyses. The descriptive statistics and the number of animals for each studied trait are presented in Table 2.

A total of 265,943 animals were included in the pedigree, which represented more than 10 generations. A total of 8,992 animals were genotyped using the PorcineSNP10K (8,652 SNPs for 886 animals), PorcineSNP50K (50,549 SNPs for 5,706 animals), PorcineSNP60K (57,019 SNPs for 865 animals), and PorcineSNP80K (64,577 SNPs for 1,535 animals) Bead Chips (Illumina, San Diego, CA, United States). In each SNP panel, animals with genotype call rate smaller than 0.90 were removed. Genotype imputation was performed using the FImpute v3 software (Sargolzaei et al., 2014). The missing genotypes were imputed first from the 10K to the 50K panel and then from the 50K or 60K to the 80K panel. The detailed imputation process was described in Chen et al. (2021). Quality control (QC) of genotype data consisted of removing SNPs with call rate below 0.90, minor allele frequency lower than 0.01, and difference between observed and expected heterozygous frequencies lower than 0.15. The genomic QC was implemented in the BLUPF90 family software (Misztal et al., 2018). In the end, 55,375 SNPs for 8,686 animals (7,017 female and 1,669 male) were included in the subsequent analyses.

Public weather station records for all farms were obtained from the Local Climatological Data at the National Oceanic and Atmosphere Administration (www.ncdc.noaa.gov/cdo-web/datatools/lcd?prior%2520=%2520N). Based on the location of the farms, climatic information was collected from the nearest airports. The average distance between an airport and the farm was 30 km (ranging from 7 to 64 km). Seven ENV were evaluated, including the average of mean temperature (MeanT), average of maximum temperature (MaxT), average of minimum temperature (MinT), dew point (DewP), average relative humidity (RH), average discomfort index (DI; Thom, 1959), and average THI calculated as in the Guide to Environmental Research on Animals (NRC, 1971). The average of the raw ENV was used as the environmental covariate. Supplementary Table S1 shows the description of the ENV range for TNB.

The critical periods evaluated in this study were chosen based on physiological knowledge (Johnson et al., 2015; Gonzalez-Rivas et al., 2020; Johnson et al., 2020) of heat stress and were trait dependent. The specific critical periods are shown in Table 1. Both chronic (i.e., HS during the entire grow-finish period; average of 120 days prior to measurement date) and acute (i.e., HS in the last 30 days prior to harvest) HS might have great impact in carcass-related traits, such as OTW, MDP, and BFT and, therefore, were evaluated for the mentioned traits. In addition, pigs exposed to in utero HS have been reported to develop a variety of postnatal phenotypes that prevent profitable production and compromised health and welfare in commercial production systems (Johnson et al., 2020; Maskal et al., 2020). These responses to HS are related to alteration of postnatal response, that is, change of core body temperature. With more intense HS, lactating sows have increased intestinal epithelial hyperpermeability and reduced reproduction (Mayorga et al., 2020). In utero heat stress can also impact long-term growth performance and decrease postnatal growth rates (Johnson et al., 2020; Maskal et al., 2020). Heat-stressed lactating sows also reduce feed intake as a mechanism to minimize metabolic heat production (Renaudeau et al., 2012; Williams et al., 2013; Cabezón et al., 2016). As such, in this study, two different critical periods were evaluated for IBF, IWE, WW, and WN: the first, a period ranging from the last stage of gestation (around day 80) until the farrowing date, and the second, a critical period ranging from the last stage of gestation throughout lactation until the weaning date. Regarding TNB, NBA, and NBD, only one critical period was evaluated. Early studies showed that the first weeks of gestation have a greater impact in the number of viable embryos (Edwards et al., 1968; Omtvedt et al., 1971; Wildt et al., 1975). Therefore, for TNB, NBA, and NBD, the average ENV ranging from 20 days before breeding to 30 days into gestation was considered.

Chen et al. (2021), using data from the same population, performed an evaluation of G×E interaction based on CGe. In brief, they estimated variance components and individual breeding values for seven traits (OTW, MDP, BFT, TNB, NBA, WN, and WW), regressing their phenotype on the average effect of CG (Chen et al., 2021). The results obtained in the current study will be compared to the previous results based on the following comparison criteria: theoretical accuracy of genomic estimated breeding values (GEBV) and rank correlation between the GEBV of individuals regressed on the recommended ENV and CGe.

The best statistical model to describe each trait was defined in two steps: (1) the systematic (fixed) effects were defined based on the backwards elimination procedure using the lm function available in the software R (R Core Team, 2020) and (2) the random effects for each trait were defined using the BLUPF90 software (Masuda, 2018; Misztal et al., 2018) by comparing models, including a different combination of random effects [i.e., animal genetic effect, animal permanent environmental effect, and common environment (litter) effect]. Model comparisons were made, and the significant random effects were defined based on the Akaike Information Criterion (Akaike, 1973). Specific and detailed information regarding the definition of the models can be found in Chen et al. (2021). The final effects for each trait are shown in Table 2.

After the definition of the significant effects to be included in the statistical model, RNM, using the single-step GBLUP approach, was implemented to obtain the reaction norm (RN) for each animal considering each ENV individually. Legendre orthogonal polynomials (order = 1; Kirkpatrick et al., 1990) were used to model the trajectory of phenotypic traits across environmental conditions. Variance components for all “trait × ENV” combinations were estimated using single-trait RNM and Bayesian inference, under a Markov chain Monte Carlo framework, using the THRGIBBS1F90 software (Misztal et al., 2002). Considering the corresponding effects for each trait (Table 2), the following model was used: y ik = α +   x i ′ β + ω φ ^ k + ∑ ( n 0 i + n 1 i φ ^ k )   + e ik , where y i is the phenotypic observation of animal i, α is the intercept, x i ' is the row incidence vector for β β , β is the vector of systematic effects described in Table 2, ω is the systematic regression coefficient of y i   on the ENV, φ ^ k is the ENV vector (expressed as first-order Legendre polynomial) at the value k, n 0 i and n 1 i are the RN intercept and slope of animal i regressed on φ ^ k for the random effect n n ∈ { a , pe ,  ce } ( n ∈ { a , pe ,  ce } , a being the animal genetic effect, pe the animal permanent environmental effect, and ce the common environment effect, as described in Table 2 for each trait), and e ik is the random residual for the animal i. The assumptions regarding the random effects are as follows: [ a 0 a 1 ] ∼ N ( 0 ,   H ⊗ [ σ a 0 2 σ a 0 a 1 σ a 0 a 1 σ a 1 2 ] ) , and [ pe 0 pe 1 ce 0 ce 1 e ] ∼ N ( 0 ,   I ⊗ [ σ pe 0 2 σ pe 0 pe 1 0 0 0 σ pe 0 pe 1 σ pe 1 2 0 0 0 0 0 σ ce 1 2 σ ce 0 ce 1 0 0 0 σ ce 0 ce 1 σ ce 1 2 0 0 0 0 0 σ e 2 ] ) , where σ n 0 2 , σ n 1 2 , and σ n 0 n 1 are the variance of coefficient  n 0 i , variance of coefficient n 1 i , and covariance between n 0 i and n 1 i , respectively, where n represents the random effects described above (i.e., a, pe, and ce), e is the residual variance, A is the pedigree-based relationship matrix, and I is an identity matrix. A chain containing a total of 600,000 iterations, with thin and burn-in of 60 and 300,000, respectively, was used for all traits except IBF. For IBF, a chain containing 900,000 iterations with thin and burn-in of 60 and 600,000 was used. These parameters allowed the model convergence for all traits analyzed in this study. The convergence was verified based on both graphical analyses and Raftery and Lewis criterion (Raftery and Lewis, 1992), both available in the Bayesian Output Analysis (Smith, 2007) package of the R software (R Core Team, 2020).

The same single-trait RNM previously mentioned in this study was used to perform ssGBLUP analyses by replacing the A by the H matrix (Misztal et al., 2009; Aguilar et al., 2010; Christensen and Lund, 2010). As the direct estimation of H is computationally demanding, H ⁻¹ was calculated directly as in Aguilar et al. (2010): H − 1 = A − 1 + [ 0 0 0 τ ( α G − β A 22 ) − 1 − ω A 22 − 1 ] , where A − 1 is the inverse of the numerator relationship matrix A , A 22 − 1 is the inverse of the subset of A related to the genotyped animals, and G − 1 is the inverse of the genomic relationship matrix G [calculated using the first method proposed by VanRaden (2008)]. For the construction of H, the default parameters from the BLUPF90 software (Masuda, 2018; Misztal et al., 2018) were used (i.e., τ = 1, ω = 1, α = 0.95, and β = 0.05). The parameters used for the comparison of the results are described in Section 2.3 and Section 2.5.

Three criteria were used to select the optimal ENV to analyze HTol based on the studied traits: (1) the ENV yielding the strongest G×E estimate as measured by the parameter σ a 1 2 σ a 1 2 , (2) the ENV yielding the highest accuracy of GEBV for the slope, and (3) the ENV yielding the highest deviation of GEBV per ENV (i.e., allowing to more easily differentiate tolerant and susceptible individuals). Weights of 0.5, 0.3, and 0.2 were given to each criterion (1, 2, and 3, respectively) to facilitate the ENV selection. The climatic variable resulting in the highest final value was selected as the best indicator to evaluate HTol in Large White pigs.

After obtaining the RN components, a genetic (co)variance matrix Γ n among values of ENV was calculated for each n random effect (described above) as in Tiezzi et al. (2020): Γ n =   Ф G n Ф ′ , where G n is the estimated (co)variance matrix between the intercept and slope terms for the corresponding n effect and Ф is a matrix of the number of rows equal to the number of unique values of the ENV and two columns (a vector of “1” and the standardized ENV). The heritability of a trait at each single value k of ENV (h ² _k) was calculated as follows: h k 2 =   Γ akk Σ ( Γ nkk ) + σ e 2 , where Γ akk is the additive genetic variance for the ENV k and Γ nkk is the variance for the n effects (i.e., a, pe, and ce) for the ENV k.

The theoretical accuracy of GEBV predicted for the slope and intercept of the trait in consideration of the animal i was calculated as follows: Acc it = 1 − SD ^ it 2 ( 1 + F i ) σ ^ at 2 , where SD i ^ is the posterior standard deviation of GEBV for animal i for the RN of  t , being t the intercept or slope terms, F i is the inbreeding coefficient, and σ ^ a 2 is the estimated variance of animal additive effect for the intercept or slope (Aguilar et al., 2020).

Genotype-by-environment interaction can be confirmed, among other factors, by a genetic correlation lower than 0.70 across the range of ENV (Mulder and Bijma, 2007). The following equation was used to calculate the genetic correlation between the range of ENV: r kk ′ = σ ^ u kk ′ σ ^ u k 2 σ ^ u k ′ 2 , where the covariance of additive genetic effects between ENV k and k ′ is σ ^ u kk ′ = σ ^ a 0 2 + σ ^ a 0 a 1 φ k + σ ^ a 0 a 1 φ k ′ + σ ^ a 1 2 φ k φ k ′ .

The approximated weighted genetic correlation between the studied traits was assessed based on the correlation between the GEBV. Animals with a GEBV theoretical accuracy for the intercept and slope terms lower than 0.30 were removed from the calculations. The following approach was used (Almeida-de-Macedo et al., 2013): r ^ g ( xy ) = ∑ w i ( x i − x ¯ ) ( y i − y ¯ ) / ∑ w i ∑ w i ( x i − x ¯ ) 2 ∑ w i × ∑ w i ( y i − y ¯ ) 2 ∑ w i , where x ¯ = ∑ w i x i ∑ w i and y ¯ = ∑ w i y i ∑ w i , w i is the reliability-based weighting of animal i and calculated as ( Rel xi × Rel yi ) / Rel xi × Rel yi ( Rel = Acc 2 ) ( Rel = Acc 2 ) , and x i and y i are the GEBV of trait x and y , respectively. The standard errors of the approximated genetic correlations were estimated as: SE = ( 1 − r ^ g ( xy ) 2 ) / ( n − 2 ) , where n is the number of animals with GEBV theoretical accuracy (for intercept and slope terms) greater than 0.30. The genetic correlations estimated between ENV and the approximated genetic correlation between traits were only calculated for the selected optimal ENV for each trait (chosen according to the three criteria previously mentioned). Genotype-by-environment interactions were considered large, moderate, or weak if the genetic correlations between environments were lower than 0.50, between 0.50 and 0.80, and higher than 0.80, respectively.

The descriptive statistics after the QC are shown in Table 2. The number of observations ranged from 6,059 for WN to 172,984 for TNB. Systematic and random effects for each trait were defined and are described in Table 2. In summary, the effect of sex, birth parity, contemporary group (growth), and weaning age divided into different classes was used for MDP (four classes), BFT (five classes), WW (five classes), and OTW (four classes); farrow parity, contemporary group (reproductive), and farrowing age divided into 11 different classes were used for TNB, NBA, NBD, WN, and IWE; and seven classes were used for IBF. Among the random effects, animal additive genetic effects were considered for all traits, animal permanent environmental effects for TNB, NBA, NBD, and IWE, and common environment effects for MDP, BFT, WW, OTW, TNB, NBA, NBD, and IWE.

The heritability (h ²) estimates across the range of the selected ENV for each trait are shown in Table 3 and the complete pattern in Figure 1. The average h ² ranged from 0.04 (IBF) to 0.42 (BFT_30). A similar pattern was observed for TNB (from 0.09 to 0.12), NBA (from 0.08 to 0.12), and NBD (from 0.04 to 0.09), where a slight increase in h ² is observed as the ENV values increase. The traits OTW_30 (from 0.21 to 0.47), WN_wd (from 0.07 to 0.31), and WW_wd (from 0.05 to 0.26) also had a similar pattern, with higher h ² in extreme ENV values and lower estimates closer to the mean. An opposite pattern of average h ² was observed for MDP_30 (from 0.27 to 0.31) and BFT_30 (from 0.38 to 0.47). As the values of ENV increased, the h ² for MDP_30 also increases, and for BFT_30 it decreases. The heritability estimates within each trait (i.e., considering all the ENV and critical periods evaluated) had similar patterns and are shown in Supplementary Table S5.

Sires were selected to have at least 30 offspring so that the phenotypes were distributed under a large range in ENV (i.e., records of daughters could not be concentrated under the comfortable environmental gradients). The reaction norms for the slope term of the additive genetic effect of the five most HT and the five most heat-susceptible (HSusc) sires are reported in Figure 2. Re-rankings of animals were clearly observed for NBA, WN_wd, WW_wd, IWE, IBF, and OTW_30, while no clear re-ranking was observed for NBD, MDP_30, and BFT_30. Meanwhile, there was a trend of re-ranking observed for TNB at the best environmental conditions. For TNB (Figure 2A), HT sires appear to have higher performance than the HSusc sires even when the temperature increases. For NBA (Figure 2B), IWE (Figure 2F), IBF (Figure 2G), and OTW_30 (Figure 2H), a clear distinction between groups was observed, with HT sires performing better in the highest temperatures, while the HS performed better in the lower temperatures. For WN_wd (Figure 2D) and WW_wd (Figure 2E), the HT sires were observed to have better performance under environments with higher humidity, while the HSusc sires had a better performance under drier environments. For MDP_30 (Figure 2I) and BFT_30 (Figure 2J), no re-ranking was observed. For all traits, even though some of them do not have a clear re-ranking, it is possible to observe that there are groups of animals performing better than the average population performance (Figure 2, black lines) in challenging environments.

In addition, the Spearman (rank) correlations are shown in Supplementary Figure S2 for combinations of intercept and slope terms between each pair of traits. The greatest rank correlations were between TNB and NBA, where between intercepts was equal to 0.917 and between slopes was equal to 0.838, indicating that the GEBV of TNB tends to increase when the GEBV of NBA increases. All other trait combinations had moderate to low rank correlations. Between MDP_30 and OTW_30, IBF and TNB, IBF and NBA, and between MDP_30 and TNB, the slope terms were equal to 0.216, 0.155, 0.188, and 0.113. Negative rank correlations between the slope terms were observed for WN-wd and TNB (−0.196), WN_wd and NBA (−0.185), BFT_30 and WW_wd (−0.149), and IWE and WW_wd (−0.141). The remaining rank correlations between trait combinations for the intercept and slope terms had low rank correlations and are presented in Supplementary Figure S2.

The genetic correlations of additive genetic effects across environmental gradients are shown in Figure 3 (for weak and moderate correlations) and Supplementary Figure S1 (for strong correlations). All estimated genetic correlations between environments considered just ENV values ranging from the 10th to the 90th percentiles, therefore excluding extreme ENV values. The exclusion of the extreme values was done to remove ENV with a low number of records, which might affect the estimation of the parameters at the extremes. As expected, the genetic correlation decreased gradually for greater differences among the environmental gradients. The lowest genetic correlation was observed for WN_wd (Figure 3; average = 0.75), even reaching negative values (-0.274). In addition to WN_wd, IBF (Figure 3; average = 0.79), WW_wd (Figure 3; average = 0.80), and IWE (Figure 3; average = 0.83) seem to be largely affected by G×E interaction. Moderate interactions were observed for OTW_30 (Figure 3; average = 0.88), NBA (Figure 3; average = 0.92), and TNB (Figure 3; average = 0.96). The highest genetic correlations were observed for NBD (Supplementary Figure S1; average = 0.98), MDP_30 (Supplementary Figure S1), and BFT_30 (Supplementary Figure S1), with both of the latter having an average of 0.99. Therefore, we concluded that WN_wd, IBF, WW_wd, and IWE are largely and OTW_30, NBA, and TNB are moderately affected by the G×E interaction. However, NBD, MDP_30, and BFT_30 were only mildly affected by G×E interaction.

Individuals with GEBV theoretical accuracies higher than 0.30 were used to calculate the weighted Pearson correlation among the RN intercepts and slopes between all studied traits considering the recommended ENV for each trait. The minimum number of selected animals and the lowest average accuracies were 5,718 and 0.4224 (95% CI of 0.4210–0.4236), respectively (Supplementary Table S6). The approximate genetic correlations between intercepts, intercept and slope, slope and intercept, and slope and slope for each pair of traits (considering the recommended ENV for each trait) are shown in Figure 4.

The greatest genetic correlations were observed between the intercept terms. The greatest positive correlations were observed between TNB and NBA, with 0.928, and between BFT_30 and OTW_30 (0.581) for the RN intercept. A moderate genetic correlation was also observed between TNB and NBD (0.328), TNB and WN_wd (0.305), and NBA and WN_wd (0.305). Moderate genetic correlations between intercepts were also observed between IBF and IWE (0.375), MDP_30 and OTW_30 (0.332), and BFT_30 and MDP_30 (0.248). For WW_wd, moderate to low genetic correlations were observed with IBF (0.121), OTW_30 (0.216), MDP_30 (0.234), and BFT_30 (0.188). The combinations for the remaining traits for the intercept term had low genetic correlations and are shown in Figure 4A.

Approximate genetic correlations between the intercept and slope terms indicate what happens to the performance of an animal for a specific trait when another trait is expressed under challenging heat conditions. In this matter, the greatest positive correlation was between TNB slope and NBA intercept (0.149) and between NBD intercept and WN slope (0.135). A similar positive correlation was also observed between WW intercept and OTW slope (0.123) and between the pairs of BFT slope with NBD intercept and IBF intercept with WN slope (both equal to 0.109). The WN slope and TNB intercept and the NBD slope and WW intercept had the same correlation value of 0.105. The greatest negative correlation was observed between the NBA slope and the TNB slope with WW intercept (−0.137 and −0.125, respectively). All other combinations between intercept and slope had values lower than 0.100 (Figure 4B,C).

The smallest genetic correlations were observed between the slope terms. This correlation indicates the relationship between the traits under a stressful environmental condition. The greatest positive slope × slope genetic correlation was observed for WW_wd and IBF (0.146). A negative correlation of -0.131 was observed for NBD and WN_wd and of -0.111 between WN_wd and TNB. All the other trait combinations had correlations under the absolute value of 0.100 and are shown in Figure 4D.

As indicated in the “Materials and Methods” section, three different metrics were used as criteria for the selection and recommendation of ENV and critical period to evaluate HTol in pigs from the Large White breed (one of the main maternal-line breeds). The first criterion, based on the ENV yielding the higher G×E interaction measured by the σ² _a1 term, was also used by other studies evaluating HS effects in reproductive and carcass traits in pure breeds and crossbred animals (Tiezzi et al., 2020; Usala et al., 2021). The second criterion was based on the theoretical accuracy of GEBV. Therefore, the climatic metric recommended had to not only have the strongest G×E but also provide more reliable GEBV estimates. Furthermore, the third criterion was based on the ENV yielding the greatest deviation of GEBV within the ENV value. The latest mentioned criterion was considered because, in this study, the goal was to use an ENV that allows a clear distinction (i.e., higher dispersion of breeding values) between HT and HSusc animals. A large variability in GEBV (i.e., higher standard deviation of GEBV within each ENV value) allows better discrimination of the animals, which facilitates the selection of the desired animals. Another metric that could also be evaluated was the deviance information criterion (DIC). However, in this study, the DIC was not included as another criteria for ENV selection due to criticisms found in the literature, such as the lack of consistency, weak theoretical justification, and overfitting (Spiegelhalter et al., 2014; Clarke and Clarke, 2018; Salles et al., 2019).

Climate records from the National Climatic Data Center weather stations were used to describe the environmental conditions experienced by the sows at different stages of ovulation, pregnancy, and post-pregnancy. Weather stations used in the current study were on average within 30 km from the farm. Freitas et al. (2006) estimated a correlation of 0.9 between on-farm weather data and weather station data even for weather stations more than 300 km away from the farm. Although these climate records are usually publicly available and have been successfully used in other worldwide studies focusing on sow HTol for different reproduction, growth, and carcass traits (Tummaruk et al., 2010; Wegner et al., 2014; Tiezzi et al., 2020; Usala et al., 2021), we acknowledge that the collection of within-barn environmental measurements is important and should be done when possible. However, the variables needed for this study were not available for this current study. It is also worth noting that, although there was a large variability in the climatic variables used in this study, they might not be representative of all the possible range of climatic conditions. Therefore, additional studies should be performed in other populations (e.g., breeds) and geographical regions to confirm our results.

In the current study, acute HS had greater impact on OTW, BFT, and MDP than chronic HS, which suggests that, during the period of chronic HS, the animal might be able to recover from the thermal stress suffered. For these traits (OTW, BFT, and MDP), the variance components for the slope term were higher considering a period of 30 days before measurement than 120 days. The GEBV accuracy estimates were also greater for the acute HS. For WW and WN, the period ranging from the last stage of gestation (around day 80) until the weaning date resulted in higher G×E interaction and GEBV accuracies. Therefore, the time period ranging from the last stage of gestation until weaning date is the recommended critical period to evaluate HT for WW and WN (referred here as WW_wd and WN_wd, respectively).

Several studies implemented different climatic variables as the environmental gradient covariate—for instance, some studies used heat load calculated as THI greater than 70, temperature, humidity, or THI (Fragomeni et al., 2016; Tiezzi et al., 2020; Usala et al., 2021). However, reports of the best environmental metric to be used when evaluating HTol in pigs from the maternal line are scarce. Therefore, we evaluated several ENV based on temperature and humidity to achieve a recommendation for each of the evaluated traits. The use of an index-based environmental variable is another way of carrying out the analyses, and it has been done in plants (e.g., Resende et al., 2020). However, in the current study, single ENV was used as environmental gradient. The use of RH as the climatic variable is the recommendation for NBD, WW_wd, WN_wd, MDP_30, and BFT_30. Maximum temperature is the recommended ENV for TNB, NBA, IBF, IWE, and OTW_30. Even though for OTW_30 the use of MaxT was recommended based on our findings, THI also resulted in a very similar G×E interaction. This choice was based on THI development as most THI are originally related to the body temperatures of cattle exposed to HS and have been adapted to other species or populations (Ingraham et al., 1979; Buffington et al., 1981; Bohmanova et al., 2007). Therefore, the use of THI in pigs is viewed with caution since different physiological and anatomical characteristics (i.e., rumen, ability to sweat, body mass, etc.) are found between cattle and swine that directly influence thermotolerance. For OTW, a high (0.79) rank correlation between GEBV estimates using THI and MaxT was observed. Therefore, due to the caution in interpreting the THI values for swine and the high rank correlation between THI with MaxT, MaxT was recommended as the ENV to evaluate HTol in maternal-line pigs for OTW. For five out of the 10 traits evaluated, the ENV selected (Table 5) were built upon RH as climate variable and the other five used MaxT. The preponderance of RH, rather than other climatic variables, in challenged sows was reported by Tummaruk et al. (2010) and Tiezzi et al. (2020) who studied reproductive performance under HS conditions. This finding might indicate that the barn cooling systems can better mitigate the impact of temperature than RH.

Different environmental gradients for the RN have been used over time (CGe: Li and Hermesch, 2016; Song et al., 2020; Chen et al., 2021 and climatic variables: Fragomeni et al., 2016; Tiezzi et al., 2020; Usala et al., 2021). However, the feasibility of using the average CGe or climatic variables as the environmental gradient for assessing G×E interactions is unknown. In this regard, Chen et al. (2021), using the same population and seven traits also used in the current study (OTW, MDP, BFT, TNB, NBA, WN, and WW), calculated the GEBV of animals using CGe as the environmental gradient. In general, the estimates based on the recommended ENV of the current study resulted in higher theoretical accuracies of GEBV, except for MDP, where the accuracy based on CGe was 0.580 (i.e., indicator of G×E interaction), while the accuracy based on the recommended ENV was 0.350 for the slope term. A Spearman (rank) correlation between the GEBV of individuals regressed on the recommended ENV and CGe shows a negative moderate rank correlation for BFT (−0.57), WW (−0.47), and MDP (−0.31) and a low positive rank correlation for the other traits (Supplementary Table S7). These results suggest that, despite being used and historically providing reliable breeding value estimates (Song et al., 2020; Chen et al., 2021), the use of CGe accounts for additional variability, including nutrition practices and management. This might not be the interest when evaluating HTol. Therefore, as suggested, the use of ENV might be more useful when the objective is to evaluate the breeding values of animals under challenging climatic stress.