Transcriptome analysis reveals differential gene expression in intramuscular adipose tissues of Jinhua and Landrace pigs.

Meat is a rich source of protein, fatty acids and carbohydrates for human needs. In addition to necessary nutrients, high fat contents in pork increase the tenderness and juiciness of the meat, featuring diverse application in various dishes. This study investigated the transcriptomic profiles of intramuscular adipose tissues in Jinhua and Landrace pigs by employing advanced RNA sequencing. Results showed significant interesting to note that there were significant differences in the expression of genes. 1,632 genes showed significant differential expression, 837 genes were up-regulated and 195 genes were down-regulated. Variations in genes responsible for cell aggregation, extracellular matrix formation, cellular lipid catabolic process, and fatty acid binding strongly supported that both pig breeds feature variable fat and muscle metabolism. Certain differentially expressed genes are included in the pathway of mitogen-activated protein kinase signaling pathway, Ras signaling pathway and insulin pathway. Results from real-time quantitative polymerase chain reaction also validated the differential expression of 17 mRNAs between meats of the two pig breeds. Overall, these findings reveal significant differences in fat and protein metabolism of intramuscular adipose tissues of two pig breeds at the transcriptomic level and suggest diversification at the genetic level between breeds of the same species.

Pigs are commonly used for meat production in China, whereas their derivatives are required by millions of people [22]. Meat quality is mainly determined by marbling or intramuscular fat contents [21]. Fat contents in intramuscular adipose tissues enhance juiciness, flavor, and palatability of pig meat [4, 5]. Jinhua, a famous typical pig breed in China, is now widely used for ham production, and exhibits high intramuscular fat contents, high prolificacy, and low growth [10]. By contrast, Landrace is an improved lean pig breed characterized by high growth and body weight [13, 16]. Thus, the two porcine breeds could be considered as comparison models to investigate genetic differences and molecular mechanisms underlying the above-mentioned phenotypic differences.

Adipose is a metabolic tissue that plays an important role in fatty acid synthesis, it serves as a mediator for differentiation and metabolism of adipocytes [23]. Downstream adipose products, adipokines, actively participate in maintenance of metabolic homeostasis [1, 12]. Although the regulatory mechanisms of adipogenesis in mammals have been studied [3, 8, 15], the differences between Jinhua and Landrace pigs remain unreported.

Identification of genes along with their molecular mechanisms has been widely conducted recently through microarray, real-time quantitative PCR (qPCR), and RNA sequence (RNA-seq) technologies [11, 27]. RNA-seq technology is a transcriptome analysis tool that provides excellent advantages to explore novel genes and their related functional properties [30]. Transcriptome comparison of adipose and muscle tissues from different porcine breeds is a new approach for studying functional genes during adipogenesis [15]. RNA-seq also enables the detection of differentially expressed genes (DEGs) between different ages and breeds of mammals [19, 20].

In this study, two porcine breeds, Jinhua and Landrace, were selected as models for transcriptome profile generation by RNA-seq technology. A total of 1,632 genes showed significant differential expression between the two breeds, and those genes were then subjected to Swiss-prot, Gene Ontology (GO), Cluster of Orthologous Groups of Proteins (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to explore gene regulation changes that can influence adipose tissue generation.

MATERIALS AND METHODS

Sample preparation

The protocol for sample isolation was performed in accordance with the Animal Care and Use Statute of China, and was approved by the animal care committee of Henan Institute of Science and Technology. Two varieties of intramuscular adipose tissue (three samples for each breed) were used. These varieties were derived from Jinhua and Landrace pigs. Jinhua and Landrace pigs (three pigs from each breed, 180 days old) were reared under the same conditions, fed corn and soybean meal, and formulated with trace mineral and vitamins to meet the requirements of National Research Council for different ages. After slaughtering, pieces of longissimus dorsi muscle were immediately collected. Then, intramuscular fat tissue was separated from muscle tissue within 10–20 min [24]. All collected samples were quickly frozen in liquid nitrogen within 30 min and then stored at −80°C for total RNA extraction.

RNA quantification and qualification

Total RNA from frozen samples was extracted using RNAprep Pure Tissue Kit following the manufacturer’s instructions (Tiangen Biotech Co., Ltd., Beijing, China). RNA degradation and contamination were monitored on 1% agarose gels. Purity, concentration and integrity of extracted RNA were determined using NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, U.S.A.) at OD260/280 ≥1.8, Qubit® RNA Assay Kit in Qubit®2.0 Flurometer (Life Technologies, Carlsbad, CA, U.S.A.). Agilent analysis of isolated RNA was assessed by RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, U.S.A.), and RNA with number ≥8.0 was considered for downstream application.

Library construction and sequencing

A total of 1 µg RNA per sample was used as input material for RNA sample preparations. Library sequences were generated using NEBNext®UltraTM RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, U.S.A.) following manufacturer’s recommendations, and index codes were added in attribute sequences to each sample. Briefly, mRNA was purified using poly-T oligo-attached magnetic beads, and fragmented to small pieces by using divalent cations under elevated temperature. RNA was reverse transcribed into cDNA. After end repair, adapter ligation and fragment selection, PCR was performed on these transcripts. The transcripts were then sequenced for library quality assessment on the Agilent Bioanalyzer 2100 system.

Clustering of index-coded samples was performed on a cBot Cluster Generation System by using TruSeq PE Cluster Kit v4-cBot-HS (Illumina) in accordance with the manufacturer’s instructions. After cluster generation, library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated (Biomarker Technologies Co., Beijing, China). To obtain clean and high-quality reads, clean reads were mapped to the reference genome sequence (ftp://ftp.ensembl.org/pub/release-75/fasta/sus_scrofa/) by using Tophat2 tool software for alignment. Only mapped reads were used for analysis. Cufflinks software was used for mapped read splicing, and reads were then compared with genomic annotation information to explore new genes.

DEG analysis and gene functional annotation

Fragments per kilobase of transcript per million fragments mapped (FPKM) was used to quantify the expression level of unigenes. The equation is as follows:

Differential expression analysis was performed using EBSeq software [7]. GO enrichment analysis of DEGs was implemented by the GOseqR packages based Wallenius’ non-central hyper-geometric distribution, whereas KOBAS (Mao et al., 2005) software was used to test the statistical enrichment of DEGs in KEGG pathways [9, 29].

qPCR analysis

To validate the results of sequencing analysis, 17 mRNAs were selected for qPCR analysis. Total RNAs from each sample were extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, U.S.A.), and All-in-oneTM First-Strand cDNA Synthesis Kit was used for reverse transcription PCR in accordance with the manufacturer’s instructions (GeneCopocie, Inc., Rockville, MD, U.S.A.). qPCRs were performed using 2 × All-in-oneTMqPCR mix with 2 µl of cDNA template and a final concentration of 0.2 µl All-in-oneTMqPCR primer. The reactions were carried out as follows: 10 min at 95°C (1 cycle), 10 sec at 95°C, 20 sec at 60°C and 10 sec at 72°C (40 cycles). mRNAs primers were used in qPCR (Supplementary file, Table 1).

Statistical analysis

False Discovery Rate (FDR) <0.05 and log2 Fold change (FC) ≥2 were set as thresholds for DEGs screening. Landrace breed was taken as the control group. 2-ΔΔCt (relative value of Jinhua and Landrace breeds) was used to calculate the mRNA expression levels between two breeds.

RESULTS

Illumina sequencing

In this study, RNA-seq technology was used to investigate the transcriptome intramuscular adipose tissues isolated from two porcine breeds. First, six cDNA libraries were constructed, representing the samples from Jinhua (L07, L08 and L09) and Landrace pigs (L10, L11 and L12). The six libraries were sequenced, resulting in 94.92 Gb clean data. More than 91.82% of the data yielded a high-quality score (Q30) (Table 1). GC content from these clean data reached more than 53.22%.

Table 1.

Sample sequencing data in L07, L08, L09, L10, L11 and L12

Name	BaseSum	Mapped reads
L07	15,459,125,786	79,502,171
L08	16,618,345,512	83,768,791
L09	14,058,501,046	70,956,817
L10	18,617,119,084	94,045,412
L11	16,344,667,718	86,992,608
L12	13,824,949,520	72,051,261

BaseSum: Total number of bases for clean data. Mapped reads: Number of reads aligned to reference genome.

Then, all reads were mapped onto the genome of sus_scrofa by Tophat2 software. A total of 79,502,171; 83,768,791; 70,956,817; 94,045,412; 86,992,608; and 72,051,261 reads were mapped, and they were used for further analysis (Table 1). After gene assemble and comparison, 1,868 new genes were observed. Through screening with databases, including Swiss-prot, GO, COG, and KEGG, 1,063 (56.91%) novel genes can be annotated (Table 2). The majority of mRNAs (17.86%) comprised long reads with more than 3,000 bp, 0.29% of mRNAs with 200 bp, and none measured less than 200 bp (Fig. 1).

Table 2.

Number of new gene annotation

Annotated database	COG	GO	KEGG	Swissprot	All
Annotated number	114	465	388	470	1,063

Fig. 1.

Size distribution of mRNAs. The histogram shows the number of mRNAs.

Differential expression analysis

Expression of mRNAs in six libraries was evaluated by FPKM value. The expression level of genes ranged from 0 FPKM to 7264.98 FPKM. A total of 1,632 genes showed significant differential expression between two samples; 837 genes were up-regulated, and 795 genes were down-regulated (Supplementary file, Table S2). Figure 2 shows the volcano plot for the differential expression level of genes between two samples.

Fig. 2.

Volcano plot of DEGs from samples of Jinhua and Landrace pigs. X-axis values correspond to log2 (fold change), and y-axis values represent –log10 (false discovery rate). Green points represent significantly down-regulated genes, blue points represent no difference in genes, and red points represent significantly up-regulated genes.

Functional enrichment analysis

All DEGs were annotated via Blast against GO, COG and KEGG database. A total of 1,440, 502 and 959 DEGs were annotated, respectively. Figure 3 shows the GEGs annotated to the GO database. Among the 61 functional classes, chemorepellent activity, cell aggregation, and extracellular matrix part presented the largest differences between DEG unigene and all unigenes (Fig. 3). Highly down-regulated DEGs were annotated by GO database to processes including negative regulation of growth, cellular lipid catabolic process, fatty acid binding, and lipid binding. DEGs were then annotated to the COG database, which is constructed based on the evolutionary relationship of bacteria, algae, and eukaryotes. Twenty-five classes, except nuclear structure and extracellular structures, were included, and the top three enriched genes were represented in general functional prediction only, signal transduction mechanisms and replication, recombination, and repair (Fig. 4). Among the 25 functional groups, lipid transport and metabolism accounted for 5.15% of total annotated DEGs. KEGG enrichment analysis showed that cancer pathway was the highest enrichment pathway, and 32, 21 and 20 DEGs are involved in mitogen-activated protein kinase (MAPK), Ras, and insulin pathways, respectively (Fig. 5).

Fig. 3.

GO biological processes enrichment of target genes for DEGs. DEG Unigene: unigene of differentially expressed mRNAs in each secondary function; All Unigene: unigenes of all mRNAs in each secondary function.

Fig. 4.

COG category enrichment of DEGs. X-axis represents 25 function classes shown on the right side of the figure, and y-axis represents the frequency of each class in the total class.

Fig. 5.

KEGG pathways category enrichment in target genes of DEGs. The bar represents the number of unigene ratios for relevant pathway.

Analysis of expression patterns by qPCR

qPCR results showed that the expression patterns of 17 genes from Jinhua pigs varied and were higher than those of Landrace pigs (Fig. 6). CPT1A gene was the most statistically significant expressed gene, which was 762.12 times more than that of Landrace sample group (Supplementary file, Table S3). SLC27A6 gene was the lowest DEG.

Fig. 6.

qPCR-detected expression levels of 17 selected mRNA.

DISCUSSION

RNA-seq is an effective method to identify novel genes and their potential functional properties. Recent studies have confirmed the contribution of transcriptome analysis at the molecular level [18, 19, 25]. The current results indicated the accuracy and high quality of sequencing data, which were then used for subsequent research. Differences in gene expression profiles can unveil mechanisms underlying biological activities of various genes. A total of 1,632 DEGs were observed, and the numbers of genes that were up- or down-regulated were similar (837 versus 795). Several DEGs have been reported to play a crucial role in regulating lipid homeostasis; for example, peroxisome proliferator-activated receptors can promote fatty acid transport and lipid droplet formation via binding to the promoter of key genes [17]. In the present study, the expressions of lipid synthesis-associated genes and acetyl-CoA carboxylase alpha and fatty acid synthase were also markedly and differentially expressed.

GO database comprises three main groups, including biological process, molecular function, and cellular component. We observed that parts of negative regulation of growth, cellular lipid catabolic process, fatty acid binding, and lipid binding feature highly down-regulated DEGs, which were annotated by GO database. Our results suggest considerable changes in gene expression in fat metabolism between the two porcine breeds. Similar phenomena have been reported for fat and muscle development [20, 28]. Annotation to COG database showed high enrichment in the groups of lipid transport and metabolism. This finding suggests the potential involvement of DEGs in lipid metabolism. However, this deduction requires further investigation. KEGG enrichment analysis will aid in exploring the relation between DEGs and pathways. Enriched MAPK signaling, Ras signaling, and insulin pathways accounted for a certain part of DEGs. MAPK signaling, Ras signaling, and insulin pathways are accompanied within lipid metabolism and adipocyte proliferation [6, 26]. These results suggest that the functional activities of these DEGs may be correlated with the aforementioned pathways.

To validate the expression data of DEGs, 17 selected genes involved in fat metabolism were analyzed from Jinhua and Landrace pigs. CPT1A was the most DEG. CPT1A gene can catalyze the transfer of the acyl group of long-chain fatty acid-CoA conjugates to carnitine and plays an important role in triglyceride metabolism [2, 14]. Up-regulation of CPT1A gene in intramuscular adipose from Jinhua pigs may play an important role in intramuscular adipose metabolism through the PPAR signaling pathway.

The present study validated the application of RNA-seq as a reliable tool to determine gene expression profiles in various tissues. By using RNA-seq, we showed the DEGs responsible for fat and protein metabolism in two porcine breeds. However, further studies must be performed to determine the differences between genes responsible for other cell signaling and metabolic pathways.