Integrating large-scale meta-GWAS and PigGTEx resources to decipher the genetic basis of 232 complex traits in pigs

Abstract

Understanding the molecular and cellular mechanisms underlying complex traits in pigs is crucial for enhancing genetic gain via artificial selection and utilizing pigs as models for human disease and biology. Here, we conducted comprehensive genome-wide association studies (GWAS) followed by a cross-breed meta-analysis for 232 complex traits and a within-breed meta-analysis for 12 traits, using 28.3 million imputed sequence variants in 70 328 animals across 14 pig breeds. We identified 6878 quantitative trait loci (QTL) for 139 complex traits. Leveraging the Pig Genotype-Tissue Expression resource, we systematically investigated the biological context and regulatory mechanisms behind these trait-QTLs, ultimately prioritizing 14 829 variant-gene-tissue-trait regulatory circuits. For instance, rs344053754 regulates UGT2B31 expression in the liver and intestines, potentially by modulating enhancer activity, ultimately influencing litter weight at weaning in pigs. Furthermore, we observed conservation of certain genetic and regulatory mechanisms underlying complex traits between humans and pigs. Overall, our cross-breed meta-GWAS in pigs provides invaluable resources and novel insights into the genetic regulatory and evolutionary mechanisms of complex traits in mammals.

Keywords: PigGTEx; complex traits; meta-GWAS; molQTL; pig.

Figure 1.

Overall study design and summary of genotypes and phenotypes. (a) Overview of study design. WGS: whole genome sequence. GWAS: genome-wide association study. TWAS: transcriptome-wide association study. SMR: summary data-based Mendelian randomization. (b and c) Principal component (PC) analysis of Pig Genomics Reference Panel (PGRP) (b) and GWAS populations (c), which were conducted based on 57 600 individuals (samples with raw genotype data) and a total of 1603 shared array SNPs using PLINK (v1.90) [84]. The first two PCs were plotted using the geom_point function from the ggplot2 (v3.3.6) package in R (v4.1.2). (d) The imputation accuracy (concordance rate between the imputed and observed genotypes) of PGRP in GWAS data. This was evaluated on chromosome 6 with 20 times 5-fold cross-validation using Chip array data from 60 720 pigs in 53 GWAS populations that had individual-level genotype data. Error bars indicate the standard deviation. (e) The total sample size was collected for each trait. Traits were classified into five main trait-categories.

Figure 2.

Summary and validation of quantitative trait loci (QTL) for pig complex traits. (a) The total number of QTL detected in 12 sub trait-categories. The number of traits included in each trait-category is shown in parentheses. (b) A Fuji-plot summarizes the 6878 lead SNPs (P < 5 × 10⁻⁸) identified in 169 within-breed and cross-breed meta-GWAS analyses. The plot was created using the Fuji-plot script developed by Kanai et al. [85]. The innermost ring (ring 1) indicates the number of traits associated with each lead SNP. Rings 2–170 indicate the 169 meta-GWAS. The order of traits is shown in Table S5 (starting with the innermost ring). The points indicate the genomic position of the 6878 lead SNPs in chromosomes 1–18. The colors represent the trait-category, consistent with (a). (c) Pearson correlation between sample size and the number of lead SNPs in 169 meta-GWAS. (d) Pearson correlation between heritability and the number of lead SNPs in 27 meta-GWAS (sample size > 15 000). Heritability was estimated using linkage disequilibrium score regression (LDSC) [86]. The Pearson correlation coefficient in (c and d) was calculated by the cor.test function in R. (e) Results of genomic predictions for individuals from several pig breeds in the PGRP with large phenotype differences, based on a linear mixed model and effect information from suggestive lead variants (P < 1 × 10⁻⁵) and the same number of random non-significant variants in the total number born alive. The x-axis labels indicate the different pig breeds. The y-axis labels indicate the genomic estimated breeding values (GEBVs). The error bars are the standard errors of GEBVs. (f) The number of different categories of QTL detected in individual GWAS and meta-GWAS. Both individual- and meta-GWAS QTL refer to the QTL that exhibit physical position overlap in both individual GWAS and meta-GWAS for the same trait. (g) rs320375241 associated with average daily gain (ADG) in individual GWAS and cross-breed meta-GWAS, respectively. (a), (b) and (c) were three random populations for ADG.

Figure 3.

Identified QTL across breeds. (a) The number of QTL from the within-breed meta-analysis for 12 complex traits in Duroc, Landrace and Yorkshire. Breed-specific QTL refers to the QTL that does not exhibit physical overlap with QTL from other breeds for the same trait. Breed-shared QTL refers to the QTL that exhibits physical overlap with QTL from the same trait in at least one breed. (b) Manhattan plots of within-breed meta-analysis for backfat thickness in Duroc (left), Landrace (middle) and Yorkshire (right). The blue (n = 68), purple (n = 192) and red (n = 211) points represent the breed-specific QTL identified in Duroc, Landrace and Yorkshire, respectively. (c and d) The distribution of allele frequencies (c) and effect sizes (d) of lead SNPs of breed-specific QTL for backfat thickness in the three breeds. Boxplots depict the median value as the center, and first and third quartiles as box boundaries. Sample sizes for each breed are provided in Table S5. (e) The Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) enrichment results of breed-specific highly expressed genes in three breeds identified from muscle RNA samples of PigGTEx [34]. The enrichment was conducted using KOBAS [87]. The sample sizes of muscle for detecting breed-specific highly expressed genes in Duroc, Landrace and Yorkshire were 157, 49 and 119, respectively. The x-axis represents the false discovery rate (FDR)-corrected P-value. (f) Enrichment of breed-specific QTL and non-breed-specific QTL for loin muscle depth in the breed-specific highly expressed genes in muscle. The error bar represents the standard error of the enrichment fold. (g) Local Manhattan plots of the cross-breed meta-analysis of ADG in all breeds (left 1), as well as a within-breed meta-analysis of Yorkshire (left 2), Landrace (left 3) and Duroc (left 4) on chromosome 1. The color indicates the magnitude of the linkage disequilibrium (LD) between rs320375241 and other SNPs. (h) Clustering for a 36 within-breed meta-analysis of 12 complex traits based on the Z-score of breed shared lead SNPs (lead SNPs with m-values > 0.9 in all the three breeds). Clustering was conducted using the umap package in R.

Figure 4.

Exploiting the PigGTEx resource to decipher regulatory mechanisms of GWAS loci. (a) Results of annotation and enrichment of lead SNPs and genome-wide-level significant SNPs from within- and cross-breed meta-GWAS with a sample size greater than 1000 in different categories of genomic regions, conserved elements and regulatory elements. The enrichment analyses were conducted using resampled SNPs that matched the MAF (within 0.02) and LD (within 0.1) of the significant SNPs. The red dots and bars indicate the proportion and the enrichment fold of associated SNPs in the category. Significance of enrichment was indicated by *, ** and *** for P < 0.05, 0.01 and 0.001, respectively. (b) Annotation results of the lead SNPs identified in this study (y-axis) and downloaded from the human GWAS catalog (https://www.ebi.ac.uk/gwas/) (x-axis) in each category of SNPeff and chromatin states (https://figshare.com/articles/dataset/6_type_of_regulator_hg19_zip/13480425). Each dot represents a genomic category. The Pearson correlation coefficient and significance were calculated by the cor.test function in R. (c) The heritability enrichment for 5 types of molecular cis-QTL (molQTL) in 14 complex traits with large sample size. The dashed line represents enrichment fold = 1. The error bar reflects the standard error of the enrichment fold. cis-eQTL: gene expression QTL, cis-sQTL: splicing QTL, cis-eeQTL: exon expression, cis-lncQTL: lncRNA expression QTL, cis-enQTL: enhancer expression QTL. (d) The estimated total SNP heritability contributed by independent molQTL in the liver and muscle for 232 complex traits. Significance is indicated by *** for P < 0.001, which was obtained by the Wolcox.test of ggsignif package in R. (e) The heritability enrichment for the genes of 7 tissue-sharing gradients in 94 complex traits with a sample size > 1000. The red color represents the scaled heritability enrichment fold. The bright blue boxes indicate the position of examples (f and g). The ‘*’ indicates heritability enrichment fold >1 and P < 0.05. Column clusters were produced by the dist function with the ‘euclidean’ method and the hclust function with the ‘complete’ method in R. The heatmap was plotted by ggplot2 package (v3.3.2) in R (v4.2.1). (f) Expression (left, black) and the overall tissue-sharing pattern at LFSR obtained by MashR (v0.2–6) (right, red) of LGALS13 in 34 tissues. (g) Expression (left, black) and the overall tissue-sharing pattern at LFSR obtained by MashR (v0.2–6) (right, red) of ABCC10 in 34 tissues. The red dashed line in (f and g) indicates LFSR < 5%.

Figure 5.

Tissue-specific regulation of GWAS loci. (a) Enrichment results for significantly associated SNPs of 19 cross-breed meta-GWAS in tissue-specific functional regions (the top 1000 tissue-specific highly expressed genes per tissue and their upstream and downstream 100 kb regions) in 34 tissues. Colors indicate enrichment fold. Rows indicate traits and columns indicate tissues. The bright blue box indicates the position of examples (b–g). Enrichment for trait-tissue pairs E_T = p_T (proportion of significant SNPs for trait Tr located in functional regions of tissue Ti)/q_T (proportion of all SNPs located in functional regions of tissue Ti), as calculated with BEDTools v2.25.0 [88]. Associated SNPs were resampled 1000 times with MAF within 0.02 and LD within 0.1, matched to calculate enrichment significance. E_T > 1 and P < 0.05 indicates the significant trait-tissue pair. Significance of enrichment is indicated by *, ** and *** for P < 0.05, 0.01 and 0.001, respectively. (b) Expression of gene UGT2B31 in 34 tissues. (c) A Manhattan plot representing the gene-based GWAS results of litter weight (weaning). (d) A Manhattan plot representing the single-tissue TWAS results of litter weight (weaning) in the liver. (e) The single-tissue TWAS results of UGT2B31 for litter weight (weaning) from S-PrediXcan [89]. (f) The summary data-based Mendelian randomization (SMR) results of litter weight (weaning) GWAS (x-axis) and cis-eQTL (y-axis) of UGT2B31 on chromosome 8 in the liver. (g) The chromatin states around UGT2B31 on chromosome 8.

Figure 6.

Comparison of complex trait genetics between humans and pigs. (a) The heritability enrichment of 126 pig GWAS in 136 human GWAS. The x-axis indicates traits and trait categories of human GWAS. The y-axis indicates the significance of heritability enrichment. The color of the dots indicates the trait category of pig GWAS. The size of the dots indicates the fold of heritability enrichment. The dashed line indicates the significant threshold (P = 0.05/126*136 pig traits). (b) An alluvium-stratum plot showing the correlation between human and pig GWAS. Colors indicate trait categories. (c) Spearman's correlations of traits between humans and pigs, which were estimated by the absolute Z score of homologous variants from GWAS summary statistics. The x-axis indicates the significance of Spearman's correlations. The red dashed line indicates the significant threshold (P = 0.05/136 human traits). The top trait pairs are labeled. (d–f) Similar regulatory mechanisms between body fatness rate (BFR) in humans and ADG in pigs. (d) A local Manhattan plot of GWAS for ADG in pigs (top) and BFR in humans (bottom). The red triangles represent homozygous variants of humans (rs11877146) and pigs (rs322242884). The color indicates the magnitude of LD between the homologous variants and other SNPs. (e) The effects of homozygous variants in (d) on the expression of homozygous gene NPC1 in the muscle of pigs (top) and humans (bottom). (f) The effects of homologous variants in (d) on the expression of the homologous gene TMEM241 in the brain of pigs (top) and humans (bottom). The black error bars in (e and f) represent the standard errors. The significance tests in (e and f) were performed with the wilcox.test function of the ggsignif package in R (v4.2.1). The effect size of eQTL in humans and pigs in (e and f) was derived from the GTEx portal (https://www.gtexportal.org/) and PigGTEx portal (https://piggtex.farmgtex.org/), respectively.