Extensive genome-wide association analyses identify genotype-by-environment interactions of growth traits in Simmental cattle

Background Understanding the genetic basis of genotype-by-environment interactions (GxE) is crucial to understand environmental adaptation in mammals and improve the sustainability of agricultural production. In addition, GxE information could also be useful to predict the vulnerability of populations to climate change. Results Here, we present an extensive study investigating the interaction of genome-wide SNP markers with a vast assortment of environmental variables and searching for SNPs controlling phenotypic variance (vQTL) using a large beef cattle dataset. We showed that GxE contribute 10%, 4%, and 3% of the phenotypic variance of birth weight, weaning weight, and yearling weight, respectively. GxE genome-wide association analysis (GWAA) detected a large number of GxE loci affecting growth traits, which the traditional GWAA did not detect, showing that functional loci may have non-additive genetic effects between genotype classes regardless of differences in genotypic means. We also showed that variance-heterogeneity GWAA can detect loci enriched with GxE effects without requiring prior knowledge of the interacting environmental factors. Functional annotation and pathway analysis of GxE genes revealed biological mechanisms by which cattle respond to changes in their environment, such as neural signaling, metabolic, hypoxia-induced, and immune system pathways. Knowledge of these pathways will be important as climate change becomes a burden on animal health and productivity. In addition, ecoregion-specific GxE SNPs detected in this study may play a crucial role in identifying resilient and adapted beef cattle across divergent environments. Conclusions We detected novel trait associations with large GxE effects for birth weight, weaning weight, and yearling weight. Functional annotation and pathway analysis uncovered genomic regions involved in response to environmental stimuli. We unraveled the relevance and complexity of the genetic basis of GxE underlying growth traits, providing new insights into how different environmental conditions interact with specific genes influencing adaptation and productivity in beef cattle and potentially across mammals

and has been previously associated with stature and growth traits in many beef cattle 3 0 1 breeds [71][72][73]. In addition, we detected a GxG vQTL on chromosome 3 at 106 Mb 3 0 2 within TRIT1, a gene that was also reported to harbor a vQTL involved in epistatic  Using z 2 from MADE, 15 vQTLs were detected (four vQTLs for univariate BW; one 3 0 5 vQTL for univariate WW; nine vQTLs for univariate YW; and four vQTLs using 3 0 6 multivariate analyses, of which three vQTLs were also identified by the univariate 3 0 7 analyses), which likely identify GxE. Based on GO information, genes near or within  Using z 2 as dependent variables in the analyses allowed the estimation of the 3 1 9 proportion of residual variance that is explained by the SNPs in GEMMA (PVE). For 3 2 0 BW z 2 , WW z 2 and YW z 2 , SNPs explained 4.9% (SD = 0.8), 4.7% (SD = 0.8), and 5.9% 3 2 1 (SD = 1.1), respectively, showing that residual variance is heritable and could be changed 3 2 2 by selection, supporting previous studies [80][81][82]. Residual variability may reflect 3 2 3 differences between animals in uniformity for a certain trait, so the reduction of variance 3 2 4 could improve economic merit [80]. Rönnegard et al. [81] reported advantages in using selecting for reduced residual variance could be more beneficial for low heritability traits. Meta-analysis of ecoregion-specific GWAA 3 3 0 We performed a meta-analysis of ecoregion-specific GWAA (GWAA for each ecoregion 3 3 1 separately) to look for significant differences (Cochran's Q statistics) in effect size of 3 3 2 SNPs between the ecoregions, which are interpreted as GxE [83]. In addition, we 3 3 3 calculated the m-value statistic for the significant SNPs (Cochran's Q statistic's P-value < 3 3 4 1e-5), which is the posterior probability of an effect being present in an ecoregion given included in the meta-analysis due to their small sample sizes. Meta-analysis of ecoregion-3 3 8 specific GWAA detected 12 GxE SNPs for BW located on chromosome 6 (Additional  respectively). However, large differences in the SNP effect size (beta) were noticed 3 5 0 between ecoregions (Additional file 1: Figure S21B and C, respectively). We also performed traditional GWAA using univariate and multivariate linear mixed 3 5 4 models in order to compare with the GxE GWAA results. We detected major QTLs on  In addition, significant SNPs (P < 1e-5) were identified on chromosomes 3 (53.6 Mb) also identified in the univariate GWAA on chromosomes 5, 6, 7, 14, 17, and 20 in 3 6 1 addition to on chromosomes 4 (24.5 Mb), 9 (57.9 Mb), 13 (63.3 Mb), 19 (50.5 Mb), and 3 6 2 26 (45.4 Mb) (Additional file 1: Figure S22D and Additional file 11). All genes detected 3 6 3 by GWAA are included in the Additional file 12.

6 4
Out of 948 total significant SNPs using multivariate GWAA analysis for the growth 3 6 5 traits, 531 SNPs were located on chromosome 6 (from 26.8 to 48.8 Mb). Likewise, 3 6 6 several SNPs interacting with Forested Mountains ecoregion, elevation, precipitation and, 3 6 7 mainly, with the temperatures were also detected in the same positions on chromosome 6. Given the importance of chromosome 6 underlying the growth traits, we repeated the 3 6 9 multivariate GWAA including significant SNPs (e.g. rs109849093, rs464458177, 3 7 0 rs110305942, rs43459303, rs109278547, and rs137209027 located at 37,211,057; 3 7 1 37, 234,136; 37,418,164; 37,695,352; 38,139,495; and 41,084,765  including six significant SNPs (not in LD) simultaneously as fixed effects in the model, a  Comparing traditional GWAA with all GxE GWAA results, we verified that 312 SNPs with additive effects (QTL) also have GxE effects (mainly on chromosome 6, but also on 3 8 1 chromosomes 7 and 20). Some SNPs showed GxE effects even greater than their additive has an additive effect of -0.43 kg for BW in the total dataset, but its GxE effects is even 3 8 4 more impactful (-0.67 kg) on the BW of animals from the Forested Mountains ecoregion. The SNP rs109868461 (chromosome 6, at 36,902,704 bp) has additive effects of -0.69 kg 3 8 6 for BW in the total dataset, however it has a GxE effect of 1.11 kg on the BW of animals  CHCHD7, was also found as a direct effect QTL for BW, meaning that this locus affects GxE GWAA, which were found to interact with multiple continuous environmental suggesting that vQTL effects are much more enriched with GxE effects than additive 3 9 7 effects. However, GxE GWAA detected a greater number of GxE loci when compared  Only few GxE SNPs were detected using the meta-analysis of ecoregion-specific 4 0 0 GWAA (Cochran's Q statistic's P-value < 1e-5), which were also identified by the GxE  Prairie ecoregions could not be included in the analysis due to their small sample sizes, 4 0 3 suggesting that this was not an optimal approach to detect GxE associations for our 4 0 4 dataset. However, creating PM-plots (contrast between environment-specific P-value and  large GxE SNP effects indicates that there is considerable opportunity to improve 4 2 0 environmental resistance and performance in cattle. Our study also revealed the first 4 2 1 evidence of vQTLs controlling phenotypic variation in body weights due to GxG or GxE 4 2 2 effects in beef cattle. In addition, we uncovered several biological pathways affected by GxE that likely drive environmental adaptation. One of our major findings is that the ecoregions showed specific GxE with large 4 2 5 effects on the traits. In addition, ecoregion GxE GWAA detected GxE loci that were not identified by the continuous GxE GWAA. One possible explanation for these results is 4 2 7 2 0 that when analyzing ecoregions, we take into account not only climate and topography, 4 2 8 but also unmeasured variables like forage quality and local pathogens, as well as the 4 2 9 combination of them. This information allowed us to identify environmentally sensitive magnitude. We also found that major additive QTLs also carry GxE effects on growth  together may be lower than the power of the variance heterogeneity test, which may 4 5 1 explain why vGWAA identified GxE loci that were not detected by direct GxE GWAA.

5 2
We also observed that some vQTLs have differences in means across ecoregions 4 5 3 (Additional file 1: Figure S27), however, not all do. This suggests that some vQTLs supporting that phenotypic variability may be an adaptive evolutionary solution to an environmental factor of interest has been measured on all the genotyped individuals, 4 6 2 the direct GxE GWAA will be able to detect a larger number of GxE loci.

6 3
We found that environmental conditions modulate growth traits through genes that 4 6 4 participate in a variety of biological mechanisms that are activated to help the body cope reactions to deal with pathogens [94]. We also identified many metabolic processes that and have been reported to be involved in resilience [12,92,94], and adaptation processes  genetic effects between genotype classes regardless of differences in genotypic means.

8 8
We revealed biological mechanisms by which beef cattle respond to changes in their 4 8 9 environment. Neural signaling, metabolic, hypoxia-induced, and immune system 4 9 0 pathways were highlighted, indicating that climate change may become a burden on 4 9 1 animal health and productivity. In addition, ecoregion-specific GxE SNPs detected in this 4 9 2 study may play a crucial role in resilience and adaptation of beef cattle to divergent improving genetic gain, it only considers additive effects and, according to our findings, ecoregion-specific genomic predictions should be created to identify animals best suited 4 9 6 for given environments. This could help producers in making optimal decisions given 4 9 7 their geographical location. Our results revealed novel trait associations and alternative  Phenotypes were pre-adjusted for fixed effect of sex estimated using '--reml-est-fix' 5 1 2 option, and random effect of contemporary groups (CG) was predicted by the BLUP 5 1 3 method using '--reml-pred-rand' flag while controlling for population stratification using also included in the WW adjustment.

1 8
The DNA samples were genotyped with various low-density assays and were imputed SNP markers for 13,427, 11,847, and 8,546 animals for BW, WW and YW, respectively.

2 9
Eight continuous environmental variables were used for the GxE GWA analyses, Rainforest ecoregion were removed from the dataset, and those that were reared in the environmental variable measurements is shown in Additional file 2: Table S2. Variance components for pre-adjusted BW, WW, and YW were estimated using multi- incidence matrices relating observations to animals. Assuming that where # denotes the Hadamard product operation.

7 0
The variance of GxE interaction was also estimated considering the ecoregions as an Genome-wide association analyses (GWAA) and YW, respectively.

8 7
We also fitted a multivariate linear mixed model for the pre-adjusted BW, WW and  The GxE GWAA were performed using GEMMA through two approaches: using eight using multivariate linear mixed model were also performed. We performed vGWAA to detect locus affecting the difference in variance between standardized to z score by the rank-based inverse-normal transformation and squared (z 2 ), which is a measure of phenotypic variance [84,116]. The vGWAA were performed for 6 2 5 BW z 2 , WW z 2 , and YW z 2 using univariate and multivariate linear mixed models 6 2 6 implemented in GEMMA software as denoted in "Genome-wide association analyses" by SNP markers with r 2 greater than 0.6 [117]. univariate within-ecoregion GWAA into a single meta-analysis for each phenotype using 6 3 5 METASOFT software [118]. Only ecoregions with more than 1000 animals were 6 3 6 analyzed. In addition, the statistic m-value for the significant SNPs (Cochran's Q 6 3 7 statistic's P-value < 1e-5) were calculated, which is an estimative of the posterior performed (Additional file 2: Table S6), based on [120], as: where P is the nominal significant threshold (e.g., 1e-5); A is the number of SNP that 6 4 7 were significant at the nominal significant threshold; and T is the total number of SNPs number of significant SNPs is smaller than expected by chance, meaning that all SNPs  We explored 10 kb and 100 kb sequence windows that flanked the significant SNPs (P < 6 5 4 1e-5) to scan for GxE genes located in their vicinity and to identify possible regulatory 6 5 5 elements, respectively, based on the ARS-UCD1.2 bovine reference genomic positions. [122].     Increased responsiveness to intravenous lipopolysaccharide challenge in steers grazing         transcriptome. PLoS One. 2015;10:e0131459.  with Mycobacterium avium ssp. paratuberculosis infection status. BMC Genet.            association of changes in swine feeding behaviour due to heat stress. Genet Sel Evol.  embryos indicates species-specific regulation of genome activation. bioRxiv.   imprinting. Hum Mutat. 1997;10:155-9.     and VEGFR-2 in Bovine Placentomes from Implantation Until Term. Placenta.    Reference-based phasing using the Haplotype Reference Consortium panel. Nat Genet.           100 kb from significant SNPs (P < 1e-5).