Integrating multi-layered biological priors to improve genomic prediction accuracy in beef cattle

Abstract

Background

Integrating multi-layered information can enhance the accuracy of genomic prediction for complex traits. However, the improvement and application of effective strategies for genomic prediction (GP) using multi-omics data remains challenging.

Methods

We generated 11 feature sets for sequencing variants from genomics, transcriptomics, metabolomics, and epigenetics data in beef cattle, then we assessed the contribution of functional variants using genomic restricted maximum likelihood (GREML). We next estimated and ranked variant scores for 43 economically important traits, and compared the prediction accuracy of the top and bottom sets using genomic best linear unbiased prediction (GBLUP) and BayesB model. In addition, we annotated the variants from GWAS with functional feature sets and performed enrichment analysis.

Results

We observed significant enrichments for 32 functional categories in 11 feature sets. The evolutionary related sets (conservation regions and selection signatures) contributed significantly to heritability (31.78-fold and 14.48-fold enrichment), while metabolomics and transcriptomics showed low heritability enrichments. We observed a significant increase in prediction accuracy using the top feature set variants compared to whole-genome sequencing (WGS) data. The prediction accuracy based on the top 10% variant set showed an average increase of 11.6% and 7.54% using BayesB and GBLUP across traits, respectively. Notably, the greatest increase of 31.52% was obtained for spleen weight (SW) using BayesB. Also, we found that the top 10% of variants show strong enrichment with weight related QTLs based on the Cattle QTL database.

Conclusions

Our findings suggest that integrating biological prior information from multiple layers can enhance our understanding of the genetic architecture underlying complex traits and further improve genomic prediction in beef cattle.

Background

Genomic prediction (GP), an effective approach for enhancing selection and promoting breeding efficiency [1, 2], has been widely applied in the fields of plant and animal breeding [3, 4]. Many parametric and nonparametric statistical methods have been proposed to predict Genomic Estimated Breeding Values (GEBVs) [5]. GBLUP constructs a genetic relationship matrix to facilitate the estimation of GEBVs [6]. The Bayesian alphabet assumes a priori that the variances of effects for many single nucleotide polymorphisms (SNPs) are zero, while the effect of SNPs follow Student’s t distribution [7]. Such models formulate distinct hypotheses regarding the distribution of marker effects and their impact on genetic variation [8, 9]. GP has been primarily applied based on SNP arrays [10, 11]. With the cost of whole-genome sequencing (WGS) decreasing, the application of WGS for GP has been widely applied in farm animals [12, 13]. GP using WGS data may promote prediction accuracy because it covers more SNPs across the genome than SNP arrays [14, 15]. A previous simulation study showed the superiority of BayesB over GBLUP using WGS, while the accuracy of GEBV increased when compared to SNP array [16, 17]. GP based on WGS data, was more accurate when using BayesB than using GBLUP [18].

The use of WGS data for GP may be limited by the substantial presence of linkage disequilibrium (LD) and diverse genomic functional regions. This scenario reduces the signal-to-noise ratio when employing WGS data directly for GP without a biological prior [19,20,21], thus many studies have been conducted to incorporate genomic information into statistical models by controlling for LD and annotating variants based on different functional classes [22, 23]. Moreover, the availability of multi-omics information (e.g., genomics, transcriptomics, proteomics and metabolomics) bridges a vital link between genotypes and phenotypes which provides a biological prior for genomic prediction. The Genotype-Tissue Expression (GTEx) Project was initiated with the objective of collating genetic influences on gene expression in human tissues and facilitating an enhanced understanding of the dynamics of regulatory genetic variation by elucidating the molecular mechanisms underlying genetic correlations with complex diseases and traits [24]. Similarly, the Farm Animal Genotype-Tissue Expression (FarmGTEx) Project serves as an extensive public repository, facilitating the discovery of tissue-specific genetic regulatory variants and the prediction of molecular phenotypes in farm animals [25, 26]. Additionally, the Functional Annotation of Animal Genomes (FAANG) Project seeks to improve comprehension of genome functionality through comprehensive annotation efforts. This annotation information helps to refine the accuracy and sensitivity of genomic selection strategies for animal [27, 28].

Many approaches have been developed to estimate genomic values and improve genomic prediction accuracy using multi-omics information [28,29,30]. Ye et al. refined the GFBLUP model using transcriptomics information in Drosophila [31], revealing that significant variations identified through a transcriptome-wide association study (TWAS) contribute more effects than those identified through genome-wide association study (GWAS). Xu et al. applied Bayesian ridge regression (BR) to evaluate the accuracy of GP for complex traits by integrating transcriptomics, proteomics, and metabolomics data, suggesting that large biobanks could reliably and efficiently explore trait–disease associations using multi-omics data [32]. Additionally, Hu et al. designed a novel GP strategy called multilayered least absolute shrinkage and selection operator (MLLASSO) by integrating multi-omics data into a single model, their finding suggested MLLASSO can significantly improve the predictability of yield in rice [33]. Although these approaches have been successfully employed in many studies, the integration of disparate data types into comprehensive system-scale analyses represents a significant challenge [34].

In this study, we generated genomic feature sets for sequencing variants by integrating multi-layered biological priors in beef cattle, then we assessed the contribution of functional variants and estimated variant scores for 43 economically important traits. Further, we evaluated and compared the GP accuracy based on functional variants using GBLUP and BayesB (Fig. 1).

Fig. 1

Schematic overview of current study. a Data collection. We divided the full variants from WGS data into 11 feature sets (annotation, LD, allele frequency, variant density, p-variants, selection signature, conservation, eQTL, mQTL, OCR and HMRs) from genomics, transcriptomics, metabolomics, and epigenetics data. b. Calculate variants score. For each of the 43 traits, we estimated the variance explained by the random effects associated with each GRM using GREML. Each GREML analysis incorporated two random effects: one based on the targeted GRM and another based on the GRM derived from the remaining variants. We calculated the proportion of genetic variance attributed to the targeted GRM for each trait. To determine the per-variant heritability, we divided the explained variance by the number of variants in the set. Finally, we averaged this value across the 11 functional sets for each variant. c. Validation analysis. To assess the reasonableness of the scores, we established six thresholds: “top-5”, “top-10”, “top-30”, “bottom-5”, “bottom-10” and “bottom-30” We then compared the variance explained by each threshold with the accuracy of the genomic predictions. To ascertain whether there are pertinent QTL enrichments for our top variants, we conducted a QTL enrichment analysis using the Cattle QTL Database

Methods

Dataset

The measurement of phenotypes and genotypes was performed as described in our previous studies [35, 36]. The phenotypic data were generated from 1577 Huaxi (derived from Chinese Simmental beef cattle), which were born between 2008 and 2020 from Ulgai, Xilingol League, and Inner Mongolia, China. After weaning, all individuals were moved to Jinweifuren Co., Ltd. for fattening under the same feeding and management conditions. All samples were genotyped by Illumina BovineHD SNP array. The SNP positions were determined based on the ARS-UCD1.2 reference genome, the SNP imputation was carried out using Run 8 of the 1000 Bull Genomes Project and 44 representative individuals from our studied population. After filtering variants with the threshold of MAF < 0.05 and DR² < 0.8, we retained a total of 10,213,925 autosome SNPs with an average DR² of 0.93. A total of 43 traits including carcass and beef cut traits were included in this study. Detailed information for traits is presented in (Supplementary file 1, Table S1). All of the traits were adjusted by gender, year, acid remove day, enter weight and enter day. The phenotypes were adjusted by the significant factors using the glm function in R.

Annotation

The SnpEff software was used to annotate and predict the effects of genetic variants, and the bovine genome annotation (ARS-UCD1.2) was downloaded from Ensembl [37, 38]. According to genome annotation information, the bovine genome was partitioned into six genomic classes, including 1) intergenic.regions, 2) intronic.regions, 3) geneend.regions, 4) cording.related.regions, 5) regulatory.related.regions, 6)3’ untranslated regions (UTR) and 5’ UTR.

LD, allele frequency and variant density

The levels of LD, allele frequency and variant density were divided based on three quartiles. We used the GCTA software to calculate the LD score in the surrounding 50 kb region, and then we used the LD score of each variant to bin all variants by quartile [39, 40]. Variants were unevenly distributed across the genome. VCFtools software was used to calculate the density within fixed 50 kb windows [41]. The allele frequency was divided by the minor allele frequency (MAF). The values of the quartiles were as follows: fourth quartile > third quartile > second quartile > first quartile.

Potential variants for production traits using WGS

We retrieved candidate variants related to body size and beef production traits in cattle, which are publicly available from sequence-based meta-GWAS with a larger number of animals from diverse cattle populations [42, 43]. Finally, 583,438 variants were used to define as the potential variants (p-variants).

Selection signature

The selscan software with a setting of the max-gap 800,000 bp was used to estimate the iHS score for autosomal SNPs [44]. The norm module of selscan was applied to normalize the iHS score, and single site values for iHS were averaged in nonoverlapping windows of 50 kb across the genome. Regions in the top 1% with the highest average |iHS| score and SNP numbers greater than 10 were regarded as candidate regions under positive selection.

Conserved sites

In this study, we used conserved sites to map gene regions that may be involved in basal metabolism. The conserved genome sites in cattle were transformed from humans. First, the conserved sites were based on conservation between 100 vertebrate species, and the Wiggle file was downloaded from UCSC. All the conserved sites in the human genome were 113,280,297, the sites were lifted over to the cattle genome (102,953,048) by LiftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver), and only the PhastCon score > 0.9 was chosen. Finally, a total of 192,825 variations were remained after merging the converted variants.

Meta-analysis of expression QTLs (eQTLs)

We conducted cis-eQTL mapping for three tissues (muscle, liver and adipose) according to our previous study [45]. The SNPs located within 1 Mb up/downstream of the transcription start sites (TSSs) were defined as potential cis-eQTLs [46]. In this study, each variant had an estimate of the effect and standard error (se), allowing for the three tissues to perform meta-analysis in this study by METAL software [47]. We obtained 240,683 variants as the eQTLs sites under threshold of 0.05 based on the false discovery rate(FDR) .

Meta-analysis of metabolic QTLs (mQTLs)

A total of 397 metabolites were extracted from 117 individuals. The polar metabolome extracts were analyzed using reversed-phase chromatographic separation with positive and negative ionization detection. The metabolome data were measured and corrected according to our previous study [48]. The metabolites were screened using three criteria: 1) We computed the pearson correlation coefficient between the traits and the metabolites. A total of 328 metabolites correlated with at least one trait were retained for subsequent analysis ($|r|$≥ 0.25). 2) The heritability of relevant metabolites was calculated, then we obtained 74 significant heritable metabolites (0.1< ${h}^{2}$  < 0.9 and $P$<0.05). 3) GWAS analyses were performed based on the LMM for the 74 metabolites. We considered only the GWAS of 66 metabolites with inflation factors ranging from 0.98 to 1.02 for the meta-analysis in METAL software [47]. We obtained 386,995 candidate pleiotropic variants for metabolites with < $P$9.79E-08.

Epigenetic signals

Peak calling for muscle samples was performed using Genrich with the following parameters: -m 30, -j (ATAC-seq mode), -r (remove PCR duplicates), -e MT (to exclude mitochondrial chromosome), -q 0.05 (FDR-adjusted P-value). We obtained 495,903 variants as open chromatin region (OCR). For hypomethylated regions (HMRs) detection, we chose a 10 kb window size for muscle samples using the Methpipe software with the default parameters [49]. We obtained 762,835 variants as HMRs.

Variation score construction

The GREMLwas used to estimate the variance components. First, the different genomic relationship matrixes (GRMs) were made of the target variants and the remaining variants. Target variants were identified in 11 categories (conservation, annotation, selection signature, p-variants, LD, allele frequency, density, eQTLs, mQTLs, OCRs, and HMRs). Variants in these 11 categories refer to the target variants, whereas the remaining variants refer to the non-target variants. The variance components are estimated using a linear mixed model.

$$y^{*} = \mathop \sum \limits_{c = 1}^{C} g_{c} + e$$

(1)

where ${\mathbf{y}}^{\mathbf{*}}$ is the vector of adjusted phenotypic values, ${\mathbf{g}}_{{\varvec{c}}}$ is the vector of individual polygenic effects associated with annotation group${\varvec{c}}$, $C$ is the total number of fitted annotation groups, and $\mathbf{e}$ is the random residual effect, which is assumed to follow a normal distribution of $\mathbf{e}$ ~$N\left( {0,\sigma_{e}^{2} I} \right)$. where ${\mathbf{g}}_{\mathbf{c}}$ is the genomic relationship matrix (GRM) computed using the variants present in category, which were calculated by Yang’s method [50]. Then, the GCTA software was used to calculate the variance explained by random effects described for each GRM [39]. The mean heritability of the variation within target region for each trait was calculated, the partitioned heritability were estimated via the mean heritability divided by the number of variants [51]. To avoid LD heterogeneity along the genome, the LD level was chosen to adjust the score [52, 53] (Supplementary file 1, Table S2).

$${\text{S}}_{{{\text{trait}}}} = \frac{{\mathop \sum \nolimits_{1}^{43} h_{i}^{2} }}{43}$$

$${\text{S}}_{{{\text{set}}}} = \frac{{{\text{S}}_{{{\text{trait}}}} }}{N}$$

$$\overline{{h_{LD}^{2} }} = \frac{{n_{LD1} \times \overline{{h_{LD1}^{2} }} + n_{LD2} \times \overline{{h_{LD2}^{2} }} + n_{LD3} \times \overline{{h_{LD3}^{2} }} + n_{LD4} \times \overline{{h_{LD4}^{2} }} }}{{n_{LD1} + n_{LD2} + n_{LD3} + n_{LD4} }}$$

(2)

$$\overline{{{\text{S}}_{{{\text{adj}}}} }} = {\text{S}}_{{{\text{set}}}} - \overline{{h_{LD}^{2} }}$$

$$S_{{{\text{variant}}}} = \frac{{\sum\nolimits_{1}^{{11}} {\overline{{S_{{{\text{adj}}}} }} } }}{{11}}~$$

${\text{S}}_{{{\text{trait}}}}$ is the heritability after pooling an average of 43 traits, and ${S}_{set}$ is the average heritability of the variant within each category, $N$ is the number of the functional sets. For the LD feature set, ${n}_{LD1}$, ${n}_{LD2}$, ${n}_{LD3}$ and ${n}_{LD4}$ are the numbers of variant members within the 1st, 2nd, 3rd and 4th LD score levels, respectively. $\overline{{h }_{LD1}^{2}}$, $\overline{{h }_{LD2}^{2}}$, $\overline{{h }_{LD3}^{2}}$ and $\overline{{h }_{LD4}^{2}}$ are the mean heritability of the 43 traits at the LD level. $\overline{{S }_{adj}}$ is the value after adjustment by LD. ${S}_{variant}$ is the mean heritability for 11 categories and is the combined score value of variance [51].

The enrichment of each category is quantified using the ratio of $EST$ to $EXP$. $EST$ denotes the estimated total heritability associated with the category, normalized by the estimated SNP heritability, which represents the proportion of heritability attributable to SNPs within that specific category. In contrast, $EXP$ signifies the expected contribution of the category to overall SNP heritability, which was calculated by the ratio of the number of SNPs within the category to the total number of SNPs analyzed [54].

Finally, the combined variant score among the 11 feature sets was obtained. The top 5%, 10% and 30% and bottom 5%, 10% and 30% ranked  variants were selected as the “top-5”, “top-10”, “top-30”, “bottom-5”, “bottom-10” and “bottom-30”, respectively.

GBLUP model

$$y^{*} = Xg + e$$

(3)

where ${\mathbf{y}}^{\mathbf{*}}$ is the vector of adjusted phenotypic values, $\mathbf{X}$ is the design matrix selected by different thresholds (the top 5%, 10%, and 30% and bottom 5%, 10%, and 30%) relating additive genetic values to the phenotype, $\mathbf{g}$ is the vector of genomic values captured by the genetic markers linked to target variants, which follow a normal distribution $g \sim {\mathbf{N}}\left( {0,\sigma_{g}^{2} {\mathbf{G}}} \right)$. $\mathbf{G}$ was calculated using Yang’s method, and $\mathbf{e}$ is a vector of random residual effects. The GBLUP model is implemented in the GCTA software [39].

BayesB model

In this study, we chose the BayesB model to predict the individual phenotypes, the equation:

$$y^{*} = \mathop \sum \limits_{j}^{k} z_{j} a_{j} \delta_{j} + e$$

(4)

where ${y}^{*}$ is the adjusted phenotype, ${z}_{j}$ is the vector of genotype (0,1,2) across animals for the SNP, ${a}_{j}$ is the allele substitution effect for the SNP, and ${\delta }_{j}$ is an indicator of whether the SNP was included (${\delta }_{j}$ = 1) or excluded (${\delta }_{j}$ = 0) in the model for a given Markov chain Monte Carlo (MCMC) iteration. The BayesB model is implemented in the GCTB software [55].

Evaluation of prediction performance

To avoid the impact of the number of variants on the accuracy of genomic prediction, we performed subsampling on WGS data by simulating a variant set containing 10% variants. The autosomal variants were sorted based on their base pair (BP) positions and then divided into bins, each containing 10 variants. From each bin, one variant was randomly selected to compose the variant panel [56, 57]. We then performed genomic prediction with simulated panels using the GBLUP model.

The accuracy of the predictions was assessed using a five-fold cross-validation method with five repetitions. Genomic prediction accuracy ($Acc=cor\left({y}^{*},GEBV\right)$) was determined by calculating the Pearson correlation coefficient between adjusted phenotypic values and GEBVs separately for each of the five-fold cross-validation replicates.

The GWAS significant regions mapping and QTLs enrichment

The GWAS was performed by using a mixed linear model-based association analysis in GCTA software [39]. The mixed linear model was used for GWAS:

$$y^{*} = bx + g + e$$

(5)

where ${\mathbf{y}}^{\mathbf{*}}$ is the adjusted phenotype, $\mathbf{b}$ is the additive effect (fixed effect) of the candidate SNP to be tested for association, $x$ is the SNP genotype indicator variable coded as 0, 1 or 2, $\mathbf{g}$ is the random effect and accumulated effect of all SNPs, $\mathbf{g}$~$N\left(0,{\sigma }^{2}\mathbf{G}\right)$, and $\mathbf{e}$ is a vector of random residual effects. FDR was used to determine the threshold values for single-trait GWAS.

QTL enrichment analysis was carried out using the GALLO package [58] by comparing the number of annotated QTLs within candidate regions to the total number of QTLs in the Cattle QTL database [59].

Results

Summary statistics and genetic parameter estimations of 43 traits

The detailed summary statistics of the 43 traits are presented in (Supplementary file 1, Table S1), including the mean, median, coefficient of variation (CV), and standard deviation (SD). As expected, strong phenotypic correlations were observed between traits (Fig. 2a). The heritability estimates from the GREML are displayed in Supplementary file 1, Table S1. Seventeen traits showed high heritability, 16 showed medium heritability, 10 showed low heritability. Hind shin (HS) had the highest heritability (0.66), while fat cover percentage (FCR) showed the lowest heritability (0.05) (Fig. 2b).

Fig. 2

The phenotype and genetic parameters of carcass traits. a The phenotype correlation between the 43 carcass traits is presented in the form of a color-coded box plot. The color of each box represents a positive correlation (red) or a negative correlation (blue). b The heritability of the 43 traits is presented in a similar format. Low heritability (0~0.2), medium heritability (0.2~.4), and high heritability (0.4~1)

The category of variants

In this study, we divided the full variants from WGS data into 11 feature sets (annotation, LD, allele frequency, variant density, p-variants, selection signature, conversation, eQTL, mQTL, OCR and HMRs). First, we annotated a total of 10,213,925 variants based on the ARS-UCD1.2 reference genome (Table 1). Among them, the majority of variants were located in intergenic and intronic regions, while only 0.65% were found in coding regions. The p-variant set contains 583,438 variants.

Table 1 Merged and original functional annotation of sequence variants

Selection signature set were selected based on the top 1% of regions with the highest |iHS| values (109,325). Different sets of variants were also divided according to the distribution of LD score (LD1, LD2, LD3, and LD4), MAF (Freq1, Freq2, Freq3, and Freq4) and variant density (Density1, Density2, Density3, and Density4) based on the quartile approach. In addition, conserved sites were selected according to a PhastCon score > 0.9 (Table 2). Finally, a total of 192,825 variants remained in our subsequent analysis.

Table 2 The summary of the functional annotation sets in this study

For the transcriptomic data, a meta-analysis was subsequently carried out to identify the candidate variants that influence the expression of genes. We obtained 240,683 variants as eQTL sites with a threshold (FDR < 0.05) (Supplementary file 2, Fig. S1a).

For the metabolomic analysis, metabolites with |Correlation coefficients| > 0.25 for at least one trait were retained. Heritability analysis was then performed on the remaining 328 metabolites with strict criteria ($0.1< h2 < 0.9$ and $P< 0.05$) (Supplementary file 2, Fig. S1b). Finally, 254 categories remained for subsequent analysis. GWAS was performed based on 74 metabolites, and we obtained 66 metabolites for meta-analysis after filtering based on the inflation factor. We identified 386,995 mQTL variants with a threshold of 0.05/10,213,925 (Supplementary file 2, Fig. S1c). For epigenomics data, a total of 495,903 OCRs were detected in muscle tissue by Genrich. We identified 762,835 HMRs in muscle tissue using Methpipe with a 10 kb window size.

The partitioned heritability estimation based on feature sets

We calculated the partitioned heritability for 43 traits to assess the contribution of each feature set to traits (Table 3). Regulatory-related regions accounted for 6.22% of ${S}_{trait}$ while representing only 0.11% of genome variants. Our analysis showed a decrease in the square of the average heritability ${S}_{trait}$ from the first to fourth quartile, indicating local LD and variant density influence variants effect. P-variants explained approximately 9.35% of the average genetic variance and constituted 5.71% of WGS variants. Conservation and selection signature variants accounted for 18.87% and 4.64% of genetic variance, respectively, despite representing only 1.89% and 1.07% of genome variants. Variants derived from transcriptomics and metabolomics data explained 2.74% and 5.84% of genetic variance. Epigenetically, HMRs and OCRs explained 7.43% and 9.54% of genetic variance. Overall, functional annotation sets contributed significantly to heritability (Fig. 3a), with regulatory-related SNPs (11,534) showing the greatest enrichment (172.96-fold), and coding regions providing a per-SNP predictability enrichment of 9.40. Smaller SNP counts in conservation and selection signatures yielded substantial contributions (31.78-fold and 14.48-fold enrichment). Metabolomics data provided greater heritability enrichment (5.03-fold) than transcriptomics (3.78-fold), while in the epigenetics category, OCRs outperformed HMRs with a 6.37-fold enrichment versus 3.22-fold. After correcting for LD categories, the ranking of variant sets from per-variant heritability showed highly correlated (R² = 0.91) with that of unranking variant sets(Fig. 3b).

Table 3 The partitioned heritability enrichment for the functional annotation sets

Fig. 3

The ranking analysis of 11 feature sets. a The heritability enrichment ratio of functional variant sets averaged across beef cattle. SLSN represents selection signature b The correlation between the ranking of variant set based on the LD adjusted per-variant $\overline{{S }_{adj}}$ (y-axis) and the ranking of variant set based on the unadjusted (observed) per-variant ${S}_{set}$ (x-axis)

Performance of GP based on different sets in the GBLUP model

We examined the heritability and compared the prediction accuracy of different variant sets and WGS data using GBLUP (Supplementary file 1, Table S3 and S4). The top sets performed significantly better than the bottom sets. We found that the top sets (averages of 0.3740, 0.4026, and 0.4267 for top-5, top-10, and top-30, respectively) had significantly higher heritability estimates than the bottom sets (averages of 0.2420, 0.2578, and 0.3164 for bottom-5, bottom-10, and bottom-30, respectively) (Fig. 4a), and the top sets (averages of 0.2793, 0.2865, and 0.2841 for top-5, top-10, and top-30, respectively) achieved higher accuracies than the bottom sets (averages of 0.2376, 0.2364, and 0.2487 for bottom-5, bottom-10, and bottom-30, respectively).

Fig. 4

Genomic prediction based on different thresholds using the GBLUP model a The heritability among 6 thresholds and WGS data in the GBLUP model b The accuracy of genomic prediction among 6 thresholds and WGS data in GBLUP model

Across all sets and WGS data, the heritability of the top-30 exceeded that of GBLUP, with an average increase of 6.54%. Specifically, SW, knuckle (KK), and bone weight (BW) increased by 13.76%, 10.01%, and 9.36%, respectively. For the top-10 set, the LRW trait exhibited 20.67% higher heritability compared to WGS data. Regarding prediction accuracy, the top-10 set significantly outperformed the WGS data, with an average improvement of 7.54%. Notably, the LRW, SW, and hind leg length (HLL) accuracies increased by 29.40%, 11.30%, and 7.29%, respectively (Fig. 4b).

In this study, we randomly selected the 10% variants from the WGS data and compared with the full WGS data, we observed enhancements in prediction accuracy (4.18%, 2.67%, and 0.78%) for the GBLUP model for average daily gain (ADG), carcass weight (CW), and live weight (LW), respectively. The top-10 set demonstrated higher accuracy (4.18%, 2.76%, and 7.95%) compared to the randomly selected variants in ADG, LW, and KK (Supplementary file 2, Fig. S2).

Performance of GP based on different sets using the BayesB model

In this study, we found that the top-5 set achieved higher accuracy compared to the bottom-5 (average ~22.43%) using BayesB, the top-10 set outperformed the bottom-10 set (average ~24.75%), and the top-30 set surpassed the bottom-30 set (average ~ 14.50%). These findings are consistent with the performance observed in the GBLUP model (Fig. 5a). Furthermore, our analysis indicated that the top-10 demonstrated superior predictive accuracy when utilizing both the GBLUP and BayesB models. Specifically, we observed that the top-10 outperformed the top-5 by an average of approximately 2.88% and surpassed the top-30 by an average of approximately 2.83% (Supplementary file 1, Table S3 and S4).

Fig. 5

Genomic prediction based on different thresholds using the BayesB model a The accuracy of genomic prediction among 6 thresholds in BayesB model. b The heritability among 6 thresholds in the BayesB model

Also, the heritability estimates of the studied traits using the BayesB were lower than GBLUP (Fig. 5b, Supplementary file 2, Fig. S3a), and BayesB significantly outperformed GBLUP in terms of genomic prediction accuracy. We found a significant increase in the prediction accuracy using BayesB. For instance, the prediction accuracy for the SW trait using the top-30 set showed an increase of 34.48% over the WGS data, and it increased by 31.52% for the SW using the top-10 set. Compared with the other methods, BayesB exhibited more improvements based on the top-10 set than others. For example, LRW, BW, and KK achieved the highest improvements in accuracy (24.90%, 12.95%, and 9.31%, respectively). Compared to GBLUP, BayesB showed higher predicted accuracy based on the top feature sets for all traits except LRW (Supplementary file 2, Fig. S3b). Particularly, the SW showed an approximately 18.17% increase based on the top-10 set compared to GBLUP, while BayesB model improved by 3.85% over GBLUP on average. The unbiasedness of genomic predictions (from ~0.84 to ~1.17) based on different thresholds were also estimated using the GBLUP and BayesB models (Supplementary file 2, Fig. S4).

Overlap of GWAS variants with functional features and QTLs enrichment analysis

The prioritized variants in the variant score were mapped to significant regions from the GWAS results for the HS and LRW. The top variant with an FDR of 3.3E-04 in the candidate region (BTA 6: 36,816,554-37,883,636 bp) was associated with HS, exhibiting a strong LD with nearby variants. Subsequently, we annotated the variants in this region with functional feature sets and observed that several variants overlapped with selection signatures, conservation, p-variants, and OCRs (Fig. 6a). For LRW, a QTL (BTA 14:25,906,554-26,883,636) was annotated with variants from OCRs (51.65%), mQTLs (25.87%), HMRs (10.39%) and conservation regions (3.87%) (Fig. 6b).

Fig. 6

Integrative GWAS analysis of two traits. a Results of the GWAS of the HS trait; the region plot of BTA 6. The colors of the variants are based on their LD with the most significant variants. b GWAS results for the FZ trait; the region plot of BTA 14. The colors of the variants are based on their LD with the most significant variants

We carried out the overlap analysis between the top-10 variant set and the functional annotation set, revealing significant correlations between the occurrence ratio of the functional annotation set and their enrichment folds (R² = 0.56). Variants within the functional annotation set that exhibited high enrichment were more frequently represented in the top-10 set. Notably, all variants categorized as regulatory-related, conservation, selection signatures, and coding-related were identified within the top-10 set. Furthermore, our findings indicated that variants within the lower enrichment-fold functional annotation set, specifically those in LD4, Freq4, and Density4 categories, occurred less frequently compared to other sets (Supplementary file 1, Table S6). Based on the Cattle QTL database [59], the top-10 set of variants was predominantly enriched with the production, meat, and carcass categories of QTLs, representing 19.15% and 21.96% of these categories, respectively (Supplementary file 2, Fig. S5a). Further enrichment analysis revealed a significant association of the top-10 set with weight related QTLs, such as metabolic body weight, average daily gain, carcass weight, and longissimus muscle area (Supplementary 2, Fig. S5b). Notably, these QTLs showed a prominent enrichment in meat color and conformation traits, demonstrating the highest richness factor among the analyzed traits.

Discussion

In recent years, numerous GP approaches have emerged for farm animals and plants [60, 61]. To improve the prediction accuracy, many methods have been proposed by integrating biological priors from multi-layered information [62, 63]. In cattle, a previous study performed GP using the BayesRC model by integrating independent biological priors, their findings revealed that BayesRC was more effective than BayesR in identifying candidate causal variants and predicting milk traits[64]. Liang et al. utilized transcriptomic data as a T matrix and combined them with wmssGBLUP and reported that the transcriptome data has the potential to improve genomic predictions [65]. However, a comprehensive methodology for GP by integrating this multi-layered information has not yet been fully developed [66]. To address this issue, we constructed 11 feature sets from multi-omics data and evaluated the contributions of functional variants to 43 economically important traits in beef cattle.

Multi-omics analyses are useful to characterize the regulatory regions and annotated mammalian genomes [27, 67,68,69,70]. Our study integrated multi-omics data, including comprehensive evolutionary, selection signal, transcriptomic, metabolomic and epigenetic data, based on the FAETH framework [51]. By estimating the variance explained by each feature set using GREML, we selected functional variant sets based on their ranking scores. This strategy can help to improve the accuracy of genomic prediction, and understand the genetic architecture underlying complex traits.

Our findings underscore the significance of incorporating functional annotation into genomic analyses [71, 72]. Zeng et al. revealed that genomic loci displaying conservation across a wider array of species were more prone to containing variants correlated with a heightened enrichment of heritability [73]. In humans, conserved regions of the genome have been found to markedly enhance the estimation of trait heritability [74], a phenomenon often attributed to the concentration of functional elements within these conserved regions [75]. Our findings revealed that the selection signature set was capable of elucidating a significant portion of the genetic variation observed in studied traits, which can contribute larger effect for GP compared with WGS data as reported by Kemper et al. [76]. Furthermore, our analysis of the set ranking result confirmed that the LD heterogeneity of variants have a substantial impact on trait heritability, which is consistent with recent evidence for the strong influence of LD properties on complex traits [19, 53].

Several studies had showed that the incorporation of transcriptomic and metabolomic information can help to elucidate the genetic basis for complex traits[77,78,79]. Intermediate QTLs such as eQTLs and mQTLs have consistently demonstrated their significance in contributing to the regulation of complex traits [80, 81]. These intermediate QTLs act as crucial mediators, connecting genetic variants to phenotypic outcomes [82, 83]. In this study, we found that the prediction accuracy using metabolomics information surpasses that of transcriptomic information, this could be explained by that the complex gene expression patterns of the diverse tissue sampled in different developmental stages [84,85,86]. Our analysis further revealed the pivotal role of intermediate QTLs, emphasizing their importance in the genetic architecture of complex traits.

Different models may influence the predictive accuracy of GEBVs [10, 18, 87]. Using the GBLUP model, Xiang et al. predicted GEBVs by calculating the variant scores in dairy cattle, and their findings suggested that genomic prediction using high-ranking variants was more accurate than genomic prediction using low-ranking variants in most scenarios [51]. In this study, we found that the top variant sets showed higher prediction accuracy than the bottom sets using the BayesB and GBLUP, which is consistent with the findings of a previous study by Xiang et al. [51].The difference between GBLUP and BayesB was mainly caused by the different assumptions regarding variation effects. Compared with the GBLUP model, the performance of Bayesian models is superior when the studied trait are controlled by multiple QTLs [88, 89].

Our study suggested that pinpointing the functional variants on the top feature set may contribute larger effect for the genetic architecture underlying complex traits. Our approach provided a comprehensive framework for GP from multi-layered biological priors, and refined our understanding of complex trait genetics architecture. New approaches from machine learning and deep learning should advance analysis strategies for incorporating multi-layered biological datasets and promote genetic gains in animal breeding programs [66, 90, 91]. Overall, integrating multi-omics data from new approaches can further facilitate investigations of the functional impacts of variants and improve genomic prediction accuracy in farm animals.

Conclusions

Our study revealed that pinpointing the effect of variants on the top feature set can enhance  our understanding of the genetic architecture underlying complex traits. Genomic selection by integrating multi-layered biological priors can improve prediction accuracy for important traits in cattle.