Genotype-by-environment interaction in Holstein heifer fertility traits using single-step genomic reaction norm models

Abstract

Background: The effect of heat stress on livestock production is a worldwide issue. Animal performance is influenced by exposure to harsh environmental conditions potentially causing genotype-by-environment interactions (G × E), especially in highproducing animals. In this context, the main objectives of this study were to (1) detect the time periods in which heifer fertility traits are more sensitive to the exposure to high environmental temperature and/or humidity, (2) investigate G × E due to heat stress in heifer fertility traits, and, (3) identify genomic regions associated with heifer fertility and heat tolerance in Holstein cattle.

Results: Phenotypic records for three heifer fertility traits (i.e., age at first calving, interval from first to last service, and conception rate at the first service) were collected, from 2005 to 2018, for 56,998 Holstein heifers raised in 15 herds in the Beijing area (China). By integrating environmental data, including hourly air temperature and relative humidity, the critical periods in which the heifers are more sensitive to heat stress were located in more than 30 days before the first service for age at first calving and interval from first to last service, or 10 days before and less than 60 days after the first service for conception rate. Using reaction norm models, significant G × E was detected for all three traits regarding both environmental gradients, proportion of days exceeding heat threshold, and minimum temperature-humidity index. Through single-step genome-wide association studies, PLAG1, AMHR2, SP1, KRT8, KRT18, MLH1, and EOMES were suggested as candidate genes for heifer fertility. The genes HCRTR1, AGRP, PC, and GUCY1B1 are strong candidates for association with heat tolerance.

Conclusions: The critical periods in which the reproductive performance of heifers is more sensitive to heat stress are trait-dependent. Thus, detailed analysis should be conducted to determine this particular period for other fertility traits. The considerable magnitude of G × E and sire re-ranking indicates the necessity to consider G × E in dairy cattle breeding schemes. This will enable selection of more heat-tolerant animals with high reproductive efficiency under harsh climatic conditions. Lastly, the candidate genes identified to be linked with response to heat stress provide a better understanding of the underlying biological mechanisms of heat tolerance in dairy cattle.

Keywords: Genotype-by-environment interaction; Heat stress; Heifer; Reaction norm; Single-step GWAS.

Fig. 1

Reproductive events and the definition of critical period in heifers. The red rectangle represents the critical period, defined as the time period for which heifers are likely to suffer from heat stress. AFC = age at first calving, IFL = interval from first to last service, CR = conception rate of first service

Fig. 2

Heritabilities estimated based on reaction norm models with the matrix H for different traits using a prop-EG or b mTHI-EG as environmental gradient. For a, the x-axis is the proportion of days exceeding the threshold with a range of 0 to 1; while for b, the x-axis is the minimum THI with a range of 15 to 75

Fig. 3

Genetic correlations estimated by reaction norm models (RNMs) with the matrix H. The color indicates the magnitude of the genetic correlation. a Correlations between different levels of prop-EG estimated from RNM under S1. The x-axis and y-axis are the proportion of days exceeding the threshold, ranging from 0 to 1. b Correlations between different levels of prop-EG estimated from RNM under S2. The x-axis and y-axis are the proportion of days exceeding the threshold, ranging from 0 to 1. c Correlations between different levels of mTHI-EG estimated from RNM under S1. The x-axis and y-axis are the minimum THI, ranging from 15 to 75. d Correlations between different levels of mTHI-EG estimated from RNM under S2. The x-axis and y-axis are the minimum THI, ranging from 15 to 75

Fig. 4

The re-ranking plots for gEBVs of sires. The blue and red lines represent sensitive and resilient sires, respectively. a Re-ranking plots for three traits estimated using prop-EG under S1. The x-axis is the proportion of days exceeding the threshold with a range of 0 to 1 and y-axis is gEBV of sire. b Re-ranking plots for three traits estimated using prop-EG under S2. The x-axis is the proportion of days exceeding the threshold with a range of 0 to 1 and y-axis is gEBV. c Re-ranking plots for three traits estimated using mTHI-EG under S1. The x-axis is the minimum THI with a range of 15 to 75 and y-axis is gEBV. d Re-ranking plots for three traits estimated using mTHI-EG under S2. The x-axis is the minimum THI with a range of 15 to 75 and y-axis is gEBV

Fig. 5

Percentage of the intercept and slope genetic variances explained by a sliding window of 20 SNPs for three fertility traits, which were estimated under scenario one of prop-EG. The x-axis is autosome segments; the y-axis represents the proportion of explained variances; the grey horizontal lines are thresholds (top 0.5%) for candidate genomic regions; and different color sets for the less relevant genomic markers indicate different traits

Fig. 6

Percentages of the intercept and slope genetic variances explained by a sliding window of 20 SNPs for three traits, which were estimated under scenario two of prop-EG. The x-axis is autosome segments; the y-axis represents the proportion of explained variances; the grey horizontal lines are thresholds (top 0.5%) for candidate genomic regions; and different color sets for the less relevant genomic markers indicate different traits

Fig. 7

Trajectories of SNP effects changing over EGs. The x-axis is environmental gradient; the y-axis represents the SNP effects; the vertical bar is the standard deviations of SNP effects at each level of EG; the blue lines indicate scenario one; red lines indicate scenario two; and, different color sets indicate different clusters. a Trajectories of SNP effects changing over prop-EG. The x-axis is the proportion of days exceeding the threshold with a range of 0 to 1. b Trajectories of SNP effects changing over mTHI-EG. The x-axis is the minimum THI with a range of 15 to 75

Fig. 8

Number of shared candidate genes for each EG in different traits and clusters. C1 = SNP effects changes in preferential ways (decrease for AFC and IFL; increase for CR); C2 = SNP effects changes in the opposite ways (increase for AFC and IFL; decrease for CR); C3 = constant SNP effects over time