Abstract
Iberian ham production includes both purebred (IB) and Duroc-crossbred
(IBxDU) Iberian pigs, which show important differences in meat quality
and production traits, such as muscle growth and fatness. This
experiment was conducted to investigate gene expression differences,
transcriptional regulation and genetic polymorphisms that could be
associated with the observed phenotypic differences between IB and
IBxDU pigs. Nine IB and 10 IBxDU pigs were slaughtered at birth.
Morphometric measures and blood samples were obtained and samples from
Biceps femoris muscle were employed for compositional and transcriptome
analysis by RNA-Seq technology. Phenotypic differences were evident at
this early age, including greater body size and weight in IBxDU and
greater Biceps femoris intramuscular fat and plasma cholesterol content
in IB newborns. We detected 149 differentially expressed genes between
IB and IBxDU neonates (p < 0.01 and Fold-Change > 1. 5). Several were
related to adipose and muscle tissues development (DLK1, FGF21 or UBC).
The functional interpretation of the transcriptomic differences
revealed enrichment of functions and pathways related to lipid
metabolism in IB and to cellular and muscle growth in IBxDU pigs.
Protein catabolism, cholesterol biosynthesis and immune system were
functions enriched in both genotypes. We identified transcription
factors potentially affecting the observed gene expression differences.
Some of them have known functions on adipogenesis (CEBPA, EGRs), lipid
metabolism (PPARGC1B) and myogenesis (FOXOs, MEF2D, MYOD1), which
suggest a key role in the meat quality differences existing between IB
and IBxDU hams. We also identified several polymorphisms showing
differential segregation between IB and IBxDU pigs. Among them,
non-synonymous variants were detected in several transcription factors
as PPARGC1B and TRIM63 genes, which could be associated to altered gene
function. Taken together, these results provide information about
candidate genes, metabolic pathways and genetic polymorphisms
potentially involved in phenotypic differences between IB and IBxDU
pigs associated to meat quality and production traits.
Introduction
The pig is the main species for meat consumption worldwide, 43% of
total produced meat comes from pigs. Most production comes from the
modern European pig breeds, which have been extensively selected and
show optimized productivity and efficiency [[47]1]. In the
Mediterranean basin, there is also a significant production of unique
high-quality traditional pork products from local breeds. The
Mediterranean breeds, also known as fatty-pig breeds, have an ancient
origin, and have been reared in extensive conditions for centuries,
exposed therefore to harsh environments and seasonal variations in food
availability (associated with the development of a thrifty genotype
[[48]2]). These breeds are smaller in size, have not undergone intense
genetic selection and are less productive than modern breeds. As a
consequence of the industrialization of pork production, three-quarters
of the traditional breeds are extinct or marginalized [[49]3]. The
exception is the Iberian pig, the most representative Mediterranean
traditional breed, which has an important commercial value based on
high quality dry-cured products in terms of consumers’ health and
acceptance [[50]4].
Peculiarities in Iberian pig metabolism drive its valued meat
properties; Iberian pigs are characterized by higher fat deposition,
fat desaturation and food intake [[51]5, [52]6], as well as by higher
circulating leptin levels in plasma [[53]7] than lean pigs, suggesting
a syndrome of leptin resistance. Moreover, the Iberian pig is also
considered an amenable and robust biomedical model for obesity and
associated cardiometabolic diseases since, when provided high levels of
food, the animals are prone to the development of dyslipidemias,
metabolic syndrome and type-2 diabetes [[54]8]. On the other hand, as
observed in other traditional breeds, productive performance is
considerably lower than that of highly selected modern breeds. To
improve reproductive and growth performances and primal cuts yield, in
the last decades Duroc breed was introduced as terminal sire cross.
Recently, Spanish law has accepted and regulated the use of Iberian X
Duroc pigs to obtain “Iberian” products.
However, the introduction of Duroc genetics is associated with a
decrease in meat quality, mainly determined by a decrease in
intramuscular fat (IMF) and monounsaturated fatty acids (MUFA) contents
[[55]9]. Intramuscular fat content and fatty acid composition are the
main factors affecting meat quality and are highly dependent on genetic
type and diet [[56]10]. Intramuscular fat content is determined both by
number and size of adipocytes within muscle fibers. During prenatal
development and immediately after birth, preadipocyte differentiation
is a very active process that slow down with animal growth [[57]11].
Later in growth adipocyte hypertrophy is the most important issue
affecting IMF content, although hyperplasia is maintained in the adult
animal to a lesser extent [[58]12]. Thus, birth is a critical
time-point to investigate the adipocyte differentiation process. On the
other hand, IMF composition and fatty acids profile depend on
lipogenesis and fatty acids metabolism. It has been reported that breed
affects adipogenesis, lipogenesis and their timing, as well as the
expression patterns of adipocyte differentiation-related genes
[[59]13]. In this sense, Iberian pig is considered a more precocious
breed than Duroc pig [[60]14].
Due to the influence of the genetic background on productive and meat
quality traits, research in the past few decades has been focused on
understanding the genetic basis of cell growth and development,
myogenesis and metabolism [[61]15]. Recently, new interest has arisen
towards the understanding of genetic mechanisms underlying lipid
synthesis and accumulation, due to its importance in meat quality
[[62]13]. Different approaches such as candidate gene expression
studies or cDNA microarray analysis have been used to investigate
genetic aspects of target parameters. Some studies based on the
microarray technology investigated transcriptome differences among
Iberian pig and Large White or Duroc pig in endocrine tissues [[63]16]
and between Iberian and Iberian X Duroc crossbred pigs in Longissimus
dorsi muscle [[64]14].
Currently, the availability of the RNA-Seq technology has allowed the
assessment of global changes in transcriptome of a number of species
including pigs [[65]15], because of its greater accuracy and
reproducibility than microarray technology [[66]17, [67]18]. RNA-Seq
allows measuring not only gene expression, but also examining genome
structure identifying SNP and other structural variation such as indel
and splice variants. Some applications of this technology include
transcript quantification, allele-specific expression, novel transcript
discovery or single nucleotide polymorphism (SNP) discovery [[68]19].
In pigs, several RNA-Seq studies have been carried out for assessing
differences in the transcriptome of muscle, fat, liver or hypothalamus
among breeds or phenotypically extreme individuals within a breed for
characters of interest [[69]15, [70]20, [71]21]. The RNA-Seq technology
is still scarcely applied to the Iberian breed, with studies
comprehending the assessment of phenotypically extreme individuals for
fatty acids composition [[72]22] or the exploration of gonad
transcriptome in Iberian and Large White pigs [[73]23]. However, to the
best of our knowledge, there are not RNA-Seq technology-based studies
focused on genetic differences between Iberian and other breeds aimed
at improving meat quality and productive traits.
Meat quality in Iberian pigs is of special interest for carcass cuts
used in the dry curing industry such as the loin and the ham. A
previous study assessed transcriptomic differences between pure Iberian
and Duroc-crossbred Iberian pigs using microarray technology in the
loin [[74]14], but no information on ham muscles transcriptome exists.
It is well known that different muscles differ in developmental timing,
metabolic and physicochemical properties, including different responses
to exercise [[75]24]. Biceps femoris (BF) muscle is the biggest muscle
in the ham and shows higher oxidative capacity, and lower drip loss
than Longissimus dorsi muscle [[76]25, [77]26]. Moreover, important
differences exist regarding the IMF content of both muscles. Karlsson
et al. (1993) reported higher IMF content in LD muscle in Yorkshire pig
breed, whilst the opposite was reported for Iberian pigs, where BF
showed remarkably greater IMF content than several others carcass
muscles [[78]27]. Also, transcriptomic and proteomic comparisons
between muscles showed important functional differences, with 15–30% of
proteome differing between LD and BF [[79]24, [80]28].On the other
hand, transcriptomic studies performed sequentially along early
development suggest the perinatal as a critical period to study genes
affecting muscle and adipose cells growth and muscle fiber
differentiation [[81]20, [82]29], in agreement with the tissue
differentiation timing commented previously. Moreover, environmental
effects are minimized at this time point.
Hence, in agreement with previous considerations, the present study was
carried out to study the BF muscle of newborn piglets in IB and in the
IBxDU cross, aiming to: 1) Verify whether phenotypic differences are
evident from the very early developmental stages (newborns) in these
closely related populations; 2)Evaluate changes in gene expression in
BF muscle that may be responsible for the observed phenotypic
differences and identify pathways and networks in which those genes are
involved; 3) Identify transcription factors affecting gene expression
in order to establish potential new candidate genes affecting
productive parameters and meat quality; 4) Identify structural variants
in these candidate genes, potentially involved in the observed
expression differences.
These results are useful for the understanding of genetic pathways
affecting pork production and may be also of translational value for
the understanding of ethnic differences in obesity and associated
disorders in lipid metabolism in human medicine.
Materials and Methods
Ethics statement
Animal manipulations were done in compliance with the regulations of
the Spanish Policy for Animal Protection RD1201/05, which meets the
European Union Directive 86/609 about the protection of animals used in
research. The experiment was specifically assessed and approved (report
CEEA 2010/003) by the INIA Committee of Ethics in Animal Research,
which is the named Institutional Animal Care and Use Committee (IACUC)
for the INIA.
Animals and sample collection
Ten pure Iberian sows raised in the same commercial farm were employed
at their third gestation cycle. All females were managed in the same
conditions. Five sows were mated to Iberian boars and five to Duroc
boars. At birth, nine pure Iberian (IB) and 10 Iberian x Duroc (IBxDU)
male piglets were randomly selected from the ten litters (two from each
litter excepting one litter providing just one Iberian male). Blood
samples were collected from newborns in sterile heparin blood vacuum
tubes (Vacutainer Systems Europe, Meylan, France). Immediately after
recovery, the blood was centrifuged at 1500g for 15 min and the plasma
was separated and stored into polypropylene vials at −20°C until
assayed for determination of glucose and lipids metabolism-indicating
parameters. After blood collection, piglets were slaughtered. Several
body development measures were obtained with a measure-tape: total body
length (from the rostral edge of the snout to the tail insertion), ham
length (from the anterior edge of the Symphysis pubica to the
articulatio tarsi), total length of anterior and posterior limbs (from
the distal edge of the hooves to the proximal edge of the scapula or
Symphysis pubica, respectively) and thoracic, abdominal and ham
circumferences. Carcasses were weighted and samples from BF muscle were
vacuum-packed in low-oxygen permeable film and kept frozen at –20°C
until fatty acid composition analysis. Prior to fatty acid analysis,
muscle samples were freeze dried for two days in a lyophilizer
(Lyoquest, Telstar, Tarrasa, Spain) and grounded in a Mixer Mill MM400
(Retsch technology, Haan, Germany) until muscle was completely
powdered. For transcriptomic analysis, BF samples were immediately
frozen in liquid nitrogen and maintained at –80°C until RNA extraction.
The metabolic status of the newborn piglets was evaluated. Glucose,
fructosamine, triglycerides, total cholesterol, high-density
lipoprotein cholesterol (HDL-c) and low-density lipoprotein cholesterol
(LDL-c) plasmatic levels were measured with a clinical chemistry
analyzer (Saturno 300 plus, Crony Instruments s. r. l., Rome, Italy).
Tissue composition analysis
Biceps femoris muscle IMF content was quantified using the method
proposed by Segura and López-Bote [[83]30] based on gravimetrical
determination of lipid content. Fatty acid methyl esters (FAMEs) were
identified by gas chromatography as described by López-Bote et al
[[84]31] using a Hewlett Packard HP-6890 (Avondale, PA, USA) gas
chromatograph equipped with a flame ionization detector and a capillary
column (HP-Innowax, 30 m × 0.32 mm i.d. and 0.25 μm polyethylene
glycol-film thickness). Results were expressed as grams per 100 grams
of detected FAMEs.
Transcriptomic analysis
RNA extraction
A total of 12 animals were randomly selected to perform transcriptomic
analysis, representing all available litters (6 animals of each genetic
type). Total RNA was extracted from 50–100mg samples of BF muscle using
the RiboPure TM of High Quality total RNA kit (Ambion, Austin, TX, USA)
following the manufacturer’s recommendations. RNA was quantified using
a NanoDrop-100 spectrophotometer (NanoDrop Technologies, Wilmington,
DE, USA). The quality of the RNA was evaluated using the RNA Integrity
Number (RIN) value from the Agilent 2100 Bioanalyzer device (Agilent
technologies, Santa Clara, CA, USA). The RIN values ranged from 7.5 to
9.8.
Library construction and RNA sequencing
Sequencing libraries were made using the mRNA-Seq sample preparation
kit (Illumina Inc., Cat. # RS-100-0801) according to manufacturer’s
protocol. Each library was sequenced using TruSeq SBS Kit v3-HS, in
paired end mode with the read length 2x76bp on a HiSeq2000 sequence
analyzer (Illumina, Inc). Images from the instrument were processed
using the manufacturer’s software to generate FASTQ sequence files.
Mapping and assembly
Sequence reads were analyzed using CLC Bio Genomic workbench software
7.0 (CLC Bio, Aarhus, Denmark). Quality control analysis was performed
using the NGS quality control tool, which assesses sequence quality
indicators based on the FastQC-project
([85]http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).
Quality was measured taking into account sequence-read lengths and
base-coverage, nucleotide contributions and base ambiguities, quality
scores as emitted by the base caller and over-represented sequences
[[86]32]. All the samples analyzed passed all the QC parameters having
the same length (76 bp), 100% coverage in all bases, 25% of A, T, G and
C nucleotide contributions, 50% GC on base content and less than 0.1%
over-represented sequences. A hierarchical clustering of the samples
was also performed. One IBxDU pig was discarded for further analysis
because the sample deviated largely from the expected grouping in the
clustering analysis, probably due to RNA sampling or processing
problems. Sequence paired-end reads (76bp) were assembled against the
annotated Sscrofa10.2 reference genome
([87]http://www.ncbi.nlm.nih.gov/genome/?term=sus+scrofa) using the
genome, annotated genes and mRNA tracks. Data was normalized by
calculating the ‘fragments per kilo base per million mapped reads’
(FPKM) for each gene [[88]33].
Differential expression analysis
The statistical analysis was performed using the total exon reads as
expression values by the Empirical analysis of differential gene
expression tool. This tool is based on the EdgeR Bioconductor package
[[89]34] and uses count data (i.e. total exon reads) for the
statistical analysis. Genes were filtered according to two criteria: a
minimum mean group expression greater than 0.5 FPKM in at least one
group and a Fold-Change (FC) of the expression differences between IB
and IBxDU groups equal or higher to 1.5. Finally, those genes with a p
≤ 0.01, corresponding to a false discovery rate (FDR) value ≤ 0.23,
were considered as differentially expressed (DE).
Systems biology study
The biological interpretation of the DE genes observed in BF muscle was
performed using three complementary approaches, in order to identify
enriched GO terms, pathways and networks involving the DE genes, and
potential regulators causing the observed changes in gene expression.
The enrichment analysis was carried out using the Wilcoxon test tool in
the ConsensusPathDB database [[90]35], available at the Max Plank
Institute website, which provides batch enrichment analyses to
highlight the most relevant GO terms associated to a gene list. Both,
the list of genes overexpressed in IB and IBxDU were used. Functional
terms with p-values lower than 0.05 were considered enriched in the
annotation categories. The p-values are corrected for multiple testing
using the FDR method and are presented as q-values.
Additionally, Ingenuity Pathway Analysis, (IPA)(Ingenuity Systems,
Qiagen, California)software was employed to identify and characterize
biological functions, gene networks and canonical pathways affected by
the DE genes.
Regulatory transcription factors (TRF), which could potentially affect
the DE genes in the dataset were also studied by following
complementary approaches. First, Regulatory Impact Factors (RIF1 and
RIF2) metrics [[91]36, [92]37] were calculated for the whole set of DE
genes obtained conditional on genetic type (149 genes). RIF1 assigns an
extreme score to those TRF that are consistently most differentially
co-expressed with the highly abundant and highly DE genes; RIF2 assigns
an extreme score to those TRF with the most altered ability to act as
predictors of the abundance of DE genes. Candidate TRFs in pigs were
obtained from Animal TFDB
([93]http://www.bioguo.org/AnimalTFDB/BrowseAllTF.php?spe=Susscrofa). A
total of 1038 TRF were retrieved. Among them, 723 showed expression
values greater than 0.5 FPKM in at least one experimental group and
thus, were used in the RIFs metrics approach.
The RIF1 and RIF2 values were computed for the i ^th TRF as follows:
[MATH:
RIF1
i=1<
msub>nde∑j=1j=ndea^jxd^j(r1ij−r
2ij)2and :MATH]
[MATH:
RIF2
i=1<
msub>nde∑j=1j=n<
mi>de[(e1jrx1<
/mn>ij)2−<
mrow>(e2j<
/mrow>xr2ij)2
] :MATH]
where n [de] is the number of DE genes, a [j] and d [j] the estimated
average expression and differential expression of the j ^th DE gene, r1
[ij] and r2 [ij] the co-expression correlation between the i ^th TRF
and the j ^th DE gene in each one of the genetic types and beinge1 [j]
and e2 [j] the expression of the j ^th gene in each genetic type
[[94]38]. Both RIF measures for each analyzed TRF were transformed to
standardized z-scores by subtracting the mean and dividing by its
standard deviation. We identified relevant TRF as those with extreme
RIF z-scores according to the corresponding confidence intervals (CI)
calculated by bootstrap. In each iteration of bootstrapping, a set of
n[de] = 149 genes were randomly selected from the 11392 expressed
genes, and the RIF1 and RIF2 z-scores of the 723 TRF were calculated.
The procedure was repeated 10,000 times for each scenario to obtain the
corresponding 95 and 99% CI intervals of both z-scores.
Complementarily, IPA software was employed to identify and characterize
potential regulators using two different tools, the upstream regulators
and the regulators tools. Both of them identify known regulators that
may be affecting expression of the dataset of DE genes. IPA-identified
regulators include genes, but also other molecules as drugs. Thus, out
of the identified regulators, only genes that were also included in the
RIFs metrics candidate TRF list (which consisted of 723 TRF) were
considered (genes included in the animal TFDB and with expression
values higher than 0.5 FPKM in at least one experimental group). Using
the information obtained from the TRF study, an additional search for
enriched pathways and networks was carried out with IPA software
considering both, DE genes and TRF.
Structural variants analysis
A search of structural variants was performed by pooling the reads
coming from all animals in each genetic type, and comparing the
variants found in each group. The probabilistic variants detection tool
(CLC Bio Genomic workbench) was used to perform the variant calling
analysis. Single-end read alignments were not ignored. The minimum
coverage in a locus to be considered was set up as 10 and the variant
probability as 90. The variant probability parameter defines how good
the evidence has to be at a particular site for the tool to report a
variant at that location. Specifically, the variant probability
threshold set as 90 means that any candidate variation in the genome
must have a probability lower than 0.1 (1–0.9) of being the same as the
reference sequence, to be considered as a variant. Variants with an
allele frequency under 5% and/or coverage under 30 reads for the
general variant analysis and under 10 reads in the candidate genes
variant analysis were not considered. Variants were considered to be
potentially fixed (frequency greater or equal to 90%) or segregating
(frequency lower than 90%).
The variants identified in genes corresponding to transcription factors
(i.e. EGR2, FOS, FOXO1, FOXO3, IRF1, STAT5B, HOXA9, ATF4, TP53, NOR-1,
ABRA, ATF3 and PPARGC1B) were functionally evaluated using the variant
effect predictor (VEP) tool from Ensembl
([95]www.ensembl.org/info/docs/tools/vep/) which includes information
about amino acid change localization and consequences, affected
transcripts, and SIFT [[96]39] and PolyPhen [[97]40] scores.
Results validation by RT qPCR
RNA obtained from the 11 animals under study was employed to perform
the technical validation of the differential expression of some genes
that were either upregulated in IB, upregulated in IBxDU or not DE
between genetic types. This technical validation was performed by
studying the Pearson correlation between the expression values obtained
from RNA-Seq data (FPKM) and the normalized gene expression data
obtained by RT qPCR. Moreover, RNA obtained from all the available
animals (9 IB and 10 IBxDU) was used to quantify expression differences
of such genes.
First-strand cDNA synthesis was carried out with Superscript II
(Invitrogen, Life Technologies, Paisley, UK) and random hexamers in a
total volume of 20 μl containing 1 μg of total RNA and following the
supplier’s instructions.
The expression of 9 genes was quantified by qPCR. Primer pairs used for
quantification were designed using Primer Select software (DNASTAR,
Wisconsin, USA) from the available GENBANK and/or ENSEMBL sequences,
covering different exons in order to assure the amplification of the
cDNA. Sequence of primers and amplicon lengths are indicated in [98]S1
Table. Standard PCRs on cDNA were carried out to verify amplicon sizes.
Quantification was performed using SYBR Green mix (Roche, Basel,
Switzerland) in a LightCycler480 (Roche, Basel, Switzerland), following
standard procedures Data were analyzed with LightCycler480 SW1.5
software (Roche, Basel, Switzerland). All samples were run in
triplicate and dissociation curves were carried out for each individual
replicate. Single peaks in the dissociation curves confirmed the
specific amplification of the genes. For each gene, PCR efficiency was
estimated by standard curve calculation using four points of cDNA
serial dilutions. Mean Cp values were employed for the statistical
analyses of differential expression. Stability of four endogenous genes
(i. e. GAPDH, B2M, TBP and ACTB) was calculated using Genorm software
[[99]41] The TBP and ACTB genes were selected as the most stable
endogenous genes to normalize the data. The qPCR expression data
normalization was performed using normalization factors calculated with
Genorm software. Relative quantities were divided by the normalization
factors, which were the geometric means of the two reference genes
quantities.
Statistical analyses of tissue composition and qPCR expression quantification
Phenotypic data were analyzed as a completely randomized design using
the general linear model (GLM) procedure using SAS version 9.2 (SAS
Inst. Inc., Cary, NC; 2009). The mean and genetic type were considered
as systematic effects, and residual effects as random. Carcass weight
was used as covariate when it was significant and removed from the
model when it was insignificant. The animal was the experimental unit
for all analysis. The results were considered to be significant at
p-value<0.05.
Statistical analysis of gene expression data was carried out following
the method proposed by Steibel et al [[100]42] which consists of the
analysis of cycles to threshold values (Cp), for the target and
endogenous genes using a linear mixed model. The following model was
used for analyzing the joint expression of the target and control genes
in different tissues:
[MATH:
ygik
r=TGgi+Bgk+Dik+egikr :MATH]
where,
[MATH:
ygij
kr=−lo<
/mi>g2(Eg−Cpgij
kr) :MATH]
, E [g] is the efficiency of the PCR of gth gene, Cp [gikr] is the
value obtained from the thermocycler software for the gth gene from the
rth replicate in a sample collected from the kth animal of the ith
genetic type, TG [gi] is the specific effect of the ith genetic type on
the expression of gene gth, B [gjk] is specific random effect of the
kth pig on the expression of gene gth, D [ijk] is a random
sample-specific effect common to all the genes, and e [gikr] is a
residual effect.
To test differences in the expression rate of genes of interest (diff
[TG]) between classes normalized by the endogenous genes, different
contrasts were performed between the respective estimates of TG levels.
Significance of diff [TG] estimates was determined with the t
statistic. To obtain FC values from the estimated diff [TG] values, the
following equation was applied:
[MATH: FC=2<
mo>−diffTG :MATH]
P-values<0.05 were considered statistically significant.
To validate the global RNA-Seq results, the concordance correlation
coefficient (CCC) [[101]43] was calculated between the FC values
estimated in BF muscle from RNA-Seq and qPCR expression measures for
the 9 genes analyzed by the two technologies (RNA-Seq and qPCR).
Results and Discussion
Phenotypic differences between genetic types
The results obtained in the present study constitute the first
assessment of phenotypic differences between IB and IBxDU piglets at
birth. There are several studies evaluating phenotypic differences
between both genotypes at weaning or adulthood [[102]9, [103]14,
[104]44–[105]46]. Pure Iberian and crossbred piglets were slaughtered
at birth at an average of 1.2 and 1.8 kg live weight, respectively (SEM
= 0.06). Genetic type affected all the carcass phenotypic parameters:
IBxDU neonates were bigger and heavier (p < 0.001) than IB newborns
([106]Table 1), reflecting previously reported differences in the same
traits in adult animals [[107]45, [108]46]. Birth is an interesting
sampling time to study differences in metabolism and growth, in which
environmental influences are minimized, but maternal effects may have
to be considered. Unfortunately no information regarding the employed
sows’ phenotype was available, but all of them originated from the same
commercial Iberian population and were checked for breed purity by
genotyping several markers usually employed for breed origin
authentication, showing homocygosity for Iberian alleles. Moreover age
and productive cycle were also homogeneous among them as indicated in
the methods section, and thus, we don’t expect that potential sow
effects may interfere with the observed between-genotype differences
regarding piglet’s phenotype or transcriptome.
Table 1. Carcass, Biceps Femoris and metabolism phenotypic characteristics in
IB and IBxDU piglets.
Genetic type SEM[109] ^3 p-value
IBxDU[110] ^1 IB[111] ^2
Carcass characteristics
Carcass weight, kg 1.41 0.96 0.05 0.0005
Ham weight, kg 0.16 0.11 0.00 0.0008
Total body lenght, cm 40.20 35.50 0.31 0.0004
Ham lenght, cm 7.45 6.33 0.12 0.0009
Forelimb lenght, cm 12.35 10.83 0.12 0.0042
Hind limb lenght, cm 15.95 13.83 0.12 0.0016
Torax circumference, cm 25.15 22.06 0.14 0.0010
Abdomen circumference, cm 18.90 17.28 0.19 0.0486
Ham circumference, cm 12.55 10.89 0.15 0.0020
Lipid and glucose metabolism-related plasma indicators
Cholesterol, mg/dl 62.19 102.36 5.60 0.0030
Fructosamine, mg/dl 169.70 133.67 10.37 0.1009
Glucose, mg/dl 132.40 123.44 10.80 0.6839
LDL[112] ^4 , mg/dl 42.16 45.82 4.40 0.4496
HDL[113] ^5 , mg/dl 22.38 41.20 4.25 0.0176
Triglycerides, mg/dl 30.00 76.67 5.11 0.0003
Biceps femoris muscle fatty acids composition (g/100 g total fatty
acids)
IMF[114] ^6 , % 1.72 2.21 0.09 0.0142
C12:0 0.65 0.58 0.03 0.2321
C14:0 2.57 2.32 0.12 0.3189
C15:1 1.28 1.18 0.06 0.3762
C16:0 25.90 25.44 0.19 0.2379
C16:1 n-9 1.90 2.09 0.05 0.3854
C16:1 n-7 5.38 4.57 0.21 0.0773
C17:0 1.69 1.44 0.07 0.0814
C17:1 0.91 0.81 0.06 0.3858
C18:0 10.85 9.96 0.25 0.1014
C18:1 n-9 23.80 25.82 0.56 0.0921
C18:1 n-7 6.15 5.69 0.18 0.2163
C18:2 n-6 7.31 9.17 0.60 0.1395
C20:1 n-9 0.53 0.53 0.01 0.9701
C20:2 n-6 0.41 0.40 0.03 0.8331
C20:3 n-6 0.62 0.55 0.02 0.0293
C20:4 n-6 6.31 5.99 0.20 0.4264
C22:1 n-9 1.21 1.08 0.05 0.1956
C22:4 n-6 1.58 1.27 0.09 0.0925
C22:5 n-3 0.48 0.50 0.02 0.5405
C22:6 n-3 0.67 0.62 0.02 0.1959
∑SFA[115] ^7 41.66 39.74 0.42 0.0350
∑MUFA[116] ^8 41.01 41.77 0.40 0.3530
∑PUFA[117] ^9 17.34 18.49 0.50 0.2639
UI[118] ^10 96.20 97.79 0.90 0.3904
∑n-3[119] ^11 1.77 1.67 0.04 0.1940
∑n-6[120] ^12 15.56 16.82 0.48 0.2089
∑n-6/∑n-3 8.78 10.17 0.30 0.0319
∑MUFA/∑SFA 0.99 1.05 0.02 0.0659
[121]Open in a new tab
^1 IBxDU = Iberian x Duroc crossbred pigs (n = 10)
^2 IB = Purebred Iberian pigs (n = 9)
^3SEM = Standard error of the mean
^4LDL = Low density lipoproteins
^5HDL = High density lipoproteins
^6IMF = Intramuscular fat
^7ΣSFA = Sum of saturated fatty acids
^8ΣMUFA = Sum of monounsaturated fatty acids
^9ΣPUFA = Sum of polyunsaturated fatty acids
^10UI = Unsaturation index = 1 × (% monoenoics) +2 × (% dienoics) +3 ×
(% trienoics) +4 × (% tetraenoics) +5 × (% pentaenoics) +6 × (%
hexaenoics) [[122]126]
^11Σn3 = Sum of n-3 fatty acids
^12Σn6 = Sum of n-6 fatty acids
The assessment of differences in glucose and lipids metabolism
([123]Table 1) showed that purebred IB piglets have greater plasma
levels of total and HDL cholesterol, and triglyceride than IBxDU
neonates. These differences at birth are concordant with the similar
differences previously found between purebred Iberian and lean
crossbred (Large White x Landrace x Pietrain) fetuses [[124]47].
Cholesterol is of vital importance for the offspring as a key
constituent of cell membranes and the precursor of hormones and
metabolic regulators [[125]48, [126]49]. Placental and fetal tissues
have the capacity for de novo cholesterol synthesis [[127]50] but the
high demand from the fetuses makes the transport of maternal
cholesterol through the placenta necessary [[128]51, [129]52].
Triglycerides are also indispensable as a major source of energy for
the developing fetus and are also transferred from maternal circulation
[[130]53, [131]54]. Previous studies have found that adequate
availability of cholesterol and triglycerides is even more critical in
fatty-pigs breeds [[132]55], which have higher values of plasma lipids
indexes than lean breeds [[133]47].
These results reinforce that genetic differences between fatty-pigs and
lean breeds are established from prenatal stages and, together with
previous results, may also give evidence of a genetic predisposition
for lipid metabolism alterations in the Iberian breed. The same
findings regarding plasma cholesterol and triglycerides levels have
been reported in humans with familial combined hyperlipidemia, the most
common genetic form of hyperlipidemia in human [[134]56, [135]57] and
in the Rapacz familial hypercholesterolaemic swine model.
Regarding the IMF content and composition in BF ([136]Table 1), IB
showed almost 30% higher IMF content than IBxDU piglets (p = 0.014).
The genetic type affected IMF composition as well, IB pigs showing
greater ∑n-6/∑n-3 ratio (p = 0.031) (mainly due to greater proportion
of C18:2 n-6), and lower ∑SFA (saturated fatty acids) content (p =
0.035). Also, a trend for a higher oleic acid content was observed in
IB pigs (p = 0.092). As reported in previous studies, crossing with
Duroc sires decreased IMF concentration, in agreement with differences
observed in adult pigs [[137]9, [138]14].
The differences observed between IB and IBxDU in parameters such as
body size and weight, lipid metabolism-related indicators or IMF were
surprising taking in account the early stage of development. This
highlights the importance of improving the knowledge on molecular
aspects responsible for such phenotypic differences at very early ages
with the dual purpose of improving production in local breeds with
distinctive products and providing adequate models for human diseases.
Identification of differentially expressed genes by RNA-Seq analysis
An average of approximately 79 million sequence reads was obtained for
each individual sample; these were assembled and mapped to the
annotated Sscrofa10.2 genome assembly (22,861 genes). In all samples,
67–77% of the reads were categorized as mapped reads to the porcine
reference sequence.
The FPKM values were used to establish the total number of genes
expressed in muscle transcriptome (>0.5 FPKM). Approximately 50% of
total porcine annotated genes in the Sscrofa10.2 genome assembly were
expressed in the studied samples (an average of 11,392 genes out of
22,861 annotated genes).
Ninety-five genes were overexpressed in IB (FC ranging from 1.9 to 12)
and 54 genes were overexpressed in IBxDU (FC ranging from 2 to 27.2) (p
< 0.01)([139]S2 Table).
Large expression differences were observed for the genes MARCO (27.2x)
and CXCL13 (27.1x), which showed greater expression level in IBxDU than
in IB piglets. MARCO and CXCL13 genes are both related to immune
response. MARCO is also involved in cytoskeleton and cell morphology
determination of certain immune cells [[140]58] and CXCL13has also been
found to be upregulated in adipocytes when compared to preadipocytes
[[141]59].
On the other hand, another unidentified protein (ENSSSCG00000023287; FC
= 12.5x), ortholog of human MYL6B gene, and the pig genes CIART (8.4x)
and ATF3 (7.8x) were upregulated in IB piglets at birth. CIART is a
transcription repressor of the mammalian circadian clock that inhibits
the activators CLOCK and BMAL-1 [[142]60]. The mammalian circadian
clock regulates sleep–wake rhythms, body temperature, blood pressure,
hormone production, immune system or cell cycle (reviewed in
[[143]61]). It is also important for energy homeostasis regulation, as
multiple genes involved in nutrient metabolism and metabolically
related hormones such as insulin or leptin display rhythmic
oscillations [[144]62]. It has been reported that animals deficient in
BMAL-1 show altered lipid homeostasis (i.e. an increase in the levels
of circulating fatty acids, including triglycerides, free fatty acids,
and cholesterol) and metabolic syndrome [[145]63], which is in
agreement with the phenotypic results observed in IB pigs. ATF3 codes
for a transcription factor considered as an adaptation-response gene
involved in a variety of processes such as immunity, regulation of the
cell cycle and apoptosis [[146]64], and cellular stress response
[[147]65].
Some other interesting DE genes are related to muscle and adipose
tissue development, for example SLC2A4, FGF21, UBC, or ACHE
[[148]66–[149]69]. Moreover, three DE genes (DLK1, MYH10 and ZWILCH)
were also observed to be DE between IB and IBxDU in a previous study
[[150]14], where Longissimus dorsi transcriptome was compared at
weaning. DLK1 is a transmembrane protein expressed in preadipocytes but
not in mature adipocytes [[151]70], thus, greater expression in IB
neonates (2.6x) suggests greater number of undifferentiated
preadipocytes in IB than in IBxDU piglets, possibly associated with a
greater adipogenic potential. However, in a previous study, DLK1 gene
was found to be upregulated in IBxDU pigs at weaning [[152]14]. The
different age of sampling could explain the opposite results: IB
piglets may be born with higher amounts of preadipocytes that
differentiate faster than those from IBxDU pigs thus, leading to lower
preadipocyte content at weaning age. Similarly, a different pattern for
myocyte differentiation has been reported between high and low muscle
growth efficiency breeds, such as Landrace, Lantang, Pietrain or Duroc
[[153]29, [154]71]. In Landrace, myocyte differentiation develops
faster than in low efficiency breeds [[155]29, [156]71]. In addition,
the proliferation and differentiation of preadipocytes are stronger and
faster in Bamei than in Landrace (representing a fat- and lean-type pig
breed, respectively) [[157]17, [158]72]. Thus, we suggest a faster
adipocyte differentiation in IB pigs at an early age that may conclude
earlier than in IBxDU pigs.
MYH10 gene codes for a heavy-chain myosin, and was also found
upregulated in IBxDU pigs when compared to IB pigs at birth and at
weaning [[159]14], which suggests a greater development of muscular
cells in crossbreds. The gene ZWILCH is overexpressed in IB at both
ages; it is involved in cell proliferation and differentiation, and it
may also play a role in the control of adipogenesis [[160]73], thus
being an interesting candidate to explain phenotypic differences in
adipogenesis and lipogenesis.
In order to validate the results obtained from the RNA-Seq analysis,
the relative expression of some DE genes (upregulated in both genetic
types) and a few non-DE genes was assessed by qPCR in all the available
samples. Five out of the seven genes identified as DE by RNA-Seq
technology were validated by qPCR. The remaining two genes showed a
trend (p = 0.085) and a numerical difference (p = 0.105). Non-DE genes
showed similar results using both technologies ([161]S3 Table). A
concordance correlation coefficient, used to asses technical validation
in high throughput transcriptomic studies [[162]43] was calculated (CCC
= 0.93) and denoted a high general concordance between RNA-Seq and qPCR
expression values. In general good individual correlation values were
obtained ([163]S3 Table). Fold-Change and significance tended to be
greater when expression differences were analyzed by RNA-Seq
technology, in accordance with its higher sensitivity [[164]74].
Biological interpretation of the differential expression results
Different approaches were used to perform an exhaustive and robust
biological interpretation of the results obtained in the transcriptome
study. Results obtained from the Gene Ontology (GO) term
overrepresentation analysis, performed to detect active biological
processes differing in both IB and IBxDU, are shown in [165]Table 2. In
addition, IPA software was used to find biological functions
overrepresented in both genetic types ([166]S4 Table and Figs
[167]1–[168]3) and to identify pathways ([169]Table 3) and networks
([170]Fig 4) associated with the DE genes.
Table 2. Gene Ontology (GO) overrepresented terms regarding the biological
process category.
[171]^*GO term p-value q-value Term name
COMMON GO:0007165 9.54E-07 6.96E-05 Signal transduction
GO:0031323 3.81E-06 1.39E-04 Regulation of cellular metabolic process
GO:0009059 6.10E-05 6.37E-04 Macromolecule biosynthetic process
GO:0044267 1.22E-04 8.10E-04 Cellular protein metabolic process
GO:0010646 2.44E-04 1.11E-03 Regulation of cell communication
GO:0030154 2.44E-04 1.11E-03 Cell differentiation
GO:0048584 7.81E-03 1.84E-02 Positive regulation of response to
stimulus
IB[172] ^1 GO:0010467 1.22E-04 8.10E-04 Gene expression
GO:0043412 1.22E-04 8.10E-04 Macromolecule modification
GO:0036211 2.44E-04 1.11E-03 Protein modification process
GO:0009889 4.88E-04 1.98E-03 Regulation of biosynthetic process
GO:0090304 9.77E-04 3.56E-03 Nucleic acid metabolic process
GO:0019438 1.95E-03 5.48E-03 Aromatic compound biosynthetic process
GO:0016070 1.95E-03 5.48E-03 RNA metabolic process
GO:0023056 7.81E-03 1.84E-02 Positive regulation of signaling
IBxDU[173] ^2 GO:0002684 1.95E-03 1.86E-02 Positive regulation of
immune system process
GO:0050790 7.81E-03 3.49E-02 Regulation of catalytic activity
GO:0016337 7.81E-03 3.49E-02 Single organismal cell-cell adhesion
[174]Open in a new tab
^1 IB = Purebred Iberian pigs
^2 IBxDU = Iberian X Duroc crossbred pigs
*GO terms are considered either common, when they are enriched in both
genetic types and specific when the GO term is only enriched in one of
the two genetic types.
Fig 1. Enriched biological functions in Iberian (IB) pigs.
[175]Fig 1
[176]Open in a new tab
The network generated by IPA software shows enriched biological
functions in IB pigs (blue color) related to lipid metabolism and
muscle atrophy. Lipid synthesis is upregulated in IB pigs. Muscle and
protein degradation are also enriched, probably leading to a decreased
muscle growth in those pigs. The network includes several transcription
factors which may play key roles in the functional differences between
genetic types, in agreement with the regulators study (FOXO3, ATF3,
TRIM63).
Fig 3. Enriched biological functions potentially related to muscle growth in
crossbred (IBxDU) pigs.
[177]Fig 3
[178]Open in a new tab
The network generated by IPA software shows enriched biological
functions in IBxDU pigs (orange color) potentially related to muscle
growth and body size, functions regulated by a wide variety of genes,
from growth factors (FGF21) to immune system related genes (MARCO,
MSR1).
Table 3. Pathways significantly enriched in Purebred (IB) and Duroc-crossbred
(IBxDU) Iberian pigs.
IB[179] ^1 IBxDU[180] ^2
Pathway p-value Pathway p-value
PI3K Signaling in B Lymphocytes <0.001 Agranulocyte Adhesion and
Diapedesis <0.001
NRF2-mediated Oxidative Stress Response 0.0022 LXR/RXR Activation
<0.001
Glutamine Biosynthesis I 0.0031 Atherosclerosis Signaling <0.001
IL-17A Signaling in Fibroblasts 0.0051 Aldosterone Signaling in
Epithelial Cells 0.0012
April Mediated Signaling 0.0060 Granulocyte Adhesion and Diapedesis
0.0019
LXR/RXR Activation 0.0060 TREM1 Signaling 0.0050
Epoxysqualene Biosynthesis 0.0062 Protein Ubiquitination Pathway 0.0054
B Cell Activating Factor Signaling 0.0066 Citrulline-Nitric Oxide Cycle
0.0071
Thyronamine and Iodothyronamine Metabolism 0.0091 IL-12 Signaling and
Production in Macrophages 0.0155
Thyroid Hormone Metabolism I (via Deiodination) 0.0091 Superpathway of
Citrulline Metabolism 0.0195
Wnt/Ca+ pathway 0.0129 nNOS Signaling in Skeletal Muscle Cells 0.0209
Tight Junction Signaling 0.0145 Differential Regulation of Cytokine
Production in Macrophages and T Helper Cells by IL-17A and IL-17F
0.0251
AcutePhase Response Signaling 0.0151 Production of Nitric Oxide and
Reactive Oxygen Species in Macrophages 0.0263
ERK5 Signaling 0.0158 Clathrin-mediated Endocytosis Signaling 0.0275
Growth Hormone Signaling 0.0191 Actin Cytoskeleton Signaling 0.0372
Clathrin-mediated Endocytosis Signaling 0.0191 Role of p14/p19ARF in
Tumor Suppression 0.0417
Melatonin Signaling 0.0195 IL-17A Signaling in Fibroblasts 0.0468
Toll-like Receptor Signaling 0.0214
Regulation of IL-2 Expression in Activated and Anergic T Lymphocytes
0.0245
UVA-Induced MAPK Signaling 0.0295
Antioxidant Action of Vitamin C 0.0355
IGF-1 Signaling 0.0355
T Cell Receptor Signaling 0.0355
Role of IL-17A in Psoriasis 0.0389
Cholesterol Biosynthesis I 0.0389
Cholesterol Biosynthesis II (24. 25-dihydrolanosterol) 0.0389
Cholesterol Biosynthesis III (via Desmosterol) 0.0389
Vitamin-C Transport 0.0417
Glucocorticoid Receptor Signaling 0.0457
FattyAcid α-oxidation 0.0479
[181]Open in a new tab
^1 IB = Purebred Iberian pigs
^2 IBxDU = Iberian X Duroc crossbred pigs
Fig 4. Gene network containing DE genes related to Cellular Compromise,
Organismal Injury and Abnormalities and Skeletal and Muscular Disorders.
[182]Fig 4
[183]Open in a new tab
This gene network shows several genes functionally interconnected and
closely associated with the ubiquitin C (UBC), a central molecule in
the network responsible for protein degradation. These genes are
differentially expressed between genotypes and indicate the importance
of the ubiquitination process in muscle growth of newborn pigs.
Upregulated functions and pathways in Biceps femoris muscle from purebred
Iberian piglets
Enriched biological functions in IB piglets identified by IPA software
were mainly related with lipid and glucose metabolism (i.e.
Concentration of lipid, Synthesis of lipid, Insulin resistance,
Quantity of adipose tissue or Fatty acid metabolism; [184]Fig 1) and
with muscle growth ([185]S4 Table). In agreement, several DE genes with
known functions related to lipid metabolism such as FOS, FDFT1, SLC4A2,
TRIM63, EPXH1, ALOX12B, FOXO3A or ACHE were overexpressed in IB,
possibly associated to the higher amount of adipose tissue observed in
BF muscle of IB when compared to IBxDU pigs ([186]Table 1). Similar
results have been found in IB and IBxDU piglets at 28 days of age in
loin muscle [[187]14]. However, in the previous study, a different set
of DE genes related to lipid metabolism was found, probably due to the
high complexity of processes and molecular mechanisms regulating lipid
metabolism at that stage of development. In another study comparing the
muscle transcriptome of two divergent breeds for muscle and fat
deposition, several muscle metabolism-related genes were identified as
potential regulators of IMF deposition [[188]29], such as FOS and ABRA
genes, identified also in the present study.
Accordingly, several pathways enriched in IB pigs (p < 0.05)
([189]Table 3) were related to the control of energy homeostasis
(Wnt/Ca+, Glutamine Biosynthesis or Fatty Acid α-oxidation pathways)and
to protein synthesis and cell growth (Growth Hormone (GH) Signaling and
IGF-1 Signaling); IGF-1 is essential during prenatal and GH during
postnatal growth [[190]75]. The effect of GH and IGF-1 on adipose
tissue development and metabolism is controversial, as both have been
proposed to play either a positive or a negative role on adipocyte
differentiation [[191]76, [192]77]. Thus, in IB newborns, the
activation of these pathways might be associated to both muscle and
preadipocytes development and differentiation. Accordingly, the
adipogenesis pathway showed a trend for enrichment (p = 0.055).
Upregulated functions and pathways in Biceps femoris muscle from
Duroc-crossbred Iberian piglets
Enriched functions in IBxDU newborns ([193]S4 Table) were related to
cell growth and differentiation, such as Recruitment and transmigration
of cells, Hyperplasia or Transformation of fibroblast cell lines
([194]Fig 2), and to protein and muscle organization and accretion such
as Mass of organism, Quantity of protein in blood and Organization of
cytoskeleton ([195]Fig 3). Consistently, DE genes in IBxDU pigs
clustered in 3 gene networks related to Cell-To-Cell Signaling and
Interaction, Hematological System Development and Function, Immune Cell
Trafficking, Post-Translational Modification, Protein Folding, Cellular
Assembly and Organization, Cellular Movement, Cell Signaling and
Cellular Function and Maintenance.
Fig 2. Enriched biological functions related to cell growth in crossbred
(IBxDU) pigs.
[196]Fig 2
[197]Open in a new tab
The network generated by IPA software shows enriched biological
functions in IBxDU pigs (orange color) involved in cellular growth and
differentiation. Upregulation of several chemokines as well as a myosin
in IBxDU pigs, and downregulation of some transcription factors (FOXO3,
ATF3) are responsible for the enrichment of such functions in
crossbreds.
Some IBxDU upregulated genes involved in those functions were
chemokines (i.e.CCL19, CCL2 or CXCL13) and chemokine receptors (ACKR1),
which have been reported to be involved in the process of myogenesis
and muscle regeneration [[198]78]. In our study, chemokines were
associated with functions such as transmigration and recruitment of
cells and organization of cytoskeleton (Figs [199]2 and [200]3), which
could be related to muscle cell growth. Moreover, MYH10, a non-muscle
myosin that regulates actin cytoskeleton remodeling, was found
overexpressed in IBxDU and involved in cytoskeleton reorganization, a
critical process during myogenesis [[201]79]. Consistently, the Actin
Cytoskeleton Signaling pathway was also enriched in IBxDU pigs in both
the present study and at 28 day of age [[202]14].
Mass of organism, associated with the overexpression of PRSS8, was
predicted by IPA software to be enriched in IBxDU piglets, in
accordance with the phenotypic differences in body weight and size
found between pure and crossbred pigs. In agreement, a decrease in body
weight was reported in mutant rats for gene PRSS8 [[203]80].
The enrichment of protein and muscle development-related pathways (p <
0.05), such as Citrulline-Nitric Oxide Cycle, Superpathway of
Citrulline Metabolism or nNOS Signaling in Skeletal Muscle, support the
greater muscle development in IBxDU piglets. Citrulline is a
non-essential amino acid that plays a role in muscle development and
affects body composition in rats, increasing lean and decreasing
adipose tissue when it is provided in the diet [[204]81]. The nNOS
Signaling in Skeletal Muscle Cells modifies the blood flow to the
muscle [[205]82]. Thus, the activation of these pathways suggests a
greater nutrient input to muscle tissues in those pigs.
Upregulated functions and pathways in Biceps femoris muscle from both
Duroc-crossbred and purebred Iberian piglets
Development-related processes depend on the genetic background, but
also on the growth stage. At birth, when growth is a very active
process, several functions and pathways enriched in both genetic types
were observed. However, these common processes ended up in different
phenotypic consequences. This might be because developmental mechanisms
such as muscle growth or immune system development are tightly
regulated and the final output depends on both the balance of
activating and inhibiting pathways and the expression levels of genes
involved in such processes.
Regarding cellular growth, the GO terms enrichment analysis retrieved
several GO terms related to normal cell cycle and growth, such as cell
differentiation, cellular protein metabolic process or macromolecule
biosynthetic process that were common between both genetic types,
supporting the importance of these processes in the early development
of pigs. Accordingly, in 3 months-old Casertana pig (an autochthonous
Italian fatty pig breed), GO terms close to the aforementioned were
upregulated when compared to 6 months old pigs [[206]83].
The importance of protein metabolism at birth is supported by the
activation of the protein ubiquitination process in both genetic types.
The Protein Ubiquitination Pathway was enriched in crossbred piglets (p
= 0.005; [207]Table 3), where mainly members of the family of Heat
Shock Proteins (HSPH1, HSPA4L and DNAJA1) were found upregulated in
IBxDU pigs. DNAJA1, together with HSP27 have been reported to play a
role in both cellular stress and meat quality; their expression was
found to decrease IMF content [[208]84], consistently with the results
observed in IBxDU pigs. However, no enrichment in atrophy or
degradation of muscle-related functions was observed in IBxDU piglets.
Conversely, IB muscle showed activation of functions related to muscle
atrophy ([209]Fig 1), mainly driven by genes that stimulate muscle
atrophy by means of the ubiquitin proteasome system (UPS) [[210]85].
This system was also upregulated in IB when compared to IBxDU pigs at
28 days of age [[211]14]; and in Basque (a French pig breed with low
lean meat and high fat contents) when compared to Large White pigs
[[212]86]. Supporting the evidence of an activated protein
ubiquitination process, the gene coding for Ubiquitin C (UBC) and
FBXO32 (one of the three ligase enzymes, together with TRIM63 of the
UPS) were also overexpressed in IB pigs.
In accordance with these results, genes implicated in these processes
were found in the most significant gene network involving DE genes in
both genetic types. This network was associated with Skeletal and
Muscular Disorders, as well as with Cellular Compromise, Organismal
Injury and Abnormalities ([213]Fig 4) and showed several genes related
to protein catabolism (UBC, HSP, TRIM63 or FOXO3) among the most
central genes in the network, which supports the importance of these
genes in protein catabolism, a very active process in developing IBxDU
and especially in IB pigs.
Iberian pigs show low protein accretion potential, although higher
relative protein synthesis rate has been reported in IB than in lean
Landrace pigs [[214]87]. Since IB muscles are smaller than Landrace
ones, it was suggested that a higher protein turnover rate in IB pigs
should be responsible for the differences between protein synthesis and
deposition [[215]87]. In this context, it seems evident that this
increased protein turnover in IB pigs is related to the activated
muscle atrophy-related genes found in this study. Moreover, Damon et
al. (2012) suggested that the activation of the UPS might affect meat
tenderness, since proteasome would be one of the main endogenous
proteolityc systems contributing post-mortem meat tenderization.
The UPS is also interconnected to the unfolded protein response (UPR),
activated under endoplasmic reticulum stress situations. The UPR is
involved in degradation of un and misfolded proteins and recently
associated to the control of adipogenesis under situations of
endoplasmic reticulum stress [[216]88, [217]89]]. It is noteworthy that
ATF3, a transcription factor inducible by endoplasmic reticulum stress
[[218]65] was upregulated in IB pigs. This could suggest certain level
of endoplasmic reticulum stress in IB pigs that might be related to
preadipocyte differentiation.
The importance of cholesterol in early stages of development could be
the cause for the activation of pathways related to cholesterol
metabolism and biosynthesis in both genetic types. For example, the
LXR/RXR activation pathway, involved in the regulation of lipid
metabolism, inflammation and cholesterol to bile acid catabolism, was
enriched in both IB and IBxDU newborns (p = 0.006 and p < 0.001,
respectively). On the contrary, at 28 days of age, this pathway was
enriched in IB pigs [[219]14]. The atherosclerosis signaling pathway
(enriched in IBxDU pigs) or the Cholesterol Biosynthesis and the
Epoxysqualene Biosynthesis pathways (enriched in IB pigs) were also
identified as relevant pathways. However, phenotypic observations
reflect a more active cholesterol and triglycerides biosynthesis in IB
than in IBxDU piglets ([220]Table 1).This could be due to the
upregulation of different genes in each genetic type: upregulated genes
in IB pigs associated with these pathways were closely related to the
cholesterol and lipid metabolism (FDFT, TF and APOM). Accordingly,
genes related to cholesterol and lipid metabolism such as RXRG, USF1,
LPL or some apolipoproteins have been proposed as candidate genes
involved in the familial combined hyperlipidemia in human [[221]56,
[222]90, [223]91], characterized by high levels of plasma cholesterol
and triglycerides, as observed in IB pigs. Moreover, the aforementioned
upregulated genes in IB pigs were connected in a network associated
with functions such as Amino Acid Metabolism, Molecular Transport and
Small Molecule Biochemistry ([224]Fig 5). On the other hand, genes
involved in such cholesterol-related pathways that were upregulated in
IBxDU piglets where associated to the immune response (MSR1 and CCL2).
Fig 5. Gene network containing genes upregulated in Iberian (IB) pigs related
to Amino Acid Metabolism, Molecular Transport and Small Molecule
Biochemistry.
[225]Fig 5
[226]Open in a new tab
The gene network shows interactions between cholesterol and lipid
metabolism-related genes (LDL, APOM, FDFT1, SLC2A4, etc.) upregulated
in Iberian pigs.
The immune system is not fully developed at birth in pigs [[227]92],
and thus, related pathways and functions were upregulated in both in IB
and IBxDU newborns. Some of these pathways include PI3K Signaling in B
Lymphocytes, April Mediated Signaling or B Cell Activating Factor
Signaling (enriched in IB pigs) and Agranulocyte Adhesion and
Diapedesis or IL-12 Signaling and Production in Macrophages (enriched
in IBxDU pigs). Moreover, the pathway IL-17A Signaling in Fibroblasts
was enriched in both genetic types. On the other hand, the GO term
Positive regulation of immune system process was enriched in IBxDU but
not in IB piglets. This is consistent with the aforementioned
upregulation of genes such as MARCO and CXCL13 in Duroc-crossbred pigs,
suggesting a more developed immune system in IBxDU than IB newborns. It
has been previously reported that the immune system is affected by pig
breed [[228]93] and domestication [[229]94]. An overrepresentation of
genes related to the immune system was reported in wild boars when
compared to domestic pig breeds [[230]95]. However, to our knowledge
there are no previous studies assessing differences in immune
efficiency introduced by sire line.
Immune system-related genes are also involved in numerous other
biological processes, such as fat accumulation, closely associated with
inflammation [[231]96, [232]97]. Thus, the enrichment of pathways
related to the immune system could be related to lipid accumulation,
although it has also been related to lower backfat thickness in pigs
[[233]98]. It is noteworthy that no biological functions related to the
immune system were found enriched in IB pigs, which might suggest that
upregulation of immune system-related pathways might be related to
adipose tissue development rather than to the immune system development
in IB pigs.
The reported active biological functions in IB and IBxDU pigs suggest a
different predisposition for cell and tissue growth between genetic
types, being IB pig metabolism intended to energy storage and IBxDU pig
metabolism to cell growth and differentiation. These differences in
biological regulation are in agreement with phenotypic differences
found between the two genetic types.
Regulatory transcription factors
We performed a regulatory factors study to investigate the driving
molecular mechanisms responsible for the differences in gene expression
observed between genetic types.
A total of 723 TRF, obtained from the animal TFDB, showed expression
values above 0.5 FPKM in our pigs and were thus considered as
expressing TRF. Among them, 92 TRF potentially affected the gene
expression profile in IB and IBxDU muscles ([234]Table 4). We
considered TRF that were either DE (7 TRF), identified by IPA software
as regulators (45 TRF) or identified in the RIFs study (48 TRF)
([235]Table 4).
Table 4. Potential regulators affecting gene expression that are: a)
differentially expressed (DE) between IB and IBxDU; b) identified by
Ingenuity Pathways Analysis (IPA) software or c) identified by RIFs study.
Ensembl ID GENE DE regulators IPA—Regulators RIFs study
p-value FC Change Z-score p-value RIF1(z) RIF2(z)
ENSSSCG00000006036 ABRA 2.27E-06 -4.37
ENSSSCG00000008123 ARID5A 2.01
ENSSSCG00000015972 ATF2 4.40E-02
ENSSSCG00000015595 ATF3 2.48E-04 -7.76
ENSSSCG00000000084 ATF4 1.42E-02
ENSSSCG00000012241 BCOR 4.17E-03
ENSSSCG00000013397 BMAL1 -2.22
ENSSSCG00000008377 CCT4 -2.11
ENSSSCG00000002866 CEBPA -2.426 7.32E-03
ENSSSCG00000002867 CEBPG 3.60
ENSSSCG00000006752 CSDE1 2.45
ENSSSCG00000011274 CTNNB1 -2.183 1.01E-02
ENSSSCG00000014336 EGR1 -0.254 5.28E-03
ENSSSCG00000010224 EGR2 2.01E-02 2.21
ENSSSCG00000008443 EPAS1 4.01E-02
ENSSSCG00000002383 FOS 1.38E-03 -4.39 0.381 6.14E-03
ENSSSCG00000012967 FOSL1 2.76E-02
ENSSSCG00000007576 FOXK1 2.35
ENSSSCG00000009370 FOXO1 1.03E-02 3.22
ENSSSCG00000004387 FOXO3 7.04E-04 -2.72 -2.364 5.05E-06
ENSSSCG00000001619 FOXP4 2.79
ENSSSCG00000015733 GLI2 2.60
ENSSSCG00000007720 GTF2IRD1 1.66E-02
ENSSSCG00000000846 HCFC2 2.34
ENSSSCG00000014388 HDAC3 2.36E-02
ENSSSCG00000010472 HHEX 4.89E-02
ENSSSCG00000004138 HIVEP2 1.70E-04
ENSSSCG00000009327 HMGB1 3.67E-03
ENSSSCG00000009704 HMGB2 2.55
ENSSSCG00000008898 HOPX 1.66E-02
ENSSSCG00000015985 HOXD3 3.59
ENSSSCG00000005917 HSF1 9.28E-05
ENSSSCG00000004238 HSF2 6.06E-03
ENSSSCG00000014277 IRF1 2.49E-02 2.19
ENSSSCG00000003178 IRF3 3.13E-02
ENSSSCG00000012853 IRF7 3.08E-02
ENSSSCG00000008119 KCNIP3 2.34
ENSSSCG00000010928 KDM5B -2.10
ENSSSCG00000006928 LMO4 2.64
ENSSSCG00000004528 MBD2 2.44E-02
ENSSSCG00000000552 MED21 2.35
ENSSSCG00000005720 MED27 2.18
ENSSSCG00000006482 MEF2D 6.93E-03
ENSSSCG00000012534 MORF4L2 -2.11
ENSSSCG00000013507 MPND 2.18
ENSSSCG00000012114 MSL3 2.34
ENSSSCG00000017882 MYBBP1A 2.07E-02
ENSSSCG00000013375 MYOD1 7.26E-03
ENSSSCG00000015475 MYOG 8.17E-06
ENSSSCG00000010399 NCOA4 3.59
ENSSSCG00000015987 NFE2L2 -2.206 2.44E-02
Ensembl ID GEN DE regulators IPA—Regulators RIFs study
p-value FC Z-score p-value RIF1(z) RIF2(z)
ENSSSCG00000001952 NFKBIA -1.279 3.40E-02
ENSSSCG00000001703 NFKBIE 4.09E-02
ENSSSCG00000005385 NOR-1 8.17E-04 -6.56
ENSSSCG00000009856 NOS1 9.85E-03 2.29 3.78E-05
ENSSSCG00000006689 PIAS3 3.69E-02
ENSSSCG00000014437 PPARGC1B 1.65E-02
ENSSSCG00000009746 RAN -2.25
ENSSSCG00000009401 RB1 1.172 2.29E-02
ENSSSCG00000008388 REL 2.04E-03
ENSSSCG00000012981 RELA 0.749 6.16E-04
ENSSSCG00000017071 SAP30L -2.28
ENSSSCG00000011201 SATB1 -2.32
ENSSSCG00000001880 SIN3A -2.12
ENSSSCG00000004952 SMAD2 2.43
ENSSSCG00000004524 SMAD4 4.94E-02
ENSSSCG00000005232 SMARCA2 -2.47
ENSSSCG00000013629 SMARCA4 4.00E-02
ENSSSCG00000006256 SOX17 3.29E-02
ENSSSCG00000017403 STAT3 3.52E-02
ENSSSCG00000017406 STAT5B 2.22
ENSSSCG00000010638 TCF7L2 2.68E-02
ENSSSCG00000001092 TDP2 2.44
ENSSSCG00000001544 TEAD3 -2.06
ENSSSCG00000017950 TP53 -1.444 1.21E-04 3.48 1.69
ENSSSCG00000007208 TRIB3 2.84
ENSSSCG00000027684 TRIM63 7.02E-04 -2.46 1.25E-02
ENSSSCG00000015370 TWIST1 3.28E-03
ENSSSCG00000006790 WDR77 -2.18
ENSSSCG00000011953 ZBTB11 2.10
ENSSSCG00000007958 ZNF174 2.31
ENSSSCG00000003233 ZNF175 2.86
ENSSSCG00000002838 ZNF423 2.26
ENSSSCG00000003070 ZNF45 2.24
ENSSSCG00000007767 ZNF668 2.26
ENSSSCG00000001206 ZSCAN26 -2.84
ENSSSCG00000002877 ZNF181 2.50
ENSSSCG00000003244 2.05
ENSSSCG00000011620 -2.41
ENSSSCG00000011943 -2.07
ENSSSCG00000013384 1.96
ENSSSCG00000016700 3.43
[236]Open in a new tab
FC: Fold-Change
Z-score: reflects the activation state of predicted transcriptional
regulators. It is based on the experimentally observed gene expression,
and on literature‐derived regulation direction information, which can
be either “activating” or “inhibiting”
RIF1 (z) extreme scores identify those transcription factors that are
consistently most differentially co-expressed with highly abundant and
highly DE genes. Bootstrap 95% and 99% confidence intervals for RIF1
z-scores: -1.996/2.074 and -2.883/2.918, respectively.
RIF2 (z) extreme scores identify transcription factors with the most
altered ability to predict the abundance of DE genes. Bootstrap 95% and
99% confidence intervals for RIF2 z-scores: -2.036/1.953 and
-2.609/2.490, respectively
In the present study, we found 7 TRF showing differences in expression
level between genetic types. Most of them (ABRA, ATF3, FOS, FOXO3,
NOR1, TRIM63) were upregulated in IB. These genes are related to muscle
and adipose tissue development ([237]Fig 1) and most of them have been
mentioned in the previous section, but some of them may be highlighted
and deserve additional comments. For example, one of the genes with
greater expression differences (ATF3) was also an identified regulator
in the Animal TFDB, suggesting an important role in the gene expression
differences observed. NOR1 showed also an important expression
difference (6.56x) and codes for a nuclear receptor involved in a wide
array of functions such as inflammation, cell cycle regulation,
apoptosis, steroidogenesis, adipogenesis, angiogenesis and energy
metabolism. Moreover, it has been found overexpressed in obese when
compared to normal humans and its expression went back to normal values
after fat loss [[238]99]. These findings are in accordance with the
results of the present study, where animals with higher IMF content
showed greater muscular expression of NOR1.
IPA software is a potent tool to identify regulators based on previous
knowledge and bibliographic references. On the other hand, the RIFs