<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /><figcaption>Extended Data Fig. 5 | Conditionally independent molecular QTLs (molQTL).</figcaption></figure>

a, Distribution and average number (red dots, right y-axis) of conditionally independent cis-QTL per eMolecules across 34 tissues. Tissues (x-axis) are ordered by increasing sample size. b, Cumulative proportion of distance to the transcription start site (TSS) of target genes for conditionally independent

cis-eQTL in each of 34 tissues. The meanings of the colors of curved lines are the same as the color key in panel (a). c,d, Comparison of distance to TSS (c) and effect size (|log₂(aFC)|) (d) among top three independent cis-eQTL per eGene across 34 tissues. The aFC is for allelic fold change. The *P* values were obtained by the two-sided Wilcoxon rank-sum test.
<figure><img src="image_42.png" /><figcaption>Extended Data Fig. 6 | See next page for caption.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Extended Data Fig. 6 | See next page for caption.</figcaption></figure>
**Extended Data Fig. 6** | Validation of cis-eQTL. **a**, Pearson's correlation of combined summary statistics (for example, Z-score, slope and *P*-value (–log₁₀ scale)) of cis-eQTL for all the eGenes across 34 tissues between TensorQTL (linear model, LM) and fastGWA (mixed linear model, MLM). **b**, Pearson's correlation of summary statistics for each eGene in each tissue between LM and MLM. **c**, Distribution of the Pearson's correlations of Z-score between LM and MLM. **d**, Relationship between correlations of Z-score and the number of significant eQTL across all the eGenes. **e**, Correlation of *P* values derived from MLM and nominal (left) or permutation-corrected (right) *P* derived from LM for the lead eQTL of all the eGenes. **f**, Replication rates (πᵢ) of blood cis-eQTL between the PigGTEx discovery population (n = 386, Discovery) and the external datasets (n = 179). For πᵢ calculation, rows are discovery populations, and columns are replication populations. The external datasets include whole-blood-cell RNA-Seq data and SNP Chip array (Chip) from 179 animals at two developmental

stages (T1 and T2). The prefix 'RNA' and 'Chip' indicate imputed genotypes from RNA-Seq and SNP array, respectively. **g**, Spearman's correlation (ρ) of effect size (z-scores) for blood cis-eQTL among the same populations above. **h**, Replication rates (πᵢ) of PigGTEx cis-eQTL in external validation datasets of three tissues, including muscle (n<sub>PigGTEx</sub> = 1,321, n<sub>external</sub> = 100), liver (n<sub>PigGTEx</sub> = 501, n<sub>external</sub> = 100) and duodenum (n<sub>PigGTEx</sub> = 49, n<sub>external</sub> = 100). The x-axis is the nominal *P*-value of cis-eQTL detected from dataset₂ and is significant in dataset₁ (that is, dataset₁ in dataset₂). **i, j**, Spearman's correlation (ρ) of effect sizes (allelic fold change, aFC in log₂ scale) between cis-eQTL and matched allele-specific expression (ASE) loci in the liver (i) and brain (j). N indicates number of tested loci. The lines are fitted by a linear regression model using the *geom.smooth* function from *ggplot2* (v3.3.2) in R (v4.0.2). The shading represents the standard error of the fitting line. **k**, Spearman's correlation (ρ) of effect sizes between cis-eQTL and matched ASE loci across 34 tissues. Red dots indicate number of tested loci (right y-axis).
**Extended Data Fig. 6** | Validation of cis-eQTL. **a**, Pearson's correlation of combined summary statistics (for example, Z-score, slope and *P*-value (–log₁₀ scale)) of cis-eQTL for all the eGenes across 34 tissues between TensorQTL (linear model, LM) and fastGWA (mixed linear model, MLM). **b**, Pearson's correlation of summary statistics for each eGene in each tissue between LM and MLM. **c**, Distribution of the Pearson's correlations of Z-score between LM and MLM. **d**, Relationship between correlations of Z-score and the number of significant eQTL across all the eGenes. **e**, Correlation of *P* values derived from MLM and nominal (left) or permutation-corrected (right) *P* derived from LM for the lead eQTL of all the eGenes. **f**, Replication rates (πᵢ) of blood cis-eQTL between the PigGTEx discovery population (n = 386, Discovery) and the external datasets (n = 179). For πᵢ calculation, rows are discovery populations, and columns are replication populations. The external datasets include whole-blood-cell RNA-Seq data and SNP Chip array (Chip) from 179 animals at two developmental

stages (T1 and T2). The prefix 'RNA' and 'Chip' indicate imputed genotypes from RNA-Seq and SNP array, respectively. **g**, Spearman's correlation (ρ) of effect size (z-scores) for blood cis-eQTL among the same populations above. **h**, Replication rates (πᵢ) of PigGTEx cis-eQTL in external validation datasets of three tissues, including muscle (n<sub>PigGTEx</sub> = 1,321, n<sub>external</sub> = 100), liver (n<sub>PigGTEx</sub> = 501, n<sub>external</sub> = 100) and duodenum (n<sub>PigGTEx</sub> = 49, n<sub>external</sub> = 100). The x-axis is the nominal *P*-value of cis-eQTL detected from dataset₂ and is significant in dataset₁ (that is, dataset₁ in dataset₂). **i, j**, Spearman's correlation (ρ) of effect sizes (allelic fold change, aFC in log₂ scale) between cis-eQTL and matched allele-specific expression (ASE) loci in the liver (i) and brain (j). N indicates number of tested loci. The lines are fitted by a linear regression model using the *geom.smooth* function from *ggplot2* (v3.3.2) in R (v4.0.2). The shading represents the standard error of the fitting line. **k**, Spearman's correlation (ρ) of effect sizes between cis-eQTL and matched ASE loci across 34 tissues. Red dots indicate number of tested loci (right y-axis).
a

<figure><img src="image_43.png" /></figure>

b

<figure><img src="image_44.png" /></figure>

d

<figure><img src="image_45.png" /></figure>

c

<figure><img src="image_46.png" /></figure>

e

<figure><img src="image_47.png" /></figure>

f

<figure><img src="image_48.png" /></figure>

g

<figure><img src="image_49.png" /></figure>

h

<figure><img src="image_50.png" /></figure>

i

<figure><img src="image_51.png" /></figure>

Extended Data Fig. 7 | See next page for caption.
a

<figure><img src="image_41.png" /></figure>

b

<figure><img src="image_42.png" /></figure>

d

<figure><img src="image_43.png" /></figure>

c

<figure><img src="image_44.png" /></figure>

e

<figure><img src="image_45.png" /></figure>

f

<figure><img src="image_46.png" /></figure>

g

<figure><img src="image_47.png" /></figure>

h

<figure><img src="image_48.png" /></figure>

i

<figure><img src="image_49.png" /></figure>

Extended Data Fig. 7 | See next page for caption.
**Extended Data Fig. 7 | Breed sharing and interaction cis-eQTL (bieQTL).**

a, Sample size of muscle RNA-Seq data across eight breed groups. b, c, Expression levels of *NMNAT1* (b) and *COMMD10* (c) at three genotypes of *cis*-eQTL in muscle across eight breed groups. d, The *cis*-eQTL discovered in each breed group (rows) that can be replicated ($\pi_i$) across all other breed groups (columns). e, The heatmap of tissues regarding the pairwise Spearman's correlation (*ρ*) of *cis*-eQTL effect sizes. Tissues are grouped by hierarchical clustering (bottom). Violin plot (left) represents Spearman's correlation between the target group and the rest. f, Pearson's correlation (*r*) of effect size between *cis*-eQTL from the multi-breed

meta-analysis (y-axis) and those from the combined muscle population (x-axis). The *P* value was obtained from Pearson's *r* test. g, Overlap of *cis*-eQTL detected from the combined muscle population (Combined) and those detected in single-breed (Single) *cis*-eQTL mapping (shared in at least two breeds). h, i, Examples of bieQTL in muscle. Each dot in (h, CA14) and (i, ATE1) represents an individual and is colored by three genotypes. Gene expression levels and ancestry enrichment scores are inverse normal transformed. The two-sided *P* value is calculated by the linear regression bieQTL model. The lines are fitted by a linear regression model using the *geom.smooth* function from *ggplot2* (v3.3.2) in R (v4.0.2).
**Extended Data Fig. 7 | Breed sharing and interaction cis-eQTL (bieQTL).**

a, Sample size of muscle RNA-Seq data across eight breed groups. b, c, Expression levels of *NMNAT1* (b) and *COMMD10* (c) at three genotypes of *cis*-eQTL in muscle across eight breed groups. d, The *cis*-eQTL discovered in each breed group (rows) that can be replicated ($\pi_i$) across all other breed groups (columns). e, The heatmap of tissues regarding the pairwise Spearman's correlation (*ρ*) of *cis*-eQTL effect sizes. Tissues are grouped by hierarchical clustering (bottom). Violin plot (left) represents Spearman's correlation between the target group and the rest. f, Pearson's correlation (*r*) of effect size between *cis*-eQTL from the multi-breed

meta-analysis (y-axis) and those from the combined muscle population (x-axis). The *P* value was obtained from Pearson's *r* test. g, Overlap of *cis*-eQTL detected from the combined muscle population (Combined) and those detected in single-breed (Single) *cis*-eQTL mapping (shared in at least two breeds). h, i, Examples of bieQTL in muscle. Each dot in (h, CA14) and (i, ATE1) represents an individual and is colored by three genotypes. Gene expression levels and ancestry enrichment scores are inverse normal transformed. The two-sided *P* value is calculated by the linear regression bieQTL model. The lines are fitted by a linear regression model using the *geom.smooth* function from *ggplot2* (v3.3.2) in R (v4.0.2).
<figure><img src="image_52.png" /><figcaption>**Extended Data Fig. 8** | Cell-type enrichment and interaction cis-eQTL (cieQTL). a, Distribution of enrichment scores (percentage) of major cell types in samples of seven tested tissues (brain: n = 415, frontal cortex: n = 75, hypothalamus: n = 73, lung: n = 149, blood: n = 386, liver: n = 501, and spleen: n = 91). Each point and whisker indicate the mean value and standard deviation, respectively. b,c, Examples of cieQTL in blood. Each dot in (b, SCRN2) and (c, HIBADH) represents an individual and is colored by three genotypes. Gene expression levels and cell-type enrichment scores are inverse normal transformed. The two-sided *P* value was calculated by the linear regression</figcaption></figure>

cieQTL model. The lines are fitted by a linear regression model using the *geomhalys* smooth function from ggplot2 (v3.3.2) in R (v4.0.2). *d-f*, Pearson's correlation (*r*) between allele-specific expression (ASE) effect sizes (allelic fold change, aFC) and specific cell-type enrichment scores for FGD2 with monocytes (*d*), SCRN2 with CD2<sup>-</sup> γδ T cells (*e*) and HIBADH with CD4<sup>+</sup> αβ T cells in the blood (*f*). The lines are fitted by a linear regression model using the *geomhalys* smooth function from ggplot2 (v3.3.2) in R (v4.0.2). The shading represents the standard error of the fitting line. *g*, ASE validation rate (π<sub>i</sub>) of breed/cell-type interaction QTL (bieQTL and cieQTL) across tissues with ≥ 5 detectable bieQTL or cieQTL.
<figure><img src="image_50.png" /><figcaption>**Extended Data Fig. 8** | Cell-type enrichment and interaction cis-eQTL (cieQTL). a, Distribution of enrichment scores (percentage) of major cell types in samples of seven tested tissues (brain: n = 415, frontal cortex: n = 75, hypothalamus: n = 73, lung: n = 149, blood: n = 386, liver: n = 501, and spleen: n = 91). Each point and whisker indicate the mean value and standard deviation, respectively. b,c, Examples of cieQTL in blood. Each dot in (b, SCRN2) and (c, HIBADH) represents an individual and is colored by three genotypes. Gene expression levels and cell-type enrichment scores are inverse normal transformed. The two-sided *P* value was calculated by the linear regression</figcaption></figure>

cieQTL model. The lines are fitted by a linear regression model using the *geomhalys* smooth function from ggplot2 (v3.3.2) in R (v4.0.2). **d-f**, Pearson's correlation (*r*) between allele-specific expression (ASE) effect sizes (allelic fold change, aFC) and specific cell-type enrichment scores for FGD2 with monocytes (**d**), SCRN2 with CD2<sup>-</sup> γδ T cells (**e**) and HIBADH with CD4<sup>+</sup> αβ T cells in the blood (**f**). The lines are fitted by a linear regression model using the *geomhalys* smooth function from ggplot2 (v3.3.2) in R (v4.0.2). The shading represents the standard error of the fitting line. **g**, ASE validation rate (π<sub>i</sub>) of breed/cell-type interaction QTL (bieQTL and cieQTL) across tissues with ≥ 5 detectable bieQTL or cieQTL.
<figure><img src="image_53.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_60.png" /></figure>

Extended Data Fig. 9 | See next page for caption.
<figure><img src="image_51.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_58.png" /></figure>

Extended Data Fig. 9 | See next page for caption.
**Extended Data Fig. 9** | **Tissue-sharing and specificity patterns of molecular QTL (molQTL).** a–d, The heatmap of tissues regarding the pairwise Spearman's correlation (ρ) of molQTL effect sizes, that is, cis-sQTL (a), cis-eeQTL (b), cis-lncQTL (c) and cis-enQTL (d). Tissues are grouped by the hierarchical clustering (bottom). Violin plot (left) represents Spearman's correlations between the target tissue and the rest. e, Distribution of number of tissues having METASOFT activity (m-value > 0.7) for each of molQTL. MolPhe: molecular phenotype. f, Pearson's correlation (r) between number of tissues an eGene expressed in

(transcript per million, TPM > 0.1) and its *cis*-eQTL effect sizes (|aFC(log₂)|). The aFC is for allelic fold change. The line and shading indicate the median and interquartile range, respectively. g, Expression levels (adjusted TMM) of *ODF2L* at three genotypes of top *cis*-eQTL (rs329043485) in blood and testis. TMM: trimmed mean of M-value normalized expression levels. There are 337, 47 and 2 samples for A/A, A/C and C/C genotypes in blood, respectively, and 148, 34 and 2 in testis, respectively. h, Expression levels (log₂TMM) of *ODF2L* across 34 tissues. Tissues are ordered (from smallest to largest) by the median expression values.
**Extended Data Fig. 9** | **Tissue-sharing and specificity patterns of molecular QTL (molQTL).** a–d, The heatmap of tissues regarding the pairwise Spearman's correlation (ρ) of molQTL effect sizes, that is, cis-sQTL (a), cis-eeQTL (b), cis-lncQTL (c) and cis-enQTL (d). Tissues are grouped by the hierarchical clustering (bottom). Violin plot (left) represents Spearman's correlations between the target tissue and the rest. e, Distribution of number of tissues having METASOFT activity (m-value > 0.7) for each of molQTL. MolPhe: molecular phenotype. f, Pearson's correlation (r) between number of tissues an eGene expressed in

(transcript per million, TPM > 0.1) and its cis-eQTL effect sizes (IaFC(log₂I)). The aFC is for allelic fold change. The line and shading indicate the median and interquartile range, respectively. g, Expression levels (adjusted TMM) of *ODF2L* at three genotypes of top cis-eQTL (rs329043485) in blood and testis. TMM: trimmed mean of M-value normalized expression levels. There are 337, 47 and 2 samples for A/A, A/C and C/C genotypes in blood, respectively, and 148, 34 and 2 in testis, respectively. h, Expression levels (log₂TMM) of *ODF2L* across 34 tissues. Tissues are ordered (from smallest to largest) by the median expression values.
<figure><img src="image_61.png" /><figcaption>Extended Data Fig. 10 | See next page for caption.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Chromosome 13 position (Mb)</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Chromosome 15 position (Mb)</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Chromosome 2 position (Mb)</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>GOSR2
MSTRG.4694&ENSSSCT000000070888
Chromosome 12 position (Mb)</figcaption></figure>

Extended Data Fig. 10 | See next page for caption.
**Extended Data Fig. 10** | Complementarity of molecular QTL (molQTL) in interpreting GWAS loci. a, Number of GWAS loci linked to *cis*-eQTL, *cis*-sQTL, *cis*-eeQTL, *cis*-lncQTL and *cis*-enQTL in 34 tissues based on four different integrative methods, including colocalization (fastEnloc), Mendelian randomization (SMR), single-tissue transcriptome-wide association studies (TWAS, S-PrediXcan) and multi-tissue TWAS (S-MultiXcan). The bottom point-line combinations of the Upset plot represent the intersections of GWAS loci linked to eGenes by different types of molecular phenotypes. b, Distribution of rank correlations between tissue-relevance-scores derived from *cis*-eQTL and those from *cis*-sQTL, *cis*-lncQTL, *cis*-eeQTL and *cis*-enQTL across 86 GWAS traits with significant colocalizations for at least one molecular phenotype. c, Significant SMR signals (*P*<sub>SMR</sub> = 9.16 × 10<sup>-5</sup>, *P*<sub>HEIDI</sub> = 0.9) between GWAS loci of average daily gain (ADG) and *cis*-eQTL of *CFAP298-TCP10L* in colon, but not for its *cis*-sQTL or *cis*-eeQTL. The orange triangle represents the top *cis*-eQTL of *CFAP298-TCP10L*. d, Significant

SMR signals ($P_{\text{SMR}} = 1.78 \times 10^{-5}$, $P_{\text{HEIDI}} = 0.07$) between GWAS loci of the average backfat thickness (BFT) and *cis*-sQTL of *MYO7B* in the small intestine, but not for its *cis*-eQTL or *cis*-eeQTL. e, Significant SMR signals ($P_{\text{SMR}} = 1.78 \times 10^{-6}$, $P_{\text{HEIDI}} = 0.97$) between GWAS loci of litter weight (LW, piglets born alive) and *cis*-eeQTL of *FBXL12* in the uterus, but not for its *cis*-eQTL or *cis*-sQTL. f, Significant SMR signals ($P_{\text{SMR(lncQTL-GWAS)}} = 4.49 \times 10^{-7}$, $P_{\text{SMR(eQTL-GWAS)}} = 5.45 \times 10^{-5}$, $P_{\text{SMR(lncQTL-eQTL)}} = 4.62 \times 10^{-7}$) among GWAS loci of loin muscle depth (LMD), *cis*-lncQTL of *MSTRG.4694&ENSSSCT00000070888*, and *cis*-eQTL of *GOSR2* in the muscle. *MSTRG.4694&ENSSSCT00000070888* is a lncRNA gene located on the 3112 bp downstream of *GOSR2*, where the Pearson's correlation of their normalized expression levels (trimmed mean of M-value, TMM) is 0.29 in muscle. The orange and green triangles in the top GWAS Manhattan plot represent the top molQTL of *GOSR2* and *MSTRG.4694&ENSSSCT00000070888*, respectively.
**Extended Data Fig. 10** | Complementarity of molecular QTL (molQTL) in interpreting GWAS loci. a, Number of GWAS loci linked to *cis*-eQTL, *cis*-sQTL, *cis*-eeQTL, *cis*-lncQTL and *cis*-enQTL in 34 tissues based on four different integrative methods, including colocalization (fastEnloc), Mendelian randomization (SMR), single-tissue transcriptome-wide association studies (TWAS, S-PrediXcan) and multi-tissue TWAS (S-MultiXcan). The bottom point-line combinations of the Upset plot represent the intersections of GWAS loci linked to eGenes by different types of molecular phenotypes. b, Distribution of rank correlations between tissue-relevance-scores derived from *cis*-eQTL and those from *cis*-sQTL, *cis*-lncQTL, *cis*-eeQTL and *cis*-enQTL across 86 GWAS traits with significant colocalizations for at least one molecular phenotype. c, Significant SMR signals (*P*<sub>SMR</sub> = 9.16 × 10<sup>-5</sup>, *P*<sub>HEIDI</sub> = 0.9) between GWAS loci of average daily gain (ADG) and *cis*-eQTL of *CFAP298-TCP10L* in colon, but not for its *cis*-sQTL or *cis*-eeQTL. The orange triangle represents the top *cis*-eQTL of *CFAP298-TCP10L*. d, Significant

SMR signals ($P_{\text{SMR}} = 1.78 \times 10^{-5}$, $P_{\text{HEIDI}} = 0.07$) between GWAS loci of the average backfat thickness (BFT) and *cis*-sQTL of *MYO7B* in the small intestine, but not for its *cis*-eQTL or *cis*-eeQTL. e, Significant SMR signals ($P_{\text{SMR}} = 1.78 \times 10^{-6}$, $P_{\text{HEIDI}} = 0.97$) between GWAS loci of litter weight (LW, piglets born alive) and *cis*-eeQTL of *FBXL12* in the uterus, but not for its *cis*-eQTL or *cis*-sQTL. f, Significant SMR signals ($P_{\text{SMR(lncQTL-GWAS)}} = 4.49 \times 10^{-7}$, $P_{\text{SMR(eQTL-GWAS)}} = 5.45 \times 10^{-5}$, $P_{\text{SMR(lncQTL-eQTL)}} = 4.62 \times 10^{-7}$) among GWAS loci of loin muscle depth (LMD), *cis*-lncQTL of *MSTRG.4694&ENSSSCT00000070888*, and *cis*-eQTL of *GOSR2* in the muscle. *MSTRG.4694&ENSSSCT00000070888* is a lncRNA gene located on the 3112 bp downstream of *GOSR2*, where the Pearson's correlation of their normalized expression levels (trimmed mean of M-value, TMM) is 0.29 in muscle. The orange and green triangles in the top GWAS Manhattan plot represent the top molQTL of *GOSR2* and *MSTRG.4694&ENSSSCT00000070888*, respectively.
<figure><img src="image_62.png" /></figure>

## Reporting Summary

Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our *Editorial Policies* and the *Editorial Policy Checklist*.

## Statistics

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a | Confirmed

☐ ☑ The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

☐ ☑ A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

☐ ☑ The statistical test(s) used AND whether they are one- or two-sided
*Only common tests should be described solely by name; describe more complex techniques in the Methods section.*

☐ ☑ A description of all covariates tested

☐ ☑ A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

☐ ☑ A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

☐ ☑ For null hypothesis testing, the test statistic (e.g. *F*, *t*, *r*) with confidence intervals, effect sizes, degrees of freedom and *P* value noted
*Give *P* values as exact values whenever suitable.*

☐ ☑ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

☐ ☑ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

☐ ☑ Estimates of effect sizes (e.g. Cohen's *d*, Pearson's *r*), indicating how they were calculated

*Our web collection on statistics for biologists contains articles on many of the points above.*

## Software and code

Policy information about availability of computer code

Data collection

All raw data analyzed in this study are publicly available for download without restrictions from SRA (https://www.ncbi.nlm.nih.gov/sra/) and BIGD (https://bigd.big.ac.cn/bioproject/) databases using the wget function in Linux. Details of RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets can be found in Supplementary Table 1, 2, 5, 8 and 9, respectively.

Data analysis

All the computational scripts and codes (with software version) for RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets analyses, as well as the respective quality control, molecular phenotype normalization, genotype imputation, molQTL mapping, functional enrichment, colocalization, SMR and TWAS are available at the FarmGTEx GitHub website (https://github.com/FarmGTEx/PigGTEx-Pipeline-v0, https://doi.org/10.6084/m9.fishshare.24247771).

For RNA-Seq data analysis, we used Trimmomatic (v0.39), STAR (v2.7.0), Stringtie (v2.1.1), featureCounts (v1.5.2), Leafcutter (v0.2.9), GATK (v4.0.8.1), phASER (v1.1.1), and Beagle (v5.1) for quality control, mapping, gene expression quantification, alternative splicing, SNP calling, ASE analysis, and genotype imputation, respectively. For sample clustering, we used MEGA (vX) and then visualized with iTOL (v6). For tissue-specific gene expression, we used limma (v3.51.2). For gene co-expression analysis, we used WGCNA (v1.69), ICA (v1.0.2), PEER (v1.3), MEGENA (v1.3.7), and CEMiTool (v1.8.3). For gene functional enrichment analysis, we used clusterProfiler (v4.0) and visualized it using Gephi (v0.9.2).

For WGS analysis, we used Trimmomatic (v0.39), BWA-MEM (v0.7.5a-r405), Picard (v2.21.2), GATK (v4.1.4.1), and Beagle (v5.1) for quality control, mapping, marked duplicated reads, variants calling, and phasing, respectively. We used used PLINK (v1.90) to do LD pruning.

For WGBS, we used FastQC (v0.11.9), Trim Galore (v0.4.5), Bismark (v0.19.0), and SMART2 (v2.2.8), Methpipe (v4.1.1), and FastQTL (v2.184) for quality evaluation, quality control, read mapping and DNA methylation level extraction, hypomethylation region detection, allele-specific
<figure><img src="image_65.png" /></figure>

# Reporting Summary

Nature Portfolio wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Portfolio policies, see our [Editorial Policies](#) and the [Editorial Policy Checklist](#).

## Statistics

For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a | Confirmed

☐ ☑ The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement

☐ ☑ A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly

☐ ☑ The statistical test(s) used AND whether they are one- or two-sided
*Only common tests should be described solely by name; describe more complex techniques in the Methods section.*

☐ ☑ A description of all covariates tested

☐ ☑ A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons

☐ ☑ A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)

☐ ☑ For null hypothesis testing, the test statistic (e.g. *F*, *t*, *r*) with confidence intervals, effect sizes, degrees of freedom and *P* value noted
*Give *P* values as exact values whenever suitable.*

☐ ☑ For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings

☐ ☑ For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes

☐ ☑ Estimates of effect sizes (e.g. Cohen's *d*, Pearson's *r*), indicating how they were calculated

*Our web collection on statistics for biologists contains articles on many of the points above.*

## Software and code

Policy information about [availability of computer code](#)

**Data collection**

All raw data analyzed in this study are publicly available for download without restrictions from SRA (https://www.ncbi.nlm.nih.gov/sra/) and BIGD (https://bigd.big.ac.cn/bioproject/) databases using the wget function in Linux. Details of RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets can be found in Supplementary Table 1, 2, 5, 8 and 9, respectively.

**Data analysis**

All the computational scripts and codes (with software version) for RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets analyses, as well as the respective quality control, molecular phenotype normalization, genotype imputation, molQTL mapping, functional enrichment, colocalization, SMR and TWAS are available at the FarmGTEx GitHub website (https://github.com/FarmGTEx/PigGTEx-Pipeline-v0, https://doi.org/10.6084/m9.fishshare.24247771).

For RNA-Seq data analysis, we used Trimmomatic (v0.39), STAR (v2.7.0), Stringtie (v2.1.1), featureCounts (v1.5.2), Leafcutter (v0.2.9), GATK (v4.0.8.1), phASER (v1.1.1), and Beagle (v5.1) for quality control, mapping, gene expression quantification, alternative splicing, SNP calling, ASE analysis, and genotype imputation, respectively. For sample clustering, we used MEGA (vX) and then visualized with iTOL (v6). For tissue-specific gene expression, we used limma (v3.51.2). For gene co-expression analysis, we used WGCNA (v1.69), ICA (v1.0.2), PEER (v1.3), MEGENA (v1.3.7), and CEMiTool (v1.8.3). For gene functional enrichment analysis, we used clusterProfiler (v4.0) and visualized it using Gephi (v0.9.2).

For WGS analysis, we used Trimmomatic (v0.39), BWA-MEM (v0.7.5a-r405), Picard (v2.21.2), GATK (v4.1.4.1), and Beagle (v5.1) for quality control, mapping, marked duplicated reads, variants calling, and phasing, respectively. We used used PLINK (v1.90) to do LD pruning.

For WGBS, we used FastQC (v0.11.9), Trim Galore (v0.4.5), Bismark (v0.19.0), and SMART2 (v2.2.8), Methpipe (v4.1.1), and FastQTL (v2.184) for quality evaluation, quality control, read mapping and DNA methylation level extraction, hypomethylation region detection, allele-specific
methylation loci analysis and methylation QTL mapping, respectively.

For Hi-C, we used Trim Galore (v0.6.7), BWA (v0.7.17), Juicer (v1.6), Arrowhead (v1.22.01), hicConvertFormat (v3.7.1), pyGenomeTracks (v3.6) for quality control, read mapping, Hi-C contact matrix construction, TAD identification, format conversion, and visualization, respectively.

For single-cell RNA-Seq, we used Seurat (v3.0.2), Azimuth (v0.4.0) and CIBERSORTx online tool (v1) for data processing, cell type annotation and cell type deconvolution, respectively.

We removed SNPs with MAF < 0.01 and/or missing rate > 0.9 using bcftools (v1.9) and employed Beagle (v5.1) to phase the filtered variants and impute sporadically missing genotypes.

For QTL mapping, we used TensorQTL (v1.0.3), aFC (v0.3), dap-g (v1.0.0), METASOFT (v2.0.1), MashR (v0.2-6), and GCTA (v1.93.0) for cis-QTL mapping, effect size estimation, fine-mapping, meta-analysis, tissue-sharing pattern estimation, and cis-QTL mapping with mixed linear model, respectively. We estimated the genetic parameters using the restricted maximum likelihood (REML) method implemented in GCTA (v1.93.0).

We computed genotype PCs based on the filtered SNPs within each of the tissues using SNPRelate (v1.26.0). To account for technical confounders among RNA-Seq samples (e.g., hidden batch effects and other technical or biological factors), we used the Probabilistic Estimation of Expression Residuals (PEER) method, implemented in peer R package (v4.0.2), to estimate a set of latent covariates within each of the 34 tissues based on gene expression matrices. We computed the mappability of each locus in the reference genome using GenMap (v1.3.0). We removed SNPs in repeat regions annotated by the UCSC RepeatMasker track.

We first used imputed genotypes to estimate the ancestry composition of all RNA-Seq samples across tissues using ADMIXTURE (v1.3.0). We estimated the effect size (aFC) of the top ieQTL of ieGenes from ASE data using the script phaser_cis_var.py in phASER (v1.1.1).

For integrative analysis between GWAS and molQTL, we used S-PrediXcan and S-MultiXcan in MetaXcan (v0.6.11) for single-tissue and multi-tissue TWAS analysis, SMR (v1.03) for Mendelian Randomization analysis, and fastENLOC (v1.0) for colocalization. We performed a meta-analysis of molQTL across all 34 tissues using MashR (v0.2-6) and METASOFT (v2.0.1). We calculated the pairwise Rand index to measure the clustering similarity using the rand.index function in the fossil (v0.4.0) R package (v4.0.2).

we performed 2,056 separate GWAS, and conducted the meta-GWAS analysis for the same traits across different populations based on GWAS summary statistics using METAL (v2011-03-25).

For enrichment analysis, we used TORUS (v1) and ClusterProfiler (v4.0) for molQTLs and genes functional annotation, respectively.

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

## Data

Policy information about availability of data

All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets

- A description of any restrictions on data availability

- For clinical datasets or third party data, please ensure that the statement adheres to our policy

All raw data analyzed in this study are publicly available for download without restrictions from SRA (https://www.ncbi.nlm.nih.gov/sra/) and BIGD (https://bigd.big.ac.cn/bioproject/) databases. Details of RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets can be found in Supplementary Tables 1, 2, 5, 8 and 9, respectively. All WGS data generated in this study are available under CNCB GSA (https://ngdc.cncb.ac.cn/) accessions: PRJCA016120, PRJCA016130, PRJCA017284, PRJCA016012, and PRJCA016216. All processed data and the full summary statistics of molQTL mapping are available at http://piggtex.farmgtex.org/.

## Human research participants

Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender NA

Population characteristics NA

Recruitment NA

Ethics oversight NA

Note that full information on the approval of the study protocol must also be provided in the manuscript.
methylation loci analysis and methylation QTL mapping, respectively.

For Hi-C, we used Trim Galore (v0.6.7), BWA (v0.7.17), Juicer (v1.6), Arrowhead (v1.22.01), hicConvertFormat (v3.7.1), pyGenomeTracks (v3.6) for quality control, read mapping, Hi-C contact matrix construction, TAD identification, format conversion, and visualization, respectively.

For single-cell RNA-Seq, we used Seurat (v3.0.2), Azimuth (v0.4.0) and CIBERSORTx online tool (v1) for data processing, cell type annotation and cell type deconvolution, respectively.

We removed SNPs with MAF < 0.01 and/or missing rate > 0.9 using bcftools (v1.9) and employed Beagle (v5.1) to phase the filtered variants and impute sporadically missing genotypes.

For QTL mapping, we used TensorQTL (v1.0.3), aFC (v0.3), dap-g (v1.0.0), METASOFT (v2.0.1), MashR (v0.2-6), and GCTA (v1.93.0) for cis-QTL mapping, effect size estimation, fine-mapping, meta-analysis, tissue-sharing pattern estimation, and cis-QTL mapping with mixed linear model, respectively. We estimated the genetic parameters using the restricted maximum likelihood (REML) method implemented in GCTA (v1.93.0).

We computed genotype PCs based on the filtered SNPs within each of the tissues using SNPRelate (v1.26.0). To account for technical confounders among RNA-Seq samples (e.g., hidden batch effects and other technical or biological factors), we used the Probabilistic Estimation of Expression Residuals (PEER) method, implemented in peer R package (v4.0.2), to estimate a set of latent covariates within each of the 34 tissues based on gene expression matrices. We computed the mappability of each locus in the reference genome using GenMap (v1.3.0). We removed SNPs in repeat regions annotated by the UCSC RepeatMasker track.

We first used imputed genotypes to estimate the ancestry composition of all RNA-Seq samples across tissues using ADMIXTURE (v1.3.0). We estimated the effect size (aFC) of the top ieQTL of ieGenes from ASE data using the script phaser_cis_var.py in phASER (v1.1.1).

For integrative analysis between GWAS and molQTL, we used S-PrediXcan and S-MultiXcan in MetaXcan (v0.6.11) for single-tissue and multi-tissue TWAS analysis, SMR (v1.03) for Mendelian Randomization analysis, and fastENLOC (v1.0) for colocalization. We performed a meta-analysis of molQTL across all 34 tissues using MashR (v0.2-6) and METASOFT (v2.0.1). We calculated the pairwise Rand index to measure the clustering similarity using the rand.index function in the fossil (v0.4.0) R package (v4.0.2).

we performed 2,056 separate GWAS, and conducted the meta-GWAS analysis for the same traits across different populations based on GWAS summary statistics using METAL (v2011-03-25).

For enrichment analysis, we used TORUS (v1) and ClusterProfiler (v4.0) for molQTLs and genes functional annotation, respectively.

For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Portfolio guidelines for submitting code & software for further information.

## Data

Policy information about availability of data

All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:

- Accession codes, unique identifiers, or web links for publicly available datasets

- A description of any restrictions on data availability

- For clinical datasets or third party data, please ensure that the statement adheres to our policy

All raw data analyzed in this study are publicly available for download without restrictions from SRA (https://www.ncbi.nlm.nih.gov/sra/) and BIGD (https://bigd.big.ac.cn/bioproject/) databases. Details of RNA-Seq, WGS, WGBS, single-cell RNA-Seq and Hi-C datasets can be found in Supplementary Tables 1, 2, 5, 8 and 9, respectively. All WGS data generated in this study are available under CNCB GSA (https://ngdc.cncb.ac.cn/) accessions: PRJCA016120, PRJCA016130, PRJCA017284, PRJCA016012, and PRJCA016216. All processed data and the full summary statistics of molQTL mapping are available at http://piggtex.farmgtex.org/.

## Human research participants

Policy information about studies involving human research participants and Sex and Gender in Research.

Reporting on sex and gender NA

Population characteristics NA

Recruitment NA

Ethics oversight NA

Note that full information on the approval of the study protocol must also be provided in the manuscript.
# Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

<figure><img src="image_63.png" /></figure>

Life sciences

<figure><img src="image_64.png" /></figure>

Behavioural & social sciences

<figure><img src="image_65.png" /></figure>

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://nature.com/documents/nr-reporting-summary-flat.pdf)

# Life sciences study design

All studies must disclose on these points even when the disclosure is negative.

**Sample size**

No power calculation was needed in advance in this study. In total, we analyzed all 11,323 RNA-Seq runs (downloaded by March, 2021) from SRA (https://www.ncbi.nlm.nih.gov/sra/), and BIGD databases (https://bigd.big.ac.cn/bioproject/), yielding 9,530 unique RNA-Seq samples. After filtering the samples with low quality (see below), all the remaining samples have been used for analysis.

**Data exclusions**

Full details of data exclusions for each analysis can be found in the Methods as well.

We filtered out RNA-Seq samples with clean read counts ≤ 500K or uniquely mapping rates < 60%, resulting in 8,262 samples. We further excluded samples with obvious clustering errors (e.g., samples labeled as liver that were not clustered with other liver samples), resulting in 7,095 samples for subsequent analysis.

For cis-QTLs detection, we excluded tissues with less than 40 individuals, resulting in 34 tissues for cis-QTL mapping.

**Replication**

To validate the cis-eQTLs, we applied four distinct strategies including linear mixed model, internal validation, external validation, and ASE validation.

First, we observed that the summary statistics of cis-eQTL derived from the linear regression model in TensorQTL had a strong correlation (an average Pearson's r of 0.91 across tissues) with those from a linear mixed model.

Second, we performed an internal validation in 18 tissues with over 80 samples by randomly dividing samples into two equal groups, and then conducting cis-eQTL mapping separately in both subgroups. We observed a high replication rate (an average π1 of 0.92) for cis-eQTL discovery.

Third, we found that 92%, 74%, 73%, and 69% of cis-eQTL in blood, liver, duodenum, and muscle, respectively, were replicated in independent datasets.

Fourth, we further found that effects (allelic fold changes, aFC) derived from allele specific expression (ASE) analysis were significantly correlated with those from cis-eQTL mapping consistently across tissues. For instance, in muscle, ASE-derived effects of 4,417 SNPs were significantly correlated (Spearman's ρ = 0.76, P < 1e-300) with their cis-eQTL effects.

**Randomization**

All the datasets are from observation studies and we used all samples publicly available after data exclusions listed above. Therefore, Randomization were not relevant in this study. Samples were grouped by tissue types.

**Blinding**

In this study, we re-analyzed all the publicly available RNA-seq data using a uniform pipeline, followed by the population-based association studies and validated the findings in independent populations. The blinding study design may be not applicable in this study.

# Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

## Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>Antibodies</td></tr><tr><td>☒</td><td>Eukaryotic cell lines</td></tr><tr><td>☒</td><td>Palaeontology and archaeology</td></tr><tr><td>☒</td><td>Animals and other organisms</td></tr><tr><td>☒</td><td>Clinical data</td></tr><tr><td>☒</td><td>Dual use research of concern</td></tr></tbody></table>

## Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>ChIP-seq</td></tr><tr><td>☒</td><td>Flow cytometry</td></tr><tr><td>☒</td><td>MRI-based neuroimaging</td></tr></tbody></table>
# Field-specific reporting

Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

<figure><img src="image_66.png" /></figure>

Life sciences

<figure><img src="image_67.png" /></figure>

Behavioural & social sciences

<figure><img src="image_68.png" /></figure>

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see [nature.com/documents/nr-reporting-summary-flat.pdf](https://nature.com/documents/nr-reporting-summary-flat.pdf)

# Life sciences study design

All studies must disclose on these points even when the disclosure is negative.

## Sample size

No power calculation was needed in advance in this study. In total, we analyzed all 11,323 RNA-Seq runs (downloaded by March, 2021) from SRA (https://www.ncbi.nlm.nih.gov/sra/), and BIGD databases (https://bigd.big.ac.cn/bioproject/), yielding 9,530 unique RNA-Seq samples. After filtering the samples with low quality (see below), all the remaining samples have been used for analysis.

## Data exclusions

Full details of data exclusions for each analysis can be found in the Methods as well.

We filtered out RNA-Seq samples with clean read counts ≤ 500K or uniquely mapping rates < 60%, resulting in 8,262 samples. We further excluded samples with obvious clustering errors (e.g., samples labeled as liver that were not clustered with other liver samples), resulting in 7,095 samples for subsequent analysis.

For cis-QTLs detection, we excluded tissues with less than 40 individuals, resulting in 34 tissues for cis-QTL mapping.

## Replication

To validate the cis-eQTLs, we applied four distinct strategies including linear mixed model, internal validation, external validation, and ASE validation.

First, we observed that the summary statistics of cis-eQTL derived from the linear regression model in TensorQTL had a strong correlation (an average Pearson's r of 0.91 across tissues) with those from a linear mixed model.

Second, we performed an internal validation in 18 tissues with over 80 samples by randomly dividing samples into two equal groups, and then conducting cis-eQTL mapping separately in both subgroups. We observed a high replication rate (an average π1 of 0.92) for cis-eQTL discovery.

Third, we found that 92%, 74%, 73%, and 69% of cis-eQTL in blood, liver, duodenum, and muscle, respectively, were replicated in independent datasets.

Fourth, we further found that effects (allelic fold changes, aFC) derived from allele specific expression (ASE) analysis were significantly correlated with those from cis-eQTL mapping consistently across tissues. For instance, in muscle, ASE-derived effects of 4,417 SNPs were significantly correlated (Spearman's ρ = 0.76, P < 1e-300) with their cis-eQTL effects.

## Randomization

All the datasets are from observation studies and we used all samples publicly available after data exclusions listed above. Therefore, Randomization were not relevant in this study. Samples were grouped by tissue types.

## Blinding

In this study, we re-analyzed all the publicly available RNA-seq data using a uniform pipeline, followed by the population-based association studies and validated the findings in independent populations. The blinding study design may be not applicable in this study.

# Reporting for specific materials, systems and methods

We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

## Materials & experimental systems

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>Antibodies</td></tr><tr><td>☒</td><td>Eukaryotic cell lines</td></tr><tr><td>☒</td><td>Palaeontology and archaeology</td></tr><tr><td>☒</td><td>Animals and other organisms</td></tr><tr><td>☒</td><td>Clinical data</td></tr><tr><td>☒</td><td>Dual use research of concern</td></tr></tbody></table>

## Methods

<table><thead><tr><th>n/a</th><th>Involved in the study</th></tr></thead><tbody><tr><td>☒</td><td>ChIP-seq</td></tr><tr><td>☒</td><td>Flow cytometry</td></tr><tr><td>☒</td><td>MRI-based neuroimaging</td></tr></tbody></table>