- Research article
- Open Access
About the existence of common determinants of gene expression in the porcine liver and skeletal muscle
BMC Genomicsvolume 20, Article number: 518 (2019)
The comparison of expression QTL (eQTL) maps obtained in different tissues is an essential step to understand how gene expression is genetically regulated in a context-dependent manner. In the current work, we have compared the transcriptomic and eQTL profiles of two porcine tissues (skeletal muscle and liver) which typically show highly divergent expression profiles, in 103 Duroc pigs genotyped with the Porcine SNP60 BeadChip (Illumina) and with available microarray-based measurements of hepatic and muscle mRNA levels. Since structural variation could have effects on gene expression, we have also investigated the co-localization of cis-eQTLs with copy number variant regions (CNVR) segregating in this Duroc population.
The analysis of differential expresssion revealed the existence of 1204 and 1490 probes that were overexpressed and underexpressed in the gluteus medius muscle when compared to liver, respectively (|fold-change| > 1.5, q-value < 0.05). By performing genome scans in 103 Duroc pigs with available expression and genotypic data, we identified 76 and 28 genome-wide significant cis-eQTLs regulating gene expression in the gluteus medius muscle and liver, respectively. Twelve of these cis-eQTLs were shared by both tissues (i.e. 42.8% of the cis-eQTLs identified in the liver were replicated in the gluteus medius muscle). These results are consistent with previous studies performed in humans, where 50% of eQTLs were shared across tissues. Moreover, we have identified 41 CNVRs in a set of 350 pigs from the same Duroc population, which had been genotyped with the Porcine SNP60 BeadChip by using the PennCNV and GADA softwares, but only a small proportion of these CNVRs co-localized with the cis-eQTL signals.
Despite the fact that there are considerable differences in the gene expression patterns of the porcine liver and skeletal muscle, we have identified a substantial proportion of common cis-eQTLs regulating gene expression in both tissues. Several of these cis-eQTLs influence the mRNA levels of genes with important roles in meat (CTSF) and carcass quality (TAPT1), lipid metabolism (TMEM97) and obesity (MARC2), thus evidencing the practical importance of dissecting the genetic mechanisms involved in their expression.
The performance of GWAS in humans has revealed that most of the regions that display significant associations with complex traits are not exonic, meaning that regulatory polymorphisms might have important effects on phenotypic variation . This realization has prompted the mapping of expression QTL (eQTL), i.e. single nucleotide polymorphisms (SNP), indels or copy number variants (CNV) that explain part of the variance of gene expression phenotypes . Such studies have revealed that the majority of eQTL exert their effects in cis- (i.e. on neighboring genes) . There are also evidences that CNV-tagging SNPs are enriched in cis-eQTLs and that they often modulate multiple expression traits . By examining the patterns of expression of 22,286 genes in 9 human tissues, the GTEx Consortium has shown that approximately 50% of eQTL are shared by the nine tissues and that most of them display consistent effects across tissues .
In pigs, hundreds of eQTLs with effects on muscle [5,6,7,8,9], liver [10, 11] and backfat  gene expression have been mapped. Often, these pig eQTL studies have targeted subsets of genes either mapping to QTL [12, 13] or displaying significant expression-phenotype correlations [5, 6, 10, 11]. The broad majority of porcine eQTL studies have targeted single anatomic locations and, in consequence, they do not provide clues about the differential genetic regulation of distinct tissues and organs. Moreover, these genome scans have explored the association of gene expression with the allelic variation of SNPs or microsatellites, neglecting the potential effect of CNVs on expression phenotypes. The goals of the current work were to compare the gene expression patterns and cis-eQTL landscapes of the pig gluteus medius (GM) skeletal muscle and liver, two tissues with highly differentiated patterns of expression , as well as to investigate the co-localization of cis-eQTLs and CNVRs.
Differential expression analysis
A total of 1204 probes were found to be overexpressed in the GM muscle (|FC| > 1.5; q-value < 0.05), whereas 1490 showed overexpression in the liver (Fig. 1 and Additional file 1). The list of genes showing the highest levels of differential expression between muscle and liver, included the bridging integrator 1 (BIN1, log2FC = 5.34, q-value = 2.21E-198), the sorbitol dehydrogenase (SORD, log2FC = − 5.82, q-value = 3.47E-185), the protein phosphatase 1 regulatory inhibitor subunit 1A (PPP1R1A, log2FC = 5.67, q-value = 5.62E-184), the zinc finger protein 106 (ZNF106, log2FC = 4.94, q-value = 1.60E-183), the hydroxysteroid 17-beta dehydrogenase 13 (HSD17B13, log2FC = − 5.54, q-value = 3.04E-180) and apolipoprotein A1 (APOA1, log2FC = − 5.88, q-value = 1.25E-179). Noteworthy, the log2FC values are very high and the q-values highly significant, thereby evidencing that the expression profiles of the porcine liver and the skeletal muscle are strongly divergent.
Detection of expression QTL in the porcine skeletal muscle and liver
A total of 76 and 28 genome-wide significant cis-eQTLs were identified in the GM muscle and liver, respectively (Tables 1 and 2). Several genes that are cis-regulated in the muscle participate in lipid or carbohydrate metabolism, e.g. HIV-1 Tat Interactive Protein 2 (HTATIP2) , exocyst complex component 7 (EXOC7)  acyl-CoA dehydrogenase short chain (ACADS) , insulin degrading enzyme (IDE) , solute carrier family 38 member 9 (SLC38A9) , and family with sequence similarity 3 member C (FAM3C) . With regard to the liver, hydroxysteroid dehydrogenase like 2 (HSDL2)  and lipase A (LIPA)  genes have also important roles in lipid metabolism. The level of overlap between the muscle and the liver cis-eQTL data sets was considerable, with 12 cis-eQTLs shared by both tissues (42.8% of the cis-eQTLs detected in the liver were also detected in the GM muscle). Four examples of shared cis-eQTLs are shown in Fig. 2. There was a strong consistency in the direction of the effects of the shared cis-eQTLs in both tissues. For instance, the cis-eQTL regulating ARL2BP and TMEM97 showed positive allelic effects in both tissues, whilst consistent negative allelic effects were estimated for the cis-eQTL modulating the hepatic and muscular expression of BSDC1, POLR3D and TAP1 (Tables 1 and 2). We also detected a number of cis-eQTLs that had effects on gene expression either in the muscle or in the liver (Figs. 3 and 4). RNA-Seq data from the GM muscle from 52 of the 103 pigs with GM muscle microarray was also available. Comparison of the genotypic means obtained with both platforms and corresponding to genes that are cis-regulated in the GM muscle provided highly consistent results, i.e. the direction of the difference between the two homozygous genotypes was concordant in 82% of the comparisons and the same ordering of genotypes (from highest to lowest expression) was obtained in 70% of the comparisons (Additional file 2). However, both methods did not always yield significant results (about 50% of coincidence with regard to statistical significance). The most probable reason for this is that the sample size of the RNA-Seq data set is half the size when compared to the microarray data set, which leads to larger standard deviations for the RNA-Seq datasets. This would consequently reduce the power to detect significant differences between the genotypic means.
Co-localization of copy number variants and eQTLs
The analysis of structural variation was performed with two softwares: PennCNV and GADA. We detected 93 and 103 CNVRs with PennCNV and GADA, respectively, and 44.08% of the CNVRs detected with PennCNV were also called by GADA. These 41 common CNVRs were distributed along 13 pig chromosomes (Additional file 3). The proportions of copy gain, loss and loss/gain CNVRs were ~ 48.78%, ~ 39.02% and ~ 12.19% respectively. The size of the CNVRs ranged between 31.4 kb and 5.2 Mb, with a mean of 457.4 kb. We compared our CNVR dataset with other CNVRs previously reported in pigs [22,23,24,25,26,27,28,29], and found that 60.97% of our CNVRs had been previously reported (Additional file 4). Real time quantitative assays were designed and used to validate 4 CNVRs (CNVR 9, 15, 32 and 38) in 39 porcine samples. According to D’haene , estimates of copy number between 1.414 and 2.449 most likely correspond to a normal copy number of 2, whilst anything below or above these thresholds might represent a deletion (CN = 1) or a duplication (CN = 3), respectively. Following these criteria, the four regions under analysis showed evidence of structural variation (Fig. 5). The co-localization of CNVRs and eQTLs was also analyzed (Additional file 5). In the GM muscle, 2 CNVRs co-localized with 3 cis-eQTLs, while 2 CNVRs co-localized with 2 hepatic cis-eQTLs. This low concordance between CNVRs and cis-eQTLs was not surprising as only 10 CNVRs co-localized with the genomic coordinates of the gene data set analyzed in our experiment.
The profiles of mRNA expression of the porcine liver and skeletal muscle are highly divergent
Skeletal muscle and liver have been selected as target tissues to perform a comparative eQTL analysis because they have a key role in body energy homeostasis, a parameter that has a strong impact on growth and fat deposition in pigs. The analysis of differential expression revealed that the patterns of expression of these two tissues are highly differentiated in pigs (Fig. 1), probably as a consequence of their distinct embryonic origins and physiological functions . Skeletal muscle and liver also show different patterns of gene expression in mouse  and humans . The strongly up-regulated genes in muscle included BIN1, which mediates calcium signaling and excitation–contraction coupling , tropomyosin (TPM2), which is fundamental for muscle contraction , and ZNF106, which is required for the postnatal maintenance of myofiber innervation by motor neurons . In the liver, we detected a strong overexpression of SORD, an enzyme that converts sorbitol into fructose , HSD17B13, whose inactivation leads to chronic liver disease , and APOA1, the major structural component of high-density lipoproteins .
Identification of cis-eQTLs with consistent effects on the mRNA levels of genes expressed in the liver and skeletal muscle
In pigs, hundreds of eQTLs with effects on the muscle transcriptome have been detected, but the majority of these studies focused on target subsets of genes either mapping to QTLs  or displaying significant expression-phenotype correlations [5, 6, 10, 11]. Herewith, we have detected 76 muscle cis-eQTLs and 28 liver cis-eQTLs which attained genome-wide significance. These numbers are consistent with previous reports. Liaubet et al.  made a genome scan based on microarray measurements of longissimus lumborum gene expression in 57 pigs and identified 335 eQTLs. Of these, only 18 had cis-regulatory effects. Similarly, Cánovas et al.  identified 478 skeletal muscle genome-wide significant eQTLs but only 13% acted in cis-. Although modest sample size (N = 103) limits our ability to identify eQTLs, this sample size ought to allow detecting eQTLs with moderate to large effects on gene expression. As a valid reference, the Genotype-Tissue Expression project detected thousands of eQTLs in 44 tissues characterized by RNA-Seq, and 18 of these tissues were represented by sample sizes equal or below the one used in the current experiment . However, it is also true that sample size largely affects the statistical power of eQTL mapping  and thus, the analysis of a larger number of pigs would most likely uncover additional eQTLs that remained undetected in our study.
Our results indicate that almost 42.8% of the cis-eQTLs regulating gene expression in the liver also display significant associations with the muscle mRNA levels of the very same genes (Tables 1 and 2, Fig. 2). The liver and muscle common cis-eQTLs encompass genes previously related to the function of these two organs, such as the cathepsin F (CTSF) gene, that harbors a polymorphism that has been associated with longissimus dorsi tenderness, ham weight and fatness in Italian crossbred pigs , the transmembrane anterior posterior transformation 1 (TAPT1) gene, which was linked to carcass and eviscerated weight in a GWAS in chicken , the mitochondrial amidoxime reducing component 2 (MARC2) gene, that appears to be involved in porcine fatness , and the transmembrane protein 97 (TMEM97) gene, which contributes to the regulation of cholesterol levels .
There is also a relevant fraction of porcine cis-eQTLs that display significant effects only in one of the two tissues (Figs. 3 and 4). In the Genotype-Tissue Expression (GTEx) pilot experiment , approximately 50% of eQTLs were shared by the nine human tissues under analysis. Moreover, two main types of eQTLs were especially prevalent, i.e. those that regulate gene expression in a single tissue and those that are ubiquitously detected in all tissues. Interestingly, the GTEx pilot analysis also showed that eQTLs affecting gene expression in the skeletal muscle show a limited replicability in other tissues . In other studies, also performed in humans, the eQTL tissue-specificity ranged between 50%  and 60–80% , which implies that the effects of many regulatory mutations are modulated by tissue-associated factors.
Several of the cis-eQTLs detected in our experiment affected the expression of genes involved in lipid or carbohydrate metabolism. For instance, the down-regulation of HTATIP2 leads to elevated fatty acid synthesis and enhanced levels of lipogenic enzymes . The EXOC7 gene regulates the uptake of fatty acids by adipocytes and ACADS is involved in the ß-oxidation of fatty acids , while FAM3C can suppress hepatic gluconeogenesis . It would be interesting to investigate whether polymorphisms associated with the expression of lipid genes also display associations with fatness traits. Two of the muscle cis-eQTLs detected in our study have been previously reported by other authors. An eQTL for the BSDC1 gene was detected by Ponsuksili et al.  and the expression of this gene was also correlated with the percentage of weight loss of the longissimus dorsi muscle. Moreover, a local eQTL that regulates the expression of PNPO and which co-localizes with several meat quality retail traits (such as the percentage of fat and moisture in meat) was described by Steibel et al. . A remarkable level of heterogeneity has been observed in the genetic determinism of production traits in different porcine breeds . In consequence, we anticipated a limited positional concordance amongst eQTLs detected in different breeds. Indeed, a joint analysis of eQTLs across five human populations revealed that varying linkage disequilibrium patterns across populations results in the detection of large numbers of eQTLs with heterogeneous effects .
Limited positional concordance between cis-eQTLs and copy number variation regions
Another goal of our study was to investigate the co-localization of cis-eQTLs and CNVRs in order to make an initial assessment of the potential impact of structural variation on gene expression in pigs. With PennCNV and GADA, we consistently detected 41 CNVRs with a mean size of 457.4 kb. When we compared our CNVR dataset with CNVRs previously reported in pigs [22,23,24,25,26,27,28,29], we found that 60.97% of the CNVRs detected by us had been previously reported in the literature (Additional file 4). A moderate agreement of CNVR locations amongst studies and populations has been evidenced in several reports [27, 49,50,51,52]. Discrepancies could be due to differences in the genetic background of the populations under analysis, filtering criteria (correction factors, criteria to define CNV and CNVRs, etc.), genotyping methods, CNV calling algorithms and the use of family information [27, 53]. In the GM muscle, 2 CNVRs (14 and 38) co-localized with 2 cis-eQTLs (SSC2, 32.6–39.4 Mb and SSC14, 102.7–106.9 Mb), whilst the co-localizations between CNVRs and cis-eQTLs in the liver tissue were somewhat similar, i.e. 2 cis-eQTLs (SSC2, 5.83–5.90 Mb and SSC13, 148.90–149.17 Mb) co-localized with CNVR 11 and 36. This limited positional concordance is probably due to the fact that only 10 CNVRs co-localized with the set of probes analyzed in the current experiment. The low positional concordance between CNVRs and eQTLs could be also due to technical reasons related to the difficulties of detecting CNVs with SNP arrays and the limited sensitivity and poor annotation of porcine microarrays. Besides, the search of CNVs was performed in a population of 350 individuals, while only 103 pigs had available gene-expression data. On the other hand, variations in copy number do not necessarily involve changes in gene expression, e.g. in heterozygous individuals, the loss of one allele can be compensated by an increase in the expression of the other allele, and duplication can generate additional copies which expression is silenced . A future objective would be to investigate the association between the number of copies of genes and their expression levels, but a high throughput method allowing the accurate determination of CNV genotypes will be needed to achieve this goal (genotype determination based on 60 K SNP array data is too imprecise to carry out such analyses).
Despite the fact that there are considerable differences in the gene expression patterns of the porcine liver and skeletal muscle, we have identified a substantial proportion of common cis-eQTLs regulating gene expression in both tissues.
Materials and methods
Phenotyping and genotyping of a commercial Duroc population
We have used a commercial Duroc line of 350 Duroc barrows that were slaughtered at the age of ~ 190 days, with an approximate live weight of 122 kg. This population was generated by crossing 5 boars with ~ 400 sows and selecting one offspring/sow (50 piglets did not complete their productive cycle so we ended up with 350 barrows with valid records). The 350 barrows were bred in the experimental farm of the Pig Control Center of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA). The specific conditions of management and feeding have been previously reported [54, 55]. Pigs were slaughtered in a commercial abattoir following the guidelines of the Spanish Royal Decree 54/1995 (January 20, 1995) to preserve animal welfare. In this way, swine were stunned with carbon dioxide (70% or higher) until they lost consciousness and they were subsequently bled by making an incision in, at least, one carotid artery. At slaughter, GM muscle and liver biopsies were obtained for 103 pigs. Total RNA purification procedures have been previously reported [8, 56]. All experimental procedures were approved by the Ethical Committee of IRTA.
Genomic DNA was extracted from blood samples by following a standard phenol-chloroform protocol. Each pig was genotyped for 62,163 SNPs with the Porcine SNP60 BeadChip (Illumina, SanDiego, CA). The quality of the genotyping results was evaluated with the GenomeStudio software (Illumina). The PLINK software  was used to filter out the SNP markers with minor allele frequency below 5%, missing genotype rate above 10%, and displaying significant departures (P-value < 0.001) from the Hardy-Weinberg equilibrium. Markers which could not be mapped to the Sscrofa11.1 assembly or that mapped to either the X or Y chromosomes were also eliminated from the data set. Markers in complete linkage disequilibrium (r2 > 0.98) were also purged. The final data set used for performing eQTL analyses contained 28,493 SNPs. The filtering criteria used for the CNV analysis was different and only SNPs that did not map to the Sscrofa11.1 assembly or that were located in sex chromosomes were removed, leading to a final number of 46,537 used SNPs.
Differential expression analysis between muscle and liver tissue
Gluteus medius muscle and liver samples were collected from 103 Duroc pigs (Lipgen population) after slaughtering, and immediately frozen in liquid nitrogen. These 103 pigs were selected on the basis of a principal component analysis focused on 13 lipid and growth related traits . We chose individuals representing two different metabolic types, i.e. (i) fat pigs with high intramuscular fat (high saturated and monounsaturated fatty acid content) and also high serum lipid levels, and (ii) pigs that were lean and displayed a low level of intramuscular fat (high polyunsaturated fatty acid content) and circulating lipids . Total RNA was extracted from both GM and liver samples, and mRNA expression profiles were characterized by hybridization to the GeneChip Porcine arrays (Affymetrix Inc., Santa Clara, CA), as previously reported by Cánovas et al. . Hepatic and muscular microarray expression data were deposited in the Gene Expression Omnibus (GEO) public repository, and are accessible through GEO Series accession number GSE115484. The Robust Multi-array Average (RMA) algorithm  was employed for carrying out data pre-processing, background correction, normalization and log-transformation of expression values. Gene Intensity level of significance for detecting expressed probes was calculated with the MAS 5.0 algorithm . Control probes and those probes that did not show expression levels above the detection threshold in all samples were filtered out. Differential expression analysis between the GM muscle and liver followed the guidelines of the limma-trend pipeline [61, 62]. The limma’s empirical Bayes method incorporated a mean-variance trend, thus making possible to adequately model the relationship between variance and gene signal intensity. In our study, fold-change (FC) values reflect the mean expression levels in the GM muscle vs liver. The correspondence between differentially expressed probes (|FC| > 1.5; q-value < 0.05) and genes was based on the Affymetrix porcine annotation data assembled database  and the Biomart database (https://www.ensembl.org/biomart/martview).
Genome scan for expression QTLs
The genome scan for eQTLs targeted a total of 3326 probes that are simultaneously expressed in both tissues. These genes were selected using the following criteria: The MAS 5.0  algorithm was used to detect the probes with expression levels above the Gene Intensity level of significance. Probes were retained when expressed in all the samples analyzed in both tissues. Expressed probes where then annotated according to Biomart database available at Ensembl repositories (https://www.ensembl.org/biomart/martview) and those with no gene correspondence were removed from further analyses. We used the GEMMA software  and followed the methods described by González-Prendes et al.  to carry out the association analyses. The fixed effects and parameters assumed in the statistical model were:
where y is the vector that defines the expression of each gene in the GM muscle and liver of the ith individual; W is the matrix with a column of 1 s and the fixed effects, i.e.“batch of fattening” (with 4 categories) and “laboratory” (microarray data were generated in two different laboratories); α is a c-vector of the corresponding coefficients including the intercept; x is an n-vector of marker genotypes; β is the SNP allelic effect estimated as a regression coefficient on the corresponding x genotype (values − 1, 0, 1); u is a n-vector of random effects with a n-dimensional multivariate normal distribution MVNn (0, λ τ − 1 K) where τ− 1 is the variance of the residual errors; λ is the ratio between the two variance components; K is a known relatedness matrix derived from SNPs and ε is the vector of errors with an MVNn (0, τ − 1 I n) being In the identity matrix. The false discovery rate approach reported by Benjamini and Hochberg  was used to correct for multiple testing at the genome-wide level. In this study, cis-eQTLs were defined as these genomic regions with SNPs significantly associated with a probe located at a maximum distance of ±1 Mb of the associated SNPs. Co-localization was defined as an overlap of at least 1 base pair between the genomic locations of two eQTLs. The GM muscle cis-eQTLs detected with microarrays were validated by considering an RNA-seq data set corresponding to the same tissue and comprising 52 Duroc pigs . These 52 pigs were a subset of the 103 pigs investigated in the current work and the methods to generate GM muscle RNA-seq data are fully described in .
Detection of copy number variation
The PennCNV software  was employed to detect copy number variants (CNVs) on the basis of the information provided by 46,537 autosomal SNPs. The PennCNV software implements a hidden Markov model (HMM) to infer CNV calls for each genotyped sample using as input the intensity signal Log R Ratio (LRR) and the B Allele Frequency (BAF) information generated with the BeadStudio software (Illumina). Samples with a standard deviation of LRR > 0.30 and BAF drift > 0.01 were discarded. Besides, a wave adjustment procedure for genomic waves was carried out . With these filtering steps, 20 samples were eliminated from the data set. Only CNVs spanning three or more consecutive SNPs were taken into account. The CNV calling was performed using the default parameters of the HMM model by assuming a UF factor of 0.01. Copy number variant regions (CNVRs) were created by merging CNVs with an overlap of 80% or more. Additionally we used the Genome Alteration Detection Analysis (GADA) package  to further validate the CNVR detected with PennCNV and minimize the rate of false positives. The GADA software uses a sparse Bayesian algorithm, based on LRR and BAF values obtained along the genome with the BeadStudio software (Illumina), to identify genomic locations where copy number changes occur. Briefly, a LRR average of markers in a chromosome segment is computed and compared to the median of the respective chromosome. Finally, the copy number class is classified as gain, loss or gain/loss. In our study, copy number variations spanning three or more consecutive SNPs were taken into account and the multiple array analysis option was employed. The parameters defined for the Bayesian learning model and the backward elimination were: 0.8 for the sparseness hyperparameter and 8 for the critical value of the backward elimination.
Validation of copy number variant regions by quantitative PCR
Quantitative real time PCR (qPCR) assays were used to validate four CNVRs in 39 porcine samples randomly selected from the Lipgen population. The relative quantification (RQ) of the CNVRs was done as previously described . Primers (Additional file 6) were designed with the Primer Express Software (Applied Biosystems). Copy number variant regions were quantified in 384-well plates using SYBR Select Master Mix in a QuantStudio 12 K Flex Real-Time PCR System platform (Applied Biosystems, Inc., Foster City, CA). Reactions were performed in triplicate, and they contained 7.5 ng of genomic DNA and primers at 300 nM in a final volume of 15 μl. The thermocycling profile was: one cycle at 95 °C for 10 min plus 40 cycles of 15 s at 95 °C and 1 min at 60 °C. Moreover, a melting curve profile (95 °C for 15 s, 60 °C for 15 s and a gradual increase in temperature with a ramp rate of 1% up to 95 °C) was implemented to maximize the specificity of the amplification reactions. Relative expression values were calculated with the Qbase+ software (Biogazelle, Ghent, Belgium) by applying the 2-ΔΔCt method (after verifying that its assumptions were adequately fulfilled) . Relative expression values were calibrated using the arithmetic mean of 3–5 samples showing the lowest number of copies for each specific assay. In the specific case of CNVR32, which encompasses a deletion, the five samples chosen for calibration were those with RQ values around 2. Normalization of the expression data was done by using a previously reported assay based on the glucagon gene .
absolute value of the fold-change
Copy number variant
Copy number variant regions
expression quantitative trait loci
Genome-wide Efficient Mixed-Model Association
Gene Expression Omnibus
genome-wide association analysis
in other words
Single nucleotide polymorphism
Nica AC, Dermitzakis ET. Expression quantitative trait loci: present and future. Philos Trans R Soc Lond Ser B Biol Sci. 2013;368. https://doi.org/10.1098/rstb.2012.0362.
Albert FW, Kruglyak L. The role of regulatory variation in complex traits and disease. Nat Rev Genet. 2015;16:197–212.
Gamazon ER, Nicolae DL, Cox NJ. A study of CNVs as trait-associated polymorphisms and as expression quantitative trait loci. PLoS Genet. 2011;7:e1001292.
Ardlie KG, Deluca DS, Segre AV, Sullivan TJ, Young TR, Gelfand ET, et al. The genotype-tissue expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348:648–60.
Ponsuksili S, Jonas E, Murani E, Phatsara C, Srikanchai T, Walz C, et al. Trait correlated expression combined with expression QTL analysis reveals biological pathways and candidate genes affecting water holding capacity of muscle. BMC Genomics. 2008;9:367.
Wimmers K, Murani E, Ponsuksili S. Functional genomics and genetical genomics approaches towards elucidating networks of genes affecting meat performance in pigs. Brief Funct Genomics. 2010;9:251–8.
Steibel JP, Bates RO, Rosa GJM, Tempelman RJ, Rilington VD, Ragavendran A, et al. Genome-wide linkage analysis of global gene expression in loin muscle tissue identifies candidate genes in pigs. PLoS One. 2011;6:e16766.
Cánovas A, Pena RN, Gallardo D, Ramírez O, Amills M, Quintanilla R. Segregation of regulatory polymorphisms with effects on the gluteus medius transcriptome in a purebred pig population. PLoS One. 2012;7:e35583.
Heidt H, Cinar MU, Uddin MJ, Looft C, Jüngst H, Tesfaye D, et al. A genetical genomics approach reveals new candidates and confirms known candidate genes for drip loss in a porcine resource population. Mamm Genome. 2013;24:416–26.
Ponsuksili S, Murani E, Brand B, Schwerin M, Wimmers K. Integrating expression profiling and whole-genome association for dissection of fat traits in a porcine model. J Lipid Res. 2011;52:668–78.
Chen C, Yang B, Zeng Z, Yang H, Liu C, Ren J, et al. Genetic dissection of blood lipid traits by integrating genome-wide association study and gene expression profiling in a porcine model. BMC Genomics. 2013;14:848.
Muñoz M, Rodríguez MC, Alves E, Folch JM, Ibañez-Escriche N, Silió L, et al. Genome-wide analysis of porcine backfat and intramuscular fat fatty acid composition using high-density genotyping and expression data. BMC Genomics. 2013;14:845.
Ferraz AL, Ojeda A, López-Béjar M, Fernandes LT, Castelló A, Folch JM, et al. Transcriptome architecture across tissues in the pig. BMC Genomics. 2008;9:173.
Yin F, Sharen G, Yuan F, Peng Y, Chen R, Zhou X, et al. TIP30 regulates lipid metabolism in hepatocellular carcinoma by regulating SREBP1 through the Akt/mTOR signaling pathway. Oncogenesis. 2017;6:e347.
Inoue M, Akama T, Jiang Y, Chun TH. The exocyst complex regulates free fatty acid uptake by adipocytes. PLoS One. 2015;10:e0120289.
Ghosh S, Kruger C, Wicks S, Simon J, Kumar KG, Johnson WD, et al. Short chain acyl-CoA dehydrogenase deficiency and short-term high-fat diet perturb mitochondrial energy metabolism and transcriptional control of lipid-handling in liver. Nutr Metab. 2016;13:17.
Farris W, Mansourian S, Chang Y, Lindsley L, Eckman EA, Frosch MP, et al. Insulin-degrading enzyme regulates the levels of insulin, amyloid ß-protein, and the ß-amyloid precursor protein intracellular domain in vivo. Proc Natl Acad Sci U S A. 2003;100:4162–7.
Castellano BM, Thelen AM, Moldavski O, Feltes M, van der Welle RE, Mydock-McGrane L, et al. Lysosomal cholesterol activates mTORC1 via an SLC38A9–Niemann-pick C1 signaling complex. Science. 2017;355:1306–11.
Chen Z, Ding L, Yang W, Wang J, Chen L, Chang Y, et al. Hepatic activation of the FAM3C-HSF1-CaM pathway attenuates hyperglycemia of obese diabetic mice. Diabetes. 2017;66:1185–97.
Skogsberg J, Lundström J, Kovacs A, Nilsson R, Noori P, Maleki S, et al. Transcriptional profiling uncovers a network of cholesterol-responsive atherosclerosis target genes. PLoS Genet. 2008;4:e1000036.
Dubland JA, Francis GA. Lysosomal acid lipase: at the crossroads of normal and atherogenic cholesterol metabolism. Front Cell Dev Biol. 2015;3:3.
Fadista J, Nygaard M, Holm L-E, Thomsen B, Bendixen C. A snapshot of CNVs in the pig genome. PLoS One. 2008;3:e3916.
Fowler KE, Pong-Wong R, Bauer J, Clemente EJ, Reitter CP, Affara NA, et al. Genome wide analysis reveals single nucleotide polymorphisms associated with fatness and putative novel copy number variants in three pig breeds. BMC Genomics. 2013;14:784.
Ramayo-Caldas Y, Castelló A, Pena RN, Alves E, Mercadé A, Souza CA, et al. Copy number variation in the porcine genome inferred from a 60 k SNP BeadChip. BMC Genomics. 2010;11:593.
Schiavo G, Dolezal MA, Scotti E, Bertolini F, Calò DG, Galimberti G, et al. Copy number variants in Italian large white pigs detected using high-density single nucleotide polymorphisms and their association with back fat thickness. Anim Genet. 2014;45:745–9.
Wang Y, Tang Z, Sun Y, Wang H, Wang C, Yu S, et al. Analysis of genome-wide copy number variations in Chinese indigenous and Western pig breeds by 60 K SNP genotyping arrays. PLoS One. 2014;9:e106780.
Fernández AI, Barragán C, Fernández A, Rodríguez MC, Villanueva B. Copy number variants in a highly inbred Iberian porcine strain. Anim Genet. 2014;45:357–66.
Chen C, Qiao R, Wei R, Guo Y, Ai H, Ma J, et al. A comprehensive survey of copy number variation in 18 diverse pig populations and identification of candidate copy number variable genes associated with complex traits. BMC Genomics. 2012;13:733.
Li Y, Mei S, Zhang X, Peng X, Liu G, Tao H, et al. Identification of genome-wide copy number variations among diverse pig breeds by array CGH. BMC Genomics. 2012;13:725.
D’haene B, Vandesompele J, Hellemans J. Accurate and objective copy number profiling using real-time quantitative PCR. Methods. 2010;50:262–70.
Li B, Qing T, Zhu J, Wen Z, Yu Y, Fukumura R, et al. A comprehensive mouse transcriptomic bodymap across 17 tissues by RNA-seq. Sci Rep. 2017;7:4200.
Fu Y, Hong T. BIN1 regulates dynamic t-tubule membrane. Biochim Biophys Acta. 2016;1863:1839–47.
Seymour J, O’Brien EJ. The position of tropomyosin in muscle thin filaments. Nature. 1980;283:680–2.
Anderson DM, Cannavino J, Li H, Anderson KM, Nelson BR, McAnally J, et al. Severe muscle wasting and denervation in mice lacking the RNA-binding protein ZFP106. Proc Natl Acad Sci. U S A 2016;113:E4494–503.
El-Kabbani O, Darmanin C, Chung RPT. Sorbitol dehydrogenase: structure, function and ligand design. Curr Med Chem. 2004;11:465–76.
Abul-Husn NS, Cheng X, Li AH, Xin Y, Schurmann C, Stevis P, et al. A protein-truncating HSD17B13 variant and protection from chronic liver disease. N Engl J Med. 2018;378:1096–106.
Zamanian-Daryoush M, Lindner D, Tallant TC, Wang Z, Buffa J, Klipfell E, et al. The cardioprotective protein apolipoprotein A1 promotes potent anti-tumorigenic effects. J Biol Chem. 2013;288:21237–52.
Liaubet L, Lobjois V, Faraut T, Tircazes A, Benne F, Iannuccelli N, et al. Genetic variability of transcript abundance in pig peri-mortem skeletal muscle: eQTL localized genes involved in stress response, cell death, muscle disorders and metabolism. BMC Genomics. 2011;12:548.
Davoli R, Schivazappa C, Zambonelli P, Braglia S, Rossi A, Virgili R. Association study between single nucleotide polymorphisms in porcine genes and pork quality traits for fresh consumption and processing into Italian dry-cured ham. Meat Sci. 2017;126:73–81.
Liu R, Sun Y, Zhao G, Wang F, Wu D, Zheng M, et al. Genome-wide association study identifies loci and candidate genes for body composition and meat quality traits in Beijing-you chickens. PLoS One. 2013;8:e61172.
Kogelman LJA, Cirera S, Zhernakova DV, Fredholm M, Franke L, Kadarmideen HN. Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA sequencing in a porcine model. BMC Med Genet. 2014;7:57.
Bartz F, Kern L, Erz D, Zhu M, Gilbert D, Meinhof T, et al. Identification of cholesterol-regulating genes by targeted RNAi screening. Cell Metab. 2009;10:63–75.
Emilsson V, Thorleifsson G, Zhang B, Leonardson AS, Zink F, Zhu J, et al. Genetics of gene expression and its effect on disease. Nature. 2008;452:423–8.
Dimas AS, Deutsch S, Stranger BE, Montgomery SB, Borel C, Attar-Cohen H, et al. Common regulatory variation impacts gene expression in a cell type-dependent manner. Science. 2009;325:1246–50.
Chen Z, Wang J, Yang W, Chen J, Meng Y, Feng B, et al. FAM3C activates HSF1 to suppress hepatic gluconeogenesis and attenuate hyperglycemia of type 1 diabetic mice. Oncotarget. 2017;8:106038–49.
Ponsuksili S, Murani E, Schwerin M, Schellander K, Wimmers K. Identification of expression QTL (eQTL) of genes expressed in porcine M. longissimus dorsi and associated with meat quality traits. BMC Genomics. 2010;11:572.
Zhang W, Yang B, Zhang J, Cui L, Ma J, Chen C, et al. Genome-wide association studies for fatty acid metabolic traits in five divergent pig populations. Sci Rep. 2016;6:24718.
Wen X, Luca F, Pique-Regi R. Cross-population joint analysis of eQTLs: fine mapping and functional annotation. PLoS Genet. 2015;11:e1005176.
Paudel Y, Madsen O, Megens H-J, Frantz LAF, Bosse M, Bastiaansen JWM, et al. Evolutionary dynamics of copy number variation in pig genomes in the context of adaptation and domestication. BMC Genomics. 2013;14:449.
Xie J, Li R, Li S, Ran X, Wang J, Jiang J, et al. Identification of copy number variations in Xiang and Kele pigs. PLoS One. 2016;11:e0148565.
Wiedmann RT, Nonneman DJ, Rohrer GA. Genome-wide copy number variations using SNP genotyping in a mixed breed swine population. PLoS One. 2015;10:e0133529.
Paudel Y, Madsen O, Megens H-J, Frantz LAF, Bosse M, Crooijmans RPMA, et al. Copy number variation in the speciation of pigs: a possible prominent role for olfactory receptors. BMC Genomics. 2015;16:330.
Clop A, Vidal O, Amills M. Copy number variation in the genomes of domestic animals. Anim Genet. 2012;43:503–17.
Gallardo D, Pena RN, Amills M, Varona L, Ramírez O, Reixach J, et al. Mapping of quantitative trait loci for cholesterol, LDL, HDL, and triglyceride serum concentrations in pigs. Physiol Genomics. 2008;35:199–209.
Gallardo D, Quintanilla R, Varona L, Díaz I, Ramírez O, Pena RN, et al. Polymorphism of the pig acetyl-coenzyme a carboxylase α gene is associated with fatty acid composition in a Duroc commercial line. Anim Genet. 2009;40:410–7.
González-Prendes R, Quintanilla R, Cánovas A, Manunza A, Figueiredo-Cardoso T, et al. Joint QTL mapping and gene expression analysis identify positional candidate genes influencing pork quality traits. Sci Rep. 2017;7:39830.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.
Cánovas A, Quintanilla R, Amills M, Pena RN. Muscle transcriptomic profiles in pigs with divergent phenotypes for fatness traits. BMC Genomics. 2010;11:372.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003;4:249–64.
Gautier L, Cope L, Bolstad BM, Irizarry RA. Affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004;20:307–15.
Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, Medvedovic M. Intensity-based hierarchical Bayes method improves testing for differentially expressed genes in microarray experiments. BMC Bioinformatics. 2006;7:538.
Law CW, Chen Y, Shi W, Smyth GK. Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29.
Carlson M. Affymetrix porcine annotation data (chip porcine). R package version 3.2.3; 2016.
Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44:821–4.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.
Cardoso TF, Cánovas A, Canela-Xandri O, González-Prendes R, Amills M, Quintanilla R. RNA-seq based detection of differentially expressed genes in the skeletal muscle of Duroc pigs with distinct lipid profiles. Sci Rep. 2017;7:40005.
Wang K, Li M, Hadley D, Liu R, Glessner J, Grant SFA, et al. PennCNV: an integrated hidden Markov model designed for high-resolution copy number variation detection in whole-genome SNP genotyping data. Genome Res. 2007;17:1665–74.
Pique-Regi R, Monso-Varona J, Ortega A, Seeger RC, Triche TJ, Asgharzadeh S. Sparse representation and Bayesian detection of genome copy number alterations from microarray data. Bioinformatics. 2008;24:309–18.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCt method. Methods. 2001;25:402–8.
Revilla M, Puig-Oliveras A, Castelló A, Crespo-Piazuelo D, Paludo E, Fernández AI, et al. A global analysis of CNVs in swine using whole genome sequence data and association analysis with fatty acid composition and growth traits. PLoS One. 2017;12:e0177014.
The authors are indebted to Selección Batallé S.A. for providing the animal material. We gratefully acknowledge to J. Reixach (Selección Batallé) and J. Soler (IRTA) for their collaboration in the experimental protocols.
Part of the research presented in this publication was funded by grants AGL2013–48742-C2–1-R and AGL2013–48742-C2–2-R awarded by the Spanish Ministry of Economy and Competitivity. We also acknowledge the support of the Spanish Ministry of Economy and Competitivity for the Center of Excellence Severo Ochoa 2016–2019 (SEV-2015-0533) grant awarded to the Centre for Research in Agricultural Genomics. Ángela Cánovas was funded under the Juan de la Cierva program awarded by the Spanish Ministry of Science. Taina F Cardoso was funded with a fellowship from the CAPES Foundation-Coordination of Improvement of Higher Education, Ministry of Education (MEC) of the Federal Government of Brazil. Rayner González-Prendes was funded by a FPU Ph.D. grant from the Spanish Ministry of Education (FPU12/00860). Emilio Mármol-Sánchez was funded by a FPU Ph.D. grant from the Spanish Ministry of Education (FPU15/01733).
Animal care, management procedures and blood sampling were performed following national guidelines for the Good Experimental Practices and they were approved by the Ethical Committee of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. List of genes that are differentially expressed (DE) in the gluteus medius muscle and the liver (q-value < 0.05). (XLSX 353 kb)
Figure S1. Boxplots depicting the mRNA expression levels of cis-eQTL regulated genes measured with RNA-Seq and microarrays in the gluteus medius muscle of 52 and 103 Duroc pigs, respectively. Means were compared with a Student’s t- test: P-value > 0.05 (ns); P-value ≤ 0.05 (*); P-value ≤ 0.01 (**); P-value ≤ 0.001 (***) and P-value ≤ 0.0001 (****). (DOCX 807 kb)
Table S2. Locations, absolute frequencies, lengths and status of CNVRs detected in 350 Duroc pigs with the PennCNV and GADA softwares. (XLSX 12 kb)
Table S3. Positional coincidence of the copy number variant regions (CNVR) detected in the current work with CNVRs reported in other pig populations. (XLSX 10 kb)
Table S4. Co-localization of CNVRs and cis-eQTLs in the gluteus medius muscle (eQTL and gene positions are expressed in Mb). Table S5. Co-localization of CNVRs and cis-eQTLs in liver tissue (eQTL and gene positions are expressed in Mb). (XLSX 13 kb)
Table S6. Sequences of the primers used to validate CNVRs by qPCR. (XLSX 9 kb)