Genome-wide scan highlights the role of candidate genes on phenotypic plasticity for age at first calving in Nellore heifers

Age at first calving (AFC) plays an important role in the economic efficiency of beef cattle production. This trait can be affected by a combination of genetic and environmental factors, leading to physiological changes in response to heifers’ adaptation to a wide range of environments. Genome-wide association studies through the reaction norm model were carried out to identify genomic regions associated with AFC in Nellore heifers, raised under different environmental conditions (EC). The SNP effects for AFC were estimated in three EC levels (Low, Medium, and High, corresponding to average contemporary group effects on yearling body weight equal to 159.40, 228.6 and 297.6 kg, respectively), which unraveled shared and unique genomic regions for AFC in Low, Medium, and High EC levels, that varied according to the genetic correlation between AFC in different EC levels. The significant genomic regions harbored key genes that might play an important biological role in controlling hormone signaling and metabolism. Shared genomic regions among EC levels were identified on BTA 2 and 14, harboring candidate genes associated with energy metabolism (IGFBP2, IGFBP5, SHOX, SMARCAL1, LYN, RPS20, MOS, PLAG1, CHCD7, and SDR16C6). Gene set enrichment analyses identified important biological functions related to growth, hormone levels affecting female fertility, physiological processes involved in female pregnancy, gamete generation, ovulation cycle, and age at puberty. The genomic regions highlighted differences in the physiological processes linked to AFC in different EC levels and metabolic processes that support complex interactions between the gonadotropic axes and sexual precocity in Nellore heifers.

In tropical beef cattle, animals are raised in a variety of production systems and exposed to a wide range of nutrition levels. These differences may cause genotype-environment (GxE) interaction on various productive traits [1][2][3] . GxE interaction occurs when the genetic variance and/or the classification of animals change according to the environment 4,5 . The re-ranking of animals caused by GxE interaction is an indication that some genomic regions might have different effects according to the environments 6,7 .
Integrating genotypic information with reaction norm (RN) models lead to the evaluation of GxE interaction more accurately compared to only pedigree-based analysis 3,8,9 . In addition, RN models can be combined with genome-wide association studies (GWAS) allowing to uncover genomic regions involved in animal adaptation by estimating SNP effects across environmental conditions 8,10 . In this context, a greater understanding can be obtained on how genetic variants are associated with reproductive traits and the environmental sensitivity by physiological adaptation.
The integration of GWAS with RN models has been shown to be a powerful approach to investigate the key regulators and genomic regions that explain the phenotypic variation across environments 10,11 . Moreover, Statistical modeling. Environmental condition descriptor. The dataset used to evaluate the sensitivity of sexual precocity to environmental variation belonged to commercial herds with high diversity of management and different environmental conditions, since these commercial herds were located in regions with variations on annual precipitation from ~700 to ~3000 mm and exhibiting a dry season that may last up to 7 months. Remarkably, differences in nutritional levels by farms are common, for instance in some of the farms, the animals received protein and mineral supplementation, especially during the dry season, while in others only urea supplementation was offered. Details of the climate classification in the Midwest, Southeast, and Northeast of Brazil can been seen in Alvares et al. 14 .
The AFC genetic sensitivity to environmental changes was assessed through the reaction norm model 15,16 . In this context, the animal's response to environmental condition changes was expressed as a function of a continuous environmental condition (EC) descriptor. In the lack of environmental condition descriptor information (e.g. temperature-humidity index, level of production, etc.), descriptors based on the CG solutions from phenotypic information, can be used. In this framework, the EC descriptor used to encompass the production level was based on the best linear unbiased estimates (BLUE) solutions of the contemporary group for yearling body weight (YBW), once it condenses the different management and environmental factors in which the animals were raised. We focused on YBW because the differences in production environments affecting YBW present an important effect on heifers' early sexual puberty 2,3,17, 18 .
The environmental descriptor used on RN model referred as environmental condition (EC), was the standardized BLUE of CG effects solutions for YBW obtained using an animal model considering the single-step GBLUP (ssGBLUP) method as follows: where y is the vector of YBW; β is a vector with the fixed effects of CG (defined by animals born in the same year and season, and raised in the same farm and management group from birth to yearling) and age at recording as a linear co-variable; a is a vector of additive genetic effects assumed normally distributed N H (0, ) a 2 σ ⊗ , in which σ a 2 is the additive genetic variance and H is the pedigree-genomic relationship matrix ⊗ is the Kronecker product and e is a residual vector assumed σ ⊗ N I (0, ) e 2 , where I is the identity matrix, and σ e 2 is the residual variance. The X and Z are known as incidence matrices related to fixed and additive genetic effects, respectively.
The H is a matrix that combines pedigree and genomic information 19 and its inverse (H −1 ) is given by: where y ij is the vector of AFC information of the animal i recorded in the level j of EC, F CG is the fixed effect of CG, ω f are the f-th fixed regression coefficients (intercept and slope) on Φ ( ) j , α fi are the random regression coefficients for additive effects of intercept and slope corresponding to animal i on EC level j, and e ij is a random residual. Residual variances were considered to be heterogeneous across EC levels and five classes were determined using K-means clustering variance 21 . The residual class 1: EC level lower than −1.5; residual class 2: −1.5 ≤ EC level < −0.5; residual class 3: −0.5 ≤ EC level < 0.; residual class 4: 0.0 ≤ EC level < 1.5, and residual class 5: EC level higher than 1.5.
The ssGRN model was fitted considering the following assumptions for the random effects a = {a j } ∼ N(0, between intercept and slope, respectively, and R is a diagonal residual variance matrix considering 5 heterogeneous classes. Estimates of (co)variance components were obtained by restricted maximum likelihood using the AIREMLF90 software 22 .
Estimates of SNP effects and variance explained. Before GWAS analyses, animals with GEBV accuracies (based on prediction error variance) lower than 0.40 for reaction norm parameters (intercept and slope) were excluded. So, a total of 2,550 heifers and 1,023 sires were considered for estimating SNP effects for the slope and the intercept. Clustering analysis was performed to assess the possible occurrence of population substructure in the population, based on the genomic relationship matrix. Results did not indicate population stratification ( Supplementary Fig. S4). The SNP effects (û kEC ) were estimated across EC levels using the following equation: , where α k is the vector of intercept and slope estimates for the k-th SNP, and Φ ′ f is the transposed vector of the f-th Legendre polynomials for each EC level. The SNP effects were obtained in three EC descriptors Low (−3.0 sd), Medium (0.0 sd) and High (3.0 sd). The genetic, phenotypic and residual variances, as well as heritability estimates for AFC across EC levels obtained with the RN model, are in Supplementary Fig. S5.
The percentage of genetic variance explained by the SNPs was estimated according to Wang et al. 23  , where z kEC is the z-score for SNP markers effects for each of three EC level. P-values for the SNP effects were computed as is the cumulative function of the normal distribution for the z-score. The p-values were corrected for multiple tests using the false discovery rate (FDR) 25  www.nature.com/scientificreports www.nature.com/scientificreports/ cance threshold (p-value < 0.01) used to consider a SNP as significant. The assessment of inflation/deflation factor were calculated as λ = median (p-value)/0.456, where values of λ between 1.0 and 1.1 were considered acceptable in GWAS 26 .
The linkage disequilibrium (LD) analysis was performed for the shared chromosome region between the three EC levels, by computing the r-square (r 2 ) 27 values for pairwise top SNP marker effect with markers around this region using the Gaston R package 28 .
Gene mapping of significant SNP in each EC level. The SNP effects estimated for AFC in the three EC levels, Low (EC = −3.0), Medium (EC = 0.0) and High (EC = 3.0), were deemed significant when −log 10 (p-value) > 6.0 (5% FDR). The Ensembl gene 94 database 29 and BioMart R package 30 were used to search genes located harboring 200 kb region of each significant SNP marker, using as reference the Bos Taurus UMD 3.1 assembly 31 . The Cattle QTL database (QTLdb) 32 was used to identify previously detected quantitative trait loci (QTL) overlapping these regions, also considering the UMD v3.1 assembly sequence as the reference map. Biological mechanisms and pathways (Gene Ontology -GO) involving the candidate genes were identified using the clusterProfiler R Package 33 , separately for each EC level, considering as background the Bovine database 34 . The association of a given gene set with AFC in each EC level, with GO terms was assessed using a hypergeometric test 35 considering a false discovery rate (FDR) of 5% for multiple test. The significance of GO terms was calculated as described by Boyle et al. 35 .  Table S1). Some of these SNP markers identified were shared among environmental conditions and can explain part of the genetic correlations for AFC among Low, Medium and High EC levels ( Fig. 1B). Medium and High EC levels shared most of the significant SNP markers (20 SNPs). A total of 2 and 3 SNP markers were shared by Low x Medium and Low x High, respectively, and 8 markers were significant for all three EC levels (Fig. 1A).

Results and Discussion
The higher number of shared SNP markers for AFC between Medium and High EC levels could, in part, explain the higher genetic correlation coefficients estimated for AFC between these two EC levels ( Fig. 1). When GxE interaction occurs, genotypes react differently according to the environmental levels 3,7 . Indeed, different genomic variants showed a specific effect on AFC in different environmental conditions (Fig. 1A). The genetic correlation for AFC in the three EC levels indicated an important environmental sensitivity resulting in a potential re-ranking of heifers under restrictive environmental conditions (Fig. 1B). These environmental differences are often associated with heat stress and seasonally poor nutrition, which causes a reduction in length of the estrous cycle, progesterone concentration and the developmental capacity of oocytes 36,37 .
The SNP markers detected (−log 10 (p-value) > 6.0) for Nellore AFC showed re-ranking of their effects across EC levels, in which the effects in the Low condition were different from those in the High condition ( Fig. 2A). The SNP effects changed in magnitude and direction of their effects according to the EC level ( Fig. 2A). Some studies in dairy cattle 7,38 , pigs 8 and beef cattle 3,10 have shown that different environmental conditions can cause substantial changes in SNP effect estimates.
Higher percentages of the total genetic variance were explained by SNP markers in the High (18.13%) compared to the Low (6.63%) and the Medium (12.67%) environmental levels, as also observed for SNP marker effect dispersion (Fig. 2 A,B). The differences in the proportion of total genetic variance explained by SNP marker across the environmental levels indicated an important change of genetic variance effects according to EC levels. Hence, unraveling the importance of genomic regions through GWAS may contribute to design a more precise strategy for the selection of complex traits 38 under different environmental conditions. Including genomic regions having environment-dependent sensitivity in statistical methods represents an additional biological insight to obtain more accurate genomic predictions and thus select Nellore heifers with greater genetic potential for sexual precocity and tolerance to harsh conditions. Streit et al. 11 , in dairy cattle, and Mota et al. 10 , in beef cattle, www.nature.com/scientificreports www.nature.com/scientificreports/ concluded that significant genomic regions across environments were environmental-dependent and associated with key physiological processes affecting the evaluated trait.
Significant regions surround genes with EC level-specific effect. Inflation-factor (λ) estimates are 1.05, 1.08 and 1.09 for the Low, Medium and High EC levels ( Supplementary Fig. S6), respectively and showed that the deviation of the observed test statistics from the theoretical quantiles was acceptable according to Georgiopoulos and Evangelou 39 .
Specific genomic regions affecting AFC in each of the EC levels were identified as well as regions shared by two or three of the levels (Fig. 3). The specific regions for the Low EC level were mapped on BTA 1, 6, 15, 17 and 27, for the Medium EC level on BTA 3, 5 and 18 and for the High EC level on BTA14. The regions along the genome that were associated with AFC for each EC (Fig. 3), imply in different physiology mechanisms leading to sexual precocity in a specific environmental condition 40 .
The variants identified with a specific effect in the Low EC level overlapped the QTLs previously associated with body weight, tridecylic acid content, daughter pregnancy rate, and calving ease (Supplementary Table S2 Table S3). These specific genes associated with Low EC level are involved in energetic metabolism by action on glucose homeostasis (GO:0042593, KCNJ11 and ABCC8) and negative regulation of insulin secretion (GO:0046676, KCNJ11 and ABCC8) ( Fig. 4 and Supplementary Table S4). Such findings support the hypothesis that production systems on pasture conditions under Low EC level exhibit an important factor for delaying the precocity onset, through its impact on the concentration of circulating metabolites 18,36,41,42 .
Some candidate genes identified in Low EC level on BTA15 (KCNJ1, KCNC1, MYOD1 and SERGEF) play a key role on response to estrogen (GO:0043627), regulation of estradiol (GO:0032355), hormone levels (GO:0010817), negative regulation of secretion (GO:0046888), muscle adaptation (GO:0014888 and GO:0043501) and are directly associated with growth rate, i.e. average daily gain ( Fig. 4 and Supplementary Table S4). According to Rosales Nieto et al. 43 , sexual precocity is influenced by growth rate, which affects muscle and fat deposition. The MYOD1 gene is a component of myogenic regulatory factors (MYF-5, MYOD, myogenin and MRF4) which is associated with muscle metabolism 36 . The myogenic factors are associated with endocrine factors, e.g. GH, estrogen, and IGF, with an important role in the regulation of muscle mass, fiber size, nutrient partitioning, and reproduction 44,45 . These factors might directly affect sexual precocity due to the effect of myogenic factors action on the insulin-glucose metabolic homeostasis with major effects on metabolic and endocrine roles on lower nutritional levels 36 .
The genes ABCC8 (BTA15-35.64 Mb), DKK2 (BTA 6-19.49 Mb), and NUCB2 (BTA15-35.64 Mb) are involved in energy and protein metabolism. The NUCB2 gene is involved in hypothalamic pathways regulating feed intake and energy homeostasis, acting on leptin 46 . The ABCC8 gene is associated with glucose homeostasis (GO:0042593), negative regulation of insulin secretion (GO:0046676) and negative regulation of peptide hormone secretion (GO:0090278) ( Fig. 4 and Supplementary Table S4) with an important effect on regulating glucose metabolism by insulin secretion 47 . The DKK2 gene increases adipogenesis and insulin resistance through on inhibition of the Wnt signaling pathway 48 . These metabolic responses to Low EC level involve the catabolic process, mainly to energy mobilization and protein degradation.
When animals are raised under Low environmental conditions, genes related to metabolism are required to recover the adaptability to poor nutrition conditions. The gene set identified in Low EC levels plays an important role to keep metabolic homeostasis, leading to different growth rates according to heifer adaptability allowing that precocious heifers reach their genetic potential. Ferraz Jr et al. 18 observed that Nellore heifers selected to sexual precocity in the low nutritional level were able to attain the lowest age at puberty compared to heifers with high AFC, mainly by differences in IGF1 and Leptin hormone level. In Low EC, physiological changes affecting levels of www.nature.com/scientificreports www.nature.com/scientificreports/ insulin and glucose could affect the AFC, because the oocyte quality and development of both oocyte and embryo respond directly to these metabolic inputs 40,49 .
The SNP markers with a specific effect on AFC in Medium EC level harbor QTL regions associated with body weight, calving ease, calving interval, heifer pregnancy, and marbling score (Supplementary Table S2). A total of 4 genes (ILDR2, POGK and TADA1) were associated with these SNP markers on BTA 3 (2.11-2.12 Mb) and they are associated with important biological factors ( Fig. 5 and Supplementary Table S5).
The ILDR2 gene on BTA 3 (2.12 Mb) participate in a biological process that affects the circulating insulin and glucose levels and increases subcutaneous adipose tissue ( Fig. 5 and Supplementary Table S5). The ILDR2 gene is associated with lipid metabolism and plays a critical role in lipoprotein assembly, mainly on glucose and calcium homeostasis, in response to metabolic stress 50 .
The genes identified in Medium EC levels could reduce the AFC by their action in insulin, glucose levels, and energetic metabolism. It is likely that the improvement in EC levels affects heifers' sexual precocity through changes in metabolic substrates (glucose and insulin) and their levels rather than a direct effect on reproductive hormones 17,37,49 . Samadi et al. 17 observed that improved nutrition increases the blood concentration of insulin and glucose in Brahman heifers allowing lower AFC. The insulin and glucose provide signals in reproductive pathway regulations, affecting the oocytes quality and development, and the modulation of the gonadotropin-releasing hormone (GnRH) secretion 40 .
In High EC level, 14 SNP markers with specific effect for AFC were identified and surrounding two genes OPRK1 (BTA14 23.39 Mb) and TMEM68 (BTA 14 24.84 Mb). These genes are associated with neurons development related to the glial cell, a key component to the responses to growth factors and GnRH (gonadotropin-releasing hormone, GO:0032276 and GO:0032274) ( Fig. 6 and Supplementary Table S6). It affects www.nature.com/scientificreports www.nature.com/scientificreports/ the integral functional elements of the synapses, responding to neuronal activity and regulating synaptic transmission and plasticity causing a reduction in the release of LH (GO:0032275) from the pituitary gland 51,52 affected by OPRK1 gene. The major impact of OPRK1 in Nellore AFC is highlighted by the important biological process on negative regulation of LH secretion, which had been pointed out as a potential candidate in Brahman 53 and Nellore 54 cattle puberty. The TMEM68 gene was identified to be associated with many reproductive traits in Nellore and Brahman 55 and showed up-regulation in heifer blastocyst 56 , and is likely involved in energy metabolism and lipid turnover 57 . In addition, the sexual precocity is associated with physiological events linking major metabolic factors for the attainment of puberty 49 .    Tables S5 and 6).

Shared regions surround genes for AFC on Medium and
The genes ADAMTS4 on BTA3 (8.34 Mb) and ADAMTS18 on BTA18 (4.63 Mb) affect folliculogenesis and ovulation, and they are also required for normal gonadal morphogenesis and function in cattle 58 . The genes ADAMTS4 and ADAMTS18 are associated with follicle growth because their expression are induced by the follicle-stimulating hormone (FSH), with a key role on tissue morphogenesis during embryonic development 58,59 . The FCER1G gene on BTA3 has been associated with QTLs related to pregnancy rate at first service, services per conception, and days open (Supplementary Table S2). This gene affects the immune response and presents an important effect on reproductive functions 60 .
The gene set (ARHGAP30, USF1) identified on BTA3 showed an important effect on the control of several genes related to lipid and glucose metabolism 61,62 . Hence, these genes can cause variation in AFC in Nellore heifers due to their influence on metabolic homeostasis, by direct action on the reproductive process, e.g. ovarian follicles, oocytes quality and embryos 18,49 . The genes USP21 and UFC1 are members of the Ubiquitin family gene and they were related to the metabolic process (GO:0008152), affecting the proteolysis regulating that represents a key aspect for the cellular function 63 .
The DEDD gene on BTA3 plays an important effect on the cell-cycle regulatory process, presenting roles in SGK1 activity and AKT protein stability, with an essential effect on corporal glucose homeostasis 64 . The DEDD gene is indispensable for the support of female fertility in mice 64 and it presents an effect on biological GO involved in female pregnancy associated with response to estrogen (GO:0060135) and reproductive process (GO:0022414). Glucose metabolism has shown an association with cattle reproduction 41,65 . In this context, the highest effect of genomic regions with a key role in glucose metabolism could lead to changes in reproductive endocrine hormones, e.g. LH and FSH 66 .
The SNP marker identified on BTA 5 (9.46-9.47 Mb), surrounding the PPP1R12A gene which plays a central role in a wide variety of biological processes. The PPP1R12A gene, also known as myosin phosphatase target subunit 1 (MYPT1), is involved in insulin signaling regulation by their action on insulin receptor substrate-1 (IRS-1) 67 . The effect of the PPP1R12A gene in sexual precocity could be related to metabolic homeostasis by action in insulin and glucose with an important role in the nervous system and ovary 68 .
The SNP markers found on BTA5 (10.85-10.88 Mb; Supplementary Table S3) are involved with energy and lipid metabolism by the action of the ACSS3 gene that stimulates the utilization of acetate generating acetyl-CoA used in lipid synthesis or energy production 69 . The PPFIA2 gene on BTA5 10.86 Mb members of liprin family show interaction with members of LAR, which have an important role in axon guidance and mammary gland development 70 . These genes are associated with the integrity of the central nervous system and represent an important role in the regulation of many aspects associated with the reproduction process 70 . www.nature.com/scientificreports www.nature.com/scientificreports/ The gene TERF2IP on BTA 18 (3.02 Mb) is associated with various metabolic pathways, mainly regulating the NF-kB pathway 71 . The NF-kB pathway has a key role in the control of energy production by the regulation of glycolysis and cell respiration 72 .
The regions shared by Medium and High EC lead to differences in AFC through their action on energy metabolism. Energy metabolism is a pathway whereby higher EC levels lead to sexual precocity in Nellore cattle. The association of energy metabolism in sexual precocity occurs mainly because the glucose is the major energy source required to ovarian function and the luteinizing hormone (LH) secretion 49 . Brickell et al. 73 observed that an increase of glucose levels in Holstein-Friesian heifers reduced their age at first breeding and calving. In Brahman heifers, the improvement of metabolic status increased glucose and lipid concentrations reducing the AFC 17 .

Shared regions surround genes for Afc on low, medium and high ec levels.
Results of the GWAS for AFC pointed out to shared genomic regions on BTA 2 (105.03-105.17 Mb) and BTA 14 (24.82-25.07 Mb) with the highest peak corresponding to marker rs137780934 located at 24.94 Mb on BTA14. The genomic regions surrounding the BovineHD0200030238 (rs443442023; BTA2 -105.05 Mb) and BovineHD1400007242 (rs137780934; BTA14 -24.94 Mb) markers have been pointed as a shared region that plays a key role in genetic differences in reproductive traits. These regions affect mostly growth and reproductive pathways, insulin growth factors (IGF) and hormonal levels [74][75][76][77] . Genomic regions on BTA14 (20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) were identified by Fortes et al. 74 harboring genes with shared effects on growth, carcass and reproductive traits, and showed a putative functional mutation. In the same region, Melo et al. 55 (Fig. 7). Previous QTLs studies detected genomic regions on BTA 2 and 14 play an important role in multiple traits (Supplementary Table S2) 74,[78][79][80][81] . These genomic regions are known to act on age at puberty, calving to conception interval, the interval from first to last insemination, body weight and intramuscular fat (Supplementary  Table S2). Shared genomic regions for AFC in three EC levels harbor candidate genes on BTA2 (IGFBP2, IGFBP5, SHOX and SMARCAL1) and BTA14 (LYN, RPS20, MOS, PLAG1, CHCD7 and SDR16C6) with a striking effect in biological processes that might help explain the variability in the sexual precocity . Taking into account the gene network analyses from GWAS across the three EC levels (Low, Medium and High) gave additional insights into the complex relationship among the specific and shared genes (Fig. 8). This result highlights the key role for molecular pleiotropy associated with the shared and specific genes across the EC level, which are likely to play an important role in the genetic architecture for sexual precocity in heifers raised under harsh conditions (Fig. 8).
These genomic regions on BTA 2 and 14 have been associated with hormone secretion (GO:0048545, GO:0060416, GO:0060986 and GO:0010817) and regulation of gonadotropin release (follicle-stimulating hormone -FSH and luteinizing hormone -LH). They are associated also with negative regulation of gonadotropin secretion (GO:0032276) and ovulation cycle (GO:0042698), affecting the estrous cycle by changes in ovarian follicle growth. Additionally, both regions affect endocrine hormone secretion (GO:0060986) with a key role in regulatory pathways of sexual precocity leading to a metabolic resilience 78 . The metabolic resilience to environmental conditions could explain their relationship with sexual precocity affecting factors related to growth rate and body condition.
Genes mapped on BTA14, particularly the gene set (PLAG1, CHCHD7, LYN, MOS, PENK and RPS20), highlight the effect on responses to the endogenous stimulus (GO:0009719; response to the hormone, response to fibroblast growth factor and response to transforming growth factor-beta). Furthermore, they have been associated with growth hormone (GO:0060416; growth hormone receptor signaling pathway) and response to hormone stimulus, resulting in a change of cellular activity and insulin levels (GO:0032868; insulin receptor signaling pathway). Such regions were also related to IGF1 (insulin-like growth factor 1) hormone levels due to polymorphisms on BTA14 82 that have been significantly associated with the reduction in IGF1 levels, fat depth and beginning of puberty in females and males [80][81][82][83] . Indeed, a link between the growth and reproductive traits in beef heifers has been previously reported by a response to these biological factors found in this study 17,18,42,49,74 .
The genomic region on BTA14 (BovineHD1400007242; rs137780934) has a key role on the puberty onset through GH signaling 84 . Hence, the PLAG1 gene can affect the beginning of puberty in heifers through GH signaling and its direct effect on IGF1 levels 85 . In addition, the PLAG1 region was associated with delayed puberty by the increased growth, which affects body structure and growth rate, leading to delays in reproductive precocity in heifers until an adequate body condition to body size relation is achieved 86 .
The genomic regions on BTA 2 harbor the gene set (IGFBP2, IGFBP5, SHOX and SMARCAL1). The genes IGFBP2 and IGFBP5 have an important role in glucose metabolism of cattle, by the regulation of the bioactivity of IGF1 and IGF2 in the ovary 41 . In this sense, the IGFBP has been pointed out as a key factor to control follicle growth by sensitivity to gonadotropins 87 . Perhaps, the key function of IGFBP2 and IGFBP5 are tied to AFC across EC levels through a hormonal link that involves GH, IGF1, glucose and insulin. These metabolic signals affect the sexual precocity through direct effects on ovarian cells as well as on gonadotropins secretion 66 .
The gene set on BTA 2 (IGFBP2 and IGFBP5) and on BTA 14 (PLAG1, CHCHD7 and LYN) could have their action linked by physiological mechanisms affecting the reproductive traits by their effects on major pathways, such as insulin, glucose and IGF systems (Figs. 4-6). They are an important signal that allows associating the reproductive events to nutritional conditions. These factors are associated with the regulation of reproductive hormones and interact with ovarian activity, embryo development, oocyte production, and quality in cattle 77,[88][89][90] .
The GWAS results for AFC in different environmental conditions showed strong evidence of genomic regions affecting physiological events, with changes on circulating metabolic signals, which result in the activation of the hypothalamus-pituitary-adrenal (HPA) axis. Thus, it is possible that differences in EC levels affect physiological response on AFC. According to Rhimd 91 and Evans et al. 92 , harsh conditions and climate changes have Figure 8. Network of candidate gene identified within 200 kb for SNP markers significantly affecting AFC in three environmental condition (Low, Medium and High). The gene network was built from known proteinprotein interactions (edges) between gene products (nodes) using the string database for Bos Taurus. The node color represents the shared or specific genes across the Low, Medium and High environmental condition.
www.nature.com/scientificreports www.nature.com/scientificreports/ an important effect on neuroendocrine mechanisms, changing the patterns of reproductive efficiency affecting puberty, estrus, and ovulation.

conclusion
Combining genome-wide scan and reaction norm models helped to identify genomic regions associated with age at first calving (AFC) in Nellore heifers in different environmental conditions. These genomic regions showed strong environmental dependence. Genomic regions identified in Low EC level BTA 1, 6, 15, 17 and 27 confirmed the importance of metabolic homeostasis on sexual precocity in restricted environments. In Medium EC levels the genomic regions on BTA 3 and 18 highlighted the role of energy and lipid metabolism involved in the reproductive pathway. In the High EC level, the genomic region on BTA14 acts on metabolic substrates. The significant genomic regions suggest that common variation in AFC trait on Medium and High EC is essential by controlling a key role in the regulation of physiological mechanisms of reproductive precocity. The shared genomic regions on BTA 2 and 14 in Low, Medium and High EC levels are directly implicated in the regulation of reproductive pathways underlying endocrine parameters associated with precocity. Overall, the major genomic regions uncovered across EC levels shared or with a specific effect are related to major modulators for metabolic adaptations to different environmental conditions allowing the heifers' sexual precocity.

Data availability
The Nellore phenotypic and genotypic information are not publicly available because they belong to commercial breeding programs. The data are available for academic use from the authors upon reasonable request.