Integration of a single-step genome-wide association study with a multi-tissue transcriptome analysis provides novel insights into the genetic basis of wool and weight traits in sheep

Abstract

Background: Genetic improvement of wool and growth traits is a major goal in the sheep industry, but their underlying genetic architecture remains elusive. To improve our understanding of these mechanisms, we conducted a weighted single-step genome-wide association study (WssGWAS) and then integrated the results with large-scale transcriptome data for five wool traits and one growth trait in Merino sheep: mean fibre diameter (MFD), coefficient of variation of the fibre diameter (CVFD), crimp number (CN), mean staple length (MSL), greasy fleece weight (GFW), and live weight (LW).

Results: Our dataset comprised 7135 individuals with phenotype data, among which 1217 had high-density (HD) genotype data (n = 372,534). The genotypes of 707 of these animals were imputed from the Illumina Ovine single nucleotide polymorphism (SNP) 54 BeadChip to the HD Array. The heritability of these traits ranged from 0.05 (CVFD) to 0.36 (MFD), and between-trait genetic correlations ranged from - 0.44 (CN vs. LW) to 0.77 (GFW vs. LW). By integrating the GWAS signals with RNA-seq data from 500 samples (representing 87 tissue types from 16 animals), we detected tissues that were relevant to each of the six traits, e.g. liver, muscle and the gastrointestinal (GI) tract were the most relevant tissues for LW, and leukocytes and macrophages were the most relevant cells for CN. For the six traits, 54 quantitative trait loci (QTL) were identified covering 81 candidate genes on 21 ovine autosomes. Multiple candidate genes showed strong tissue-specific expression, e.g. BNC1 (associated with MFD) and CHRNB1 (LW) were specifically expressed in skin and muscle, respectively. By conducting phenome-wide association studies (PheWAS) in humans, we found that orthologues of several of these candidate genes were significantly (FDR < 0.05) associated with similar traits in humans, e.g. BNC1 was significantly associated with MFD in sheep and with hair colour in humans, and CHRNB1 was significantly associated with LW in sheep and with body mass index in humans.

Conclusions: Our findings provide novel insights into the biological and genetic mechanisms underlying wool and growth traits, and thus will contribute to the genetic improvement and gene mapping of complex traits in sheep.

Fig. 1

GWAS results of five wool traits and live weight in Merino sheep. a–f show the mean fibre diameter (MFD), coefficient of variation of the fibre diameter (CVFD), crimp number (CN), mean staple length (MSL), greasy fleece weight (GFW), and live weight (LW), respectively. Each dot represents one window region of 20 consecutive SNPs. The X-axis represents 26 autosomes, and the Y-axis represents the percentages of the genetic variance explained by the windows

Fig. 2

Detection of tissues and cell types related to wool traits and live weight in sheep. a Functional enrichment analysis of tissue-specific genes (the top 10% of genes based on t-statistics) for 13 organ systems (central nervous system (CNS), cardiovascular system (Cardio), skin, muscle, liver, lung, kidney, gastrointestinal (GI) tract, endocrine (Endoc), immune system, male reproductive system (Male_R), female reproductive system (Fem_R), and embryonic system). b–d Relationships between the six economic traits and the 13 organ systems, GI tract tissues and immune tissues, respectively. The colour corresponds to the enrichment degree (i.e., − log₁₀FDR), which was computed by sum-based GWAS signal enrichment analysis based on the top 10% tissue-specific genes and a 20-kb extension. *corrected-P (FDR) < 0.05

Fig. 3

Heatmap of 77 of 81 candidate genes based on the sheep expression atlas. The gene expression levels are normalized as transcripts per million (TPM). The colour corresponds to the log₁₀ (TPM + 0.25) value. The Y-axis represents the 77 candidate genes, and the X-axis represents the 13 organ systems

Fig. 4

Heatmap of 78 of 81 candidate genes based on the results of the phenome-wide association study (PheWAS). The colour corresponds to the − log₁₀P-value value from the PheWAS results. The Y-axis represents the 78 candidate genes, and the X-axis represents the 12 trait domains

Fig. 5

Expression patterns and results of the phenome-wide association study (PheWAS) for four candidate genes. a, b BNC1; c, d GHR; e, f CHRNB1; g, h SPHK1. In a, c, e and g, the Y-axis represents gene expression (TPM), and the X-axis represents the samples in 13 organ systems. In b, d, f and h, each dot is one trait. The Y-axis represents the -log₁₀P-value value from the PheWAS results, and the X-axis represents 12 trait domains