All phenotypic records were collected at the performance test (PT) station of the National Breeders Association of Rendena cattle—ANARE, Trento Italy (www.ANARE.it). All phenotypes belonged to young (on average of 1 month of age) candidate bulls. About 60 young bulls are tested every year at the PT station for a total period of 11 months, following the criteria reported in Mantovani et al. (1997). Records have been collected since 1985, when PT started, until present times. The phenotypes collected during the PT are body weight (BW), average daily gain (ADG), carcass fleshiness (CF) and dressing percentage (DP). Both CF and DP are evaluated in vivo by 3 skilled operators at the end of the PT period and averaged to obtain the final score. The CF evaluation applies the same scores of post-mortem carcass appraisal established by the European Union Council (SEUROP), where the middle class (R) is equal to 100 points and other classes (upper or lower classes) correspond to 10-points-variations. Furthermore, the evaluation also considers sub-classes (e.g., R+ and R-for the middle class) that are spaced 3.33 points from the class score. DP is a visual prediction of the post-mortem measure of DP: the operator makes a visual appraisal of the individual at the end of the performance test, offering an estimate of the expected DP—i.e., conformation—at slaughter (Mantovani et al., 1997). Average daily gain (ADG) is calculated as the linear regression of weight (BW) on age. For this study, ADG and BW were collected at different stages of PT. ADG has been divided into ADG_i and ADG_f: ADG_i covers the daily gain of the first half of the testing period (since entering the PT station until the 6th month), while ADG_f covers the daily gain of the second half (from the 6th month to the end of the period). ADG covering the entire PT test was labeled as ADG_tot. BW was split along the same timeline as ADG: body weight at the entrance to the station (BW_i), at 6 months (BW_m) and at the end of PT (BW_f). Data cleaning consisted of removing animals with a regression of weight on age showing a coefficient of determination below 0.9 (for further details, see Guzzo et al., 2019).

The biological material of the animals chosen for the genotyping resulted from salivary swab, hair (at least 30 bulbs), or ear tissue from biopsy brand, collected by ANARE on females and young candidate bulls at PT, as well as from semen of proven bulls, already subjected in the past to PT and progeny test for milk and to a large extent now eliminated. The Bovine 150K Array GGPv3 Bead Chip (HD, 138,974 SNPs), and Illumina Bovine LD GGPv3 (LD, 26,497 SNPs), were used for genotyping (Illumina, Illumina Inc., San Diego, CA, United States). The overlapping between the two panels is about 60%. The HD platform was used for 554 young bulls, while 1,416 individuals (174 males and 1,242 females) were genotyped with LD chips. To achieve a reliable genomic imputation accuracy, the 174 males were animals with at least one parent and one half-sib genotyped with HD chips. The genotyped females were individuals with a kinship of at least 0.2 with phenotyped animals.

Before proceeding with imputation, we performed a preliminary quality control removing SNPs with a minor allele frequency (MAF) < 0.01 and call rate lower than 0.90, using Plink program (Purcell et al., 2007). Only the 29 autosomal chromosomes (BTA) were used for association, and progeny conflicts were fixed using the seekparentsf90 program (Aguilar et al., 2018).

AlphaImpute2 was used for imputation (Whalen and Hickey, 2020), as it combines a population imputation algorithm (Positional Burrows Wheeler Transform) with pedigree-based imputation (iterative peeling); we used the same parameters as in Mancin et al. (2021a). The accuracy of the imputations was roughly estimated as a correlation between true and imputed SNPs. To this aim, ten rounds of cross-validation were performed: in each round the overlapping SNPs between the two panels were removed in ten animals and then imputed using the HD panel from young bulls as reference population (Supplementary Table S1). Subsequently, the correlation between the true and the imputed genotypes was calculated on these animals.

After imputation, we performed a second genomic quality control with the preGSf90 program (Aguilar et al., 2018): the SNPs with MAF lower than 0.05 and SNPs that deviated too much for the expected value of heterozygosis (i.e., Hardy-Weinberg Equilibrium) were removed. In accordance to Wiggans et al. (2012) the threshold for was set to 0.15: SNPs were deleted if | n   o f   h e t e r o z i g o i s   n − 2 p q | > 0.15 . In addition, SNPs with a call-rate < 0.90 and animals with a call rate < 0.95 were excluded. The final genomic database contained 1,690 animals (698 with both genotypic and phenotypic information), and 113,279 SNPs. Genome-wide linkage disequilibrium (LD) within chromosome was also calculated, as the squared correlation of allele counts for two SNP. Principal Components Analysis (PCA) of G matrix and LD were also calculated with pregsGSf90.

Single step genome-wide association (ssGWAS) models were used to estimate allele substitution effect. In ssGWAS, the estimation of allele substitution effects was obtained from a linear transformation of the BLUP of breeding value under ssGBLUP model (Aguilar et al., 2019). Mancin et al. (2021b) showed the advantages of this method in terms of QTL detection and control of populations structure over two-step methods in which de-regression of breeding value as pseudo phenotype is required. This issue is particularly evident in the presence of unbalanced data (i.e., sex-limited traits). In fact, the ssGWAS allows the use of both male and female genomes even when analyzing a phenotype collected only in individuals of one sex.

The ssGBLUP model used in this analysis, written in matrix form, is the following: [ X ′ X X ′ Z Z ′ X Z ′ Z + H − 1 σ e 2 σ a 2 ] [ b ^ a ^ ] = [ X ′ y Z ′ y ] (1) Where phenotypes are included in vector y, X is the incidence matrix of fixed effects (group of contemporaries, cow parity class and months of birth), b is the vector of these effects. The contemporary group has 147 levels, with each level consisting of bulls grouped together at the Performance Test because homogeneous by age (i.e., born within 1 month of each other; 82.

Animals per group on average, minimum 5 and maximum 142). The parity order of cow has four classes (first parity; second parity; third to seventh parity; above the eighth parity), and the classes of months of birth correspond to the single months, as in Guzzo et al. (2019).

Z represents the incident matrix that relates the random genetic additive effects to the phenotype, with effects represented by vector a . The vector of random residual error (e) has a normal distribution N ( 0 , I σ e 2 ) , where σ e 2 is the residual variance. In the ssGBLUP vector of additive genetic effects is distributed as N ( 0 , H σ a 2 ) , where σ a 2 is the additive genetic variance and H is the (co)variances structure which combines pedigree and genomic relationships (Aguilar et al., 2010). Its inverse, used in Eq. 1 is described as: H − 1 = A − 1 +   [ 0 0 0 G − 1 − A 22 − 1 ] (2) where A − 1 and A 22 − 1 are the inverse of the pedigree kindship matrix respectively for all animals and for only genotyped animals. Since the frequencies of current genotyped population are used to center G and pedigree and genomic matrices have different bases, G was adjusted so the average diagonal and off-diagonal matches the averages of A ²² . Pedigree kinship (sub) matrix was estimated tracing back the pedigree up to 7 generations, i.e., 6,644 animals. G matrix was built using the methods proposed by VanRaden (2008), as follows: G 0   = MM ′ 2 ∑ p i   ( 1 − p i   ) (3) where M is a matrix of SNP content centered by twice the current allele frequencies, and p i   is the allele frequency for the ith SNP (VanRaden, 2008).

Additionally, to avoid singularity problems, the final G was computed as G = λ G 0 + β I   (4) Where G is the matrix present in the Eq. 2, I is an identity matrix of the same dimensions, λ and β are two weighting coefficients, with λ = 0.99 and β = 0.01. These values were chosen due to their influence on the power of signal detection of the GWAS, and because they resulted in inflation close to optimum values. In addition, G was adjusted to a better blending with diagonal and off-diagonal of A²² as described in Vitezica et al. (2011): δ = 1 − 0.5 n 2 ( ∑ i ∑ j A 22 ( i , j ) −   ∑ i ∑ j G i , j ) (5)

Then, the vector of estimated breeding values was obtained as: g ^ = λ δ 1 2 ∑ p q M ′ G − 1 a ^ 22 (6) Where a ^ 22 is the vector of estimated breeding values of genotyped animals. The prediction error variances g ^ , necessary to calculate the p-values, were calculated following Gualdrón Duarte et al. (2014) and computed as in Aguilar et al. (2019), where: V a r ( g ^ ) = V a r ( λ δ 1 2 ∑ p q M ′ G − 1 a ^ 22 ) (7) V a r ( g ^ ) = λ δ 1 2 ∑ p q M ′ G − 1 Var ( a ^ 22 ) G − 1 M   λ δ 1 2 ∑ p q (8)

Since Var ( a ^ 22 ) is equal to PEV ( a ^ 22 ) − var ( a 22 ) ; thus Var ( a ^ 22 ) = G σ ^ a 2 − C 22 . It follows that formula Eq. 8 becomes: V a r ( g ^ ) = λ δ 1 2 ∑ p q M ′ G − 1 ( G σ ^ a 2 − C 22 ) G − 1 M   λ δ 1 2 ∑ p q (9) C 22 is a submatrix of C belonging to the genotyped animals and represents the prediction error variances of a ^ 22 . The p-values are then calculated as p - v a l u e i = 2 ( 1 − Φ ( | g ^ i sd ( g ^ i ) | ) ) (10) Where g ^ i is the allele substitution effect of SNP i and s d ( g ^ i ) represents the square root of Eq. 9, Φ (∙) is the cumulative density function (CDF) of the normal distribution.

Two thresholds were used for the association tests: a genome-wide 5% significant level of −log10(p) = 5.55 (0.05/17,766) and a suggestive association with −log10(p) = 4.29 (0.1/17,766). These are the thresholds corrected for multiple tests i.e., p n where p is the probability level of significance and n is the corresponding number of independent SNPs (n = 17,766) calculated using the “poolR” R package (https://cran.r-project.org/web/packages/poolr; R Core Team, 2021), according to Li and Ji (2005). The number of independent tests was calculated based on the number of eigenvalues. Instead of the standard approach of Cheverud (2001), we used the approach by Li and Ji (2005), a function that decomposes the eigenvalues in the integral part (Effective Number Independent Test) and the nonintegral part.

The (co)variance components have been estimated with REML using Average-Information algorithm (Gilmour et al., 1995). Approximate standard error of (co)variance components has also been estimated through Monte Carlo sampling as in Houle and Meyer (2015), in which standard deviations were calculated from Monte Carlo chains sampled from multinormal distribution with covariance being the inverse of the Average Information Matrix and the estimated variances as the expectation. Then the heritability for the 3 phenotypes was calculated under single trait models as in Eq. 1. Heritability was calculated as: h 2 =   σ a 2 ( σ a 2 + σ e 2 )   ; where σ a 2 and σ e 2 are, respectively, the additive genetic and the residual variances.

Genetic and phenotypic correlations were estimated with bi-traits models, which are equivalent to Eq. 1 except for the animal additive genetic and residual variance, assumed to follow a multivariate normal distribution with mean 0 and variances G ⊗ H, and R ⊗ I, where G = | σ a 1 2 σ a 1 a 2 σ a 1 a 2 σ a 2 2 | ;   R | σ e 1 2 σ e 1 e 2 σ e 2 e 1 σ e 2 2 | (11) where G is the matrix of additive genetic (co)variances σ² _a1, σ² _a2, σ_a1a2 of traits 1 and 2, R the matrix of residual (co)variances σ² _e1, σ² _e2 and σ_e1e2 of traits 1 and 2.

The correlation was estimated as: c o v =   σ i 1 i 2 ( σ i , 1     σ i , 2 ) where i stands for the genetic and phenotypic correlation; 1 and 2 refer to the different performance test traits, and σ i 1 i 2 is the covariance between traits 1 and traits 2, off diagonal of Eq. 11. For phenotypic (co)variance, we mean the sum of the genetic and the phenotypic (co)variances. Traits that do not include zero in their correlations Higher Posterior Density Interval (HPD) were declared significantly correlated.

All the genomic analyses were carried out with BLUPF90 family software (Aguilar et al., 2018) following the procedure described in Lourenco et al. (2020). Manhattan plots were drawn using “ggplot” R package (Wickham, 2016), as were the LD graphs.

Pathway’s enrichment analysis was conducted to identify which biological pathways and functional elements were enriched for the investigated traits. From GWAS results, we selected SNPs with nominal p-values of < 0.01 which were mapped to genes based on a distance of 15 kb from the coding region using the “biomaRt” R package (Drost and Paszkowski, 2017) and Bos taurus UMD3.1 assembly as in Pegolo et al. (2020). Functional enrichment analysis was carried out on the list of significant genes using the Cytoscape plugin ClueGo (Bindea et al., 2009). As functional categories, we used cellular component, biological process, and molecular functions within the Gene Ontology (GO, http://geneontology.org) database and the Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/). The Benjamini-Hochberg correction was applied to declare significant pathways: only pathways showing FDR < 0.05 were retained. The minimum number of genes in the pathway was set to 3; the minimum percentage of genes present in the pathway was set to 4%. To simplify the redundance of GO terms we provide figures with similar terms grouped based on their semantic similarity using the R packages “rrvgo” (Sayols, 2020). In addition, we investigated if the candidate regions declared as significant by our GWAS overlapped with QTL in animal QTLdb, identified with R package “GALLO” (Fonseca et al., 2020).

Descriptive statistics after data editing of the phenotypes are shown in Table 1. Phenotypic and genetic correlations and the heritability (h²) for the analyzed traits are reported in Table 2. Body weight traits presented an average value of h² lower than other traits: BW_i showed the lower heritability (0.130), while BW_m and BW_f had heritability of 0.220. In fact, as reported in literature, a large discrepancy of values has been observed for heritability of body weights, and generally, traits similar to birth weight or weaning weight have a slightly lower heritability than weight measured in more advanced stages (Yin and König, 2018). Average daily gain (ADG_tot) presented an intermediate heritability of 0.322 partitioned into 0.164 and 0.220 for ADG in the first and last period. As for body weight, ADG presents lower h² in first stages of the performance test, and h² values agree with what has been found in the literature (Yin and König, 2018). The highest heritabilities were found for the traits related to the carcass conformation, with a value of 0.45 and 0.47 respectively for CF and DP, close to what was observed in other local dual-purpose or beef cattle (Albera et al., 2001; Sbarra et al., 2013; Mancin et al., 2021c). These traits also appeared highly genetic correlated. All ADG traits were moderately genetically correlated with them, with a value of 0.5 on average. On the contrary, body weight measured at the beginning of the performance test was not significantly correlated with CF and DP. Interestingly, the weights measured in more advanced periods showed an increase of genetic correlation with a value close to 0.7. Body weight and ADG also presented a strong genetic correlation with body weight traits, especially for the traits measured at the final stages of the performance test. In terms of genetic correlations, the results agree with what was found in other local dual-purpose or beef breeds (Veselá et al., 2011; Filipčík et al., 2020). Phenotypic correlation followed the same trends of genetic correlation but with a lower magnitude (Table 2, under diagonal).

A homogeneous density distribution (number of SNPs per Mb) was found throughout the genome, apart from few relatively small blank areas in 12 chromosomes. For further details on SNP density on each chromosome after imputation and quality control, see Supplementary Figure S2. The new imputed panel had a SNPs density close to the one found in the young bulls genotyped with HD platforms. A value of imputation accuracy of 0.95 ± 0.05 was observed via cross-validation in the HD males (Supplementary Figure S2). Combined with the high correlation between the A and G matrix, these results confirm the reliability of the new Alphaimpute2 algorithm for this population.

The PCA scatterplots (Figure 1) illustrate a homogenous distribution of allele frequencies in individuals that comprised our study population. No stratification has been observed in the first two components, suggesting that most G matrix variance is explained by many eigenvalues with small effect. Genome-wide linkage disequilibrium and MAF have also been explored since the availability of high-density SNP platforms permits to explore the LD decay at an unprecedented resolution. In addition, MAF and LD are useful for understanding differences in population history and demography and for its impacts for genome-wide mapping studies. LD decay per each chromosome is reported in Supplementary Figure S3. As expected, most tightly linked SNPs presented strong levels of LD while it rapidly declines when the distance increases. A within-chromosome LD average value of 0.19 ± 0.12 has been observed. When the distance between markers is lower than 1 Mb, the LD squared correlation between pairs of loci across autosomes (r²) (Hill and Robertson, 1968) reached an average value of 0.17 ± 0.27, and when the distance was > 1 Mb LD decreased to 0.04 ± 0.09 (Supplementary Figure S3). Larger levels of LD have been observed for chromosome 6 (0.20), while lower levels of LD were observed for chromosome 28 (0.18). An average value of 0.29 ± 0.12 was observed for minor allele frequency; no noticeable difference has been observed along the 29 chromosomes, with MAF values ranging from 0.28 ± 0.12 (chromosome 12) to 0.30 ± 0.12 (chromosome 19). With respect to the other local Italian breeds (i.e., Fabbri et al., 2020), Rendena presents a lower level of LD. This issue implicitly underlines the reassuring demographic situation of Rendena compared with other indigenous cattle of Italy, as it demonstrates a lower risk of inbreeding depression.

The full results of GWAS are reported in Table 3. We found a total of 8 SNP significantly associated with 5 of the investigated traits, and 47 SNPs suggestively associated with all 7 investigated traits (Figure 2). Pathway analysis revealed that out of 113,279 SNPs, 77,506 were located within a 15 kb window of annotated genes; in the end, 14,380 annotated genes were used as a background for each trait. On average, 628 genes near significant SNPs (<0.01) were identified and subsequently used for pathway analysis of each trait. All traits presented an inflation factor close to optimum values of 1 (Figure 2) calculated based on the median chi-squared test. In addition, analysis on localized linkage disequilibrium (0.5 Mb form significant SNP), has been carried out (Figures 3–7), and results indicated that all significant candidate genes are extremely close to the significant SNPs, except for candidate gene ZNF784, which is situated between two significant SNP (Figure 6).

The main region of interest for both CF and DP traits was situated on a gene-rich region of BTA18, between 55 and 62 Mb, where 3 significant and 9 suggestively associated SNPs allowed to locate several candidate genes (Figure 6; Table 3). The QTL with the highest significance for CF (suggestively associated for DP) was located within candidate gene LOC513941 (Figure 7), translating into a cationic amino acid transporter 3-like. This type of transporters regulates the metabolism of cationic amino acids, a key factor for growth and beef characteristics in cattle (Liao et al., 2009). Further corroboration of the importance of this metabolic pathway for CF was the enrichment of 10 GO terms (Figure 10; Supplementary Figure S4H), within the group of “amino acid transport,” such as amino acid transmembrane transporter activity (GO:0015171), and amino acid transmembrane transport (GO:0003333).

A second SNP in the same region (significant for DP and suggestively associated for CF; Table 3) was located within gene CCDC155 (Coiled-coil domain containing 155). This gene encodes for a protein involved in dynein complex binding and actin filament organization and it has been associated with beef conformation (Lemos et al., 2016; Hardie et al., 2017). Apart from being the main component of the cytoskeleton, actin constitutes together with myosin the myofilaments, which grant muscle cells their mobility and thus ultimately their organization and dynamics. The association of actin filaments and carcass traits was again made apparent also by the number (more than 30) and diversity of enriched GO terms related to actin (Figure 10; Supplementary Figures S4G,H): for example, those related to GO:0098858 (CF), actin-based cell projection; GO:0030048 (CF and DP), actin filament-based movement; GO:0070161 (CF and DP), anchoring junction; GO:0030833 regulation of actin filament polymerization; GO:0005912 (CF and DP), adherens junction (Londoño-Gil et al., 2021). Similarly, for DP 20 terms were enriched for pathways associated with actin filament-based GO terms (Supplementary Figure S4G).

In the same region of BTA18, our analysis found two more candidate genes with a known association with size and growth traits, all with one or more suggestively associated SNPs for CF. Siglec-5 is a gene commonly found in GWAS concerning cattle size and growth traits; its over-expression indicates a deficiency of leptin, and thus longer gestation time and bigger fetuses (Hardie et al., 2017). KLK12 is a kallikrein gene, a serin protease associated with food intake and feed efficiency at the transcript level in backfat and rumen (Kern et al., 2016). LOC101904435 and ZNF784 are zinc-finger proteins: the former is suggestively associated with both CF and DP; the latter only with CF but is linked to food intake in cattle (Olivieri et al., 2016).

Finally, three more SNPs (one significant both for CF and DP and two SNPs suggestively associated for DP) were situated within NLRP2 gene (NACHT, LRR and PYD domains-containing protein 2), a key player in early embryogenesis, maternal effects, immune response, and inflammasome (Peng et al., 2012).

Taken together, these results about carcass traits have numerous substantial implications. Firstly, we highlight how the 57–62 Mb region on BTA18 can truly be considered a hotspot of genetic diversity in this breed (as it is for several others; Grigoletto et al., 2020; Purfield et al., 2020). Secondly, as expected with strongly correlated traits, CF and DP shared part of their genetic architecture, as significant SNPs for the two traits are mostly in the same region. Only another region was shared, as both traits reported two suggestively associated SNPs close to each other on BTA28 (Table 3). The region encompasses the PHYHIPL gene, which influences feed efficiency (Abo-Ismail et al., 2018), whose link with carcass traits has recently been established (Seabury et al., 2017).

CF was associated only with two more SNPs, one on BTA12 and the other on BTA14 (Table 3). While the former was more than 1 Mb far away from any annotated functional element, the latter fell within SAMD12, a gene already found to have a significant dominance signal to chuck roll and be associated with 18-months weight in Simmental (Zhuang et al., 2020). On the other hand, DP had an almost significant signal on BTA1: the gene closest to the SNP was SIM2, already known to be associated with carcass quality, differentiation of longissimus, and semimembranosus muscle (De Las Heras-Saldana et al., 2019; Edea et al., 2020). To conclude, the strongest of the remaining suggestively associated signals for DP came from BTA4, within LOC112446424, a non-coding RNA close to candidate gene SLC13A4, a cationic canal important both for muscle traits in sheep and growth and development in cattle (Carvalho et al., 2020; Kaur et al., 2020).

While, as we mentioned, results from pathway analyses (represented in Figure 10; Supplementary Figures S4G,H), and GWAS were often complementary, pathway analyses for both CF and DP resulted in the enrichment of a robust number of pathways related to neuron activity, not really pointed out by GWAS results. Such pathways referred to the regulation of neuroblast proliferation (GO:1902692 for CF), chemical synaptic transmission (GO:0007268 for CF), neurogenesis (GO:0022008 for CF and DP), neuron projection (GO:0043005 for DP), synapse (GO:0045202 for DP) and especially synaptic transmission, glutamatergic (GO:0035249 for DP and, to a lesser extent, CF).

Glutamatergic synapses guide the development of growth neurons and regulate feeding motivation in the hippocampus (Huang et al., 2017). The relation between feeding motivation and nutrient intake is crucial to maintaining energy intake and storage (Illius et al., 2002). Such relationship is complex, involving leptin (see above-mentioned gene Siglec-5), and the NPY/AgRP system, which makes food intake-stimulating peptides, which can dramatically influence metabolism and consequently carcass traits (Seabury et al., 2017; Ruud et al., 2020). Among the genes more often represented in the glutamatergic synapse network enriched in our analysis, several were linked with food intake and metabolism (for example, GRM8), eating behavior (GRIK3), insulin secretion, and lipolysis (ADCY1, Olivieri et al., 2016). In support of this hypothesis, we also found out that the enriched KEGG term for DP Glutamatergic synapse (KEGG:04724) belonged to the same group of Circadian entrainments (KEGG:04713) and Apelin signaling pathway (KEGG:04371), both also enriched. Circadian rhythm has a strong connection with feeding behavior (Mrode et al., 2019), and apelin is a peptide connected with food intake and lipid metabolism (Bertrand et al., 2015). The same was true also for CF, with KEGG term Hippo signaling pathway (KEGG:04390) appearing multiple times (Supplementary Figure S4H). This might reflect a greater role of regulatory systems of feeding motivation, nutrient intake, and storage in shaping the variability of these traits. On the other hand, glutamatergic synapses are also involved in physiological responses to stressors and environmental changes. QTLs from the QTLdb associated to our candidate regions for these two traits are reported in Supplementary Tables S2C,D.

The results of our study can help frame the genetic architecture of our between-traits correlation, including such traits that are measures of the same trait in different time points or intervals (the three BW and the three ADG). Within BW, we demonstrated how also from the genomic point of view the weight at the half of the PT was underlined by a mixture of QTLs that were also found either at the start or at the end of the PT. On the other hand, no common SNPs resulted significant both for BW_i and BW_f, and the number of enriched pathways in common was very low (Supplementary Figures S4A–C; Table 3). For what concerns ADG, there was also a deep difference between the signals found for ADG_i and ADG_f, with the latter reflecting much more closely the total ADG, and again no SNPs were shared by ADG_i and ADG_f (Table 3). Moreover, the lowest number of significant SNPs and pathways for BW was at BW_m, and for ADG was ADG_i, with these two traits sharing a temporal correspondence.

Interestingly, we found many genes in common between measures of different traits taken at the same time. For example, both SNPs on BTA7 and BTA1 were significant both for BW_m and ADG_i. Also, one SNP on BTA26 was suggestively associated both for BW_f and ADG_f (Table 3).

These results have several implications: firstly, from an economic point of view, they show that the timing of the trait measurement is crucial. Different life stages can result in different genetic signals; if used for a selection program, this can have an economic and conservation impact. While this is of course expected, given the succession of different biological processes during development, very few studies include such a time stratification in their analysis of productive traits. Even if such a process is difficult to infer, our results show that complexity—intended as the number of functional elements, their diversity, and pathways involved—might increase with age.