Integration of Transcriptome and Whole Genomic Resequencing Data to Identify Key Genes Affecting Swine Fat Deposition

Abstract

Fat deposition is highly correlated with the growth, meat quality, reproductive performance and immunity of pigs. Fatty acid synthesis takes place mainly in the adipose tissue of pigs; therefore, in this study, a high-throughput massively parallel sequencing approach was used to generate adipose tissue transcriptomes from two groups of Songliao black pigs that had opposite backfat thickness phenotypes. The total number of paired-end reads produced for each sample was in the range of 39.29–49.36 millions. Approximately 188 genes were differentially expressed in adipose tissue and were enriched for metabolic processes, such as fatty acid biosynthesis, lipid synthesis, metabolism of fatty acids, etinol, caffeine and arachidonic acid and immunity. Additionally, many genetic variations were detected between the two groups through pooled whole-genome resequencing. Integration of transcriptome and whole-genome resequencing data revealed important genomic variations among the differentially expressed genes for fat deposition, for example, the lipogenic genes. Further studies are required to investigate the roles of candidate genes in fat deposition to improve pig breeding programs.

About the Authors:

Kai Xing

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Feng Zhu

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Liwei Zhai

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Huijie Liu

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Yuan Wang

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Zhijun Wang

Affiliation: Tianjin Ninghe primary pig breeding farm, Ninghe, 301500, Tianjin, China

Shaokang Chen

Affiliation: Animal husbandry and veterinary station of Beijing, Beijing, 100125, Beijing, China

Zhuocheng Hou

* E-mail: [email protected] (ZCH); [email protected] (CDW)

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Chuduan Wang

* E-mail: [email protected] (ZCH); [email protected] (CDW)

Affiliation: National Engineering Laboratory for Animal Breeding and MOA Key Laboratory of Animal Genetics and Breeding, Department of Animal Genetics and Breeding, China Agricultural University, Beijing, 100193, China

Introduction

Pigs are an important source of meat worldwide [1]. The fatness traits of pigs impact growth rate, meat quality and reproductive performance [2].The thickness of backfat is a good indicator for fat deposition as it correlates well with over all body fat, carcass cross-sectional fat area ratios and intramuscular fat content[3]. Elucidating the role of specific genes in lipid metabolism and obesity is critically important and feasible to achieve in pigs. The similarity of pigs to humans in terms of body size and other physiological/anatomical features, including their innate tendency to over-consume food, has made pigs an important animal model for the study of the genetic basis of metabolic diseases such as obesity, type II diabetes, metabolic syndrome and atherosclerosis[4].

Adipose tissues, liver and skeletal muscles are the most important organs in animals with respect to whole body lipid metabolism [1,5]. Lipid is stored in adipose tissue to provide energy. Meanwhile the adipose tissue is the main location for de novo fatty acid synthesis, and is the major source of circulating free fatty acids (FFAs) in pigs[6]. Lipid in adipose tissue can be hydrolyzed to glycerol and FFAs. FFAs, combined with plasma albumin, can be transported and used as an energy source through oxidation [7]. Furthermore, it also produces adipocytokines, peptide hormones and resistin, and hence functions as an endocrine organ[8].

Technological advances now allow sequencing to be performed more economically and efficiently than ever before, providing excellent opportunities to investigate biological problems. The transcriptome of pig adipose tissue has recently been studied using RNA sequencing (RNA-seq). Differences in adipose tissue transcriptomes from different breeds and different growth phases were studied[9]. Functional differences were reported for the different types of adipose tissue[10],and adipose tissues in two phenotypically different individuals were compared for different fat deposition and fatty acid composition in muscle [11,12]. Meanwhile, whole genome resequencing has also been employed in livestock species such as cattle[13] and chicken[14].

Previous RNA-Seq studies to analyze the regulation of adipose deposition in pigs have either used only a very small number of adipose tissue samples or have lacked adequate control of genetic backgrounds. In our study, we applied RNA-Seq to obtain adipose tissue transcriptomes from two pair full-sibs and from one pair of unrelated pigs, the pigs in each pair having opposite backfat thickness phenotypes. We also used whole-genome resequencing technology to investigate the genetic basis for variations in expression among the differentially expressed genes.

Materials and Methods

Experimental design, animals and phenotypes

A Songliao black pig resource population (around age, 216 days which ranged from 211 to 218 days and average live weight, 103.9 kg which ranged 85.0 kg to 105.4 kg) was used in this study. All animals were housed in consistent and standard environmental conditions, given the same diet and ad libitum. Pedigree information was available for all animals. Live backfat thickness of 53 individuals was measured on the last 3rd and 4th ribs using the B mode real time ultrasound (HS1500, Honda, Japan) to choose pairs with divergent backfat phenotypes. Selection methods and standards were the same as in our previous study [15]. Based on our criteria, three pairs of pigs, two of which were full-sibs, were chosen with each pair having opposite backfat thickness phenotypes. All chosen pigs were stunned with a captive bolt, exsanguinated and slaughtered in commercial abattoir called Beijing Huadu Sunshine Food co., LTD. Slaughter house management gave the necessary permissions for slaughter. The chosen pigs were slaughtered according to guidelines of operating procedures of pig-slaughtering (GB/T 17236–2008), which was promulgated by General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China (AQSIQ) and Standardization Administration of the People’s Republic of China(SAC).All efforts to minimize animal suffering were made during the study. The whole procedure for collection of the tissue samples of all animals was by our researchers. This study was approved specifically by the Animal Welfare Committee of China Agricultural University (Permit number: DK996). Adipose samples from backfat between the last 3rd and 4th ribs were frozen in liquid nitrogen immediately after slaughter and stored at -80°C until used for RNA extraction. Liver samples were collected and stored at -20°C until used for DNA extraction.

Our experimental design is shown in Fig 1. We identified differentially expressed genes (DEGs) between Songliao Black pigs with higher backfat of live (BH) and BL Songliao Black pigs with lower backfat of live (BL) groups. We also used whole-genome resequencing to identify mutations in the two pig groups to identify and annotate genetic variation influencing fat deposition. We integrated transcriptome and genomic resequencing data to identify key genes and mutations affecting swine fat deposition.

[Figure omitted. See PDF.]

Fig 1. Experimental study design.

Schematic picture of the experimental design. Two of three pairs extreme individuals are full-sibs. BH: Songliao Black pigs with higher backfat of live group. BL:Songliao Black pigs with lower backfat of live group.

https://doi.org/10.1371/journal.pone.0122396.g001

mRNA Library construction and sequencing

One paired-end mRNA library was constructed from the adipose tissues of each of the six pigs. Purified mRNA was first fragmented using an RNA fragmentation kit (Ambion, Austin, TX, USA). Paired-end libraries were prepared following the Illumina paired-end library preparation protocol (Illumina, San Diego, California, USA). The libraries were sequenced using a multiplexed paired-ends protocol with 180bp of data collected per run on an Illumina HiSeq 2000 (Illumina). The average insert size for the paired-end libraries was 180 bp.

DNA library construction and sequencing

Sample preparation, cluster generation and sequencing were performed according to Illumina protocols. Total DNA from each pig was quantified, evenly mixed together for each group (BH and BL) and used for whole-genome resequencing. Purity and yield were checked using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, California, USA). Briefly, two multiplexed paired-end libraries were prepared and sequenced using a HiSeq2000 (Illumina). Average insert size for each library was 220 bp.

Bioinformatic analysis

The raw data were subjected todata quality evaluation, filtering and trimming before mapping on the swine reference genome. After the data quality control process, we used TopHat software for mapping short-reads and easyRNASeq software to quantify the raw reads. Differentially expressed genes (DEGs) were discovered using NOISeq [16]. NOISeq is a nonparametric approach for the identification of differentially expressed genes. Because of the noise distribution from the actual data, the NOISeq method could better adapt to the size of data set and was more effective in controlling the rate of false discoveries. A threshold of 0.8 was used for this probability, meaning that the gene is four times more likely to be differentially expressed than non-differentially expressed, as an adjusted P-value threshold of 0.001 for the other methods[17]. Functional enrichment analysis of the DEGs were performed in the DAVID system [18] and functional protein interaction networks were constructed in the STRING database[19]. Details of the RNA-Seq analysis are in our previous publication[15]. BWA software was used to map genomic resequencing data to the reference pig genome (Sscrofa10.2) with default parameters[20]. Quality score recalibration, local indel realignment and removal of duplicates from alignments was performed using GATK[21]. Single nucleotide variants (SNVs) and indel variants were called using GATK.

Quantitative PCR analysis

To validate changes in transcript levels between BH and BL groups (n = 3)(for genes with functions related to fatty acid synthesis), six genes, including ACACA, LDHA, ELOVL6, CYP1A2, PDK1 and SCD, were selected and quantified using qRT-PCR. Total RNA was extracted from adipose tissue and converted to cDNA using a Revert Aid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. The cDNA samples were then analyzed by RT-PCR using a Light Cycler 480 Real-Time PCR System (Roche, Hercules, CA, USA). The RT-PCR reactions were performed in a final volume of 20 μl using the Roche SYBR Green PCR Kit (Roche, Hercules, CA, USA) according to the manufacturer’s instructions. The RT-PCR primers designed for the six genes are listed in S1 Table. Gene expression values for qRT-PCR were normalized using the housekeeping gene, GAPDH. Triplicate qRT-PCRs were performed on each cDNA and the average Ct was used for further analysis. The relative quantification values were calculated using the 2-ΔΔCt method.

Results

Analysis of RNA Deep Sequencing Data

Similar to our previous study, live individuals in the BH group had twice the backfat thickness compared with that in the BL group[15].Other traits related to fat deposition, such as carcass backfat thickness and kidney fat weight, were also divergent in the BH and BL groups (Fig 2). Sequencing six samples yielded average 43.9 millions 90 bp paired-end reads and an around 80% of reads were mapped against the reference pig genome assembly, Sscrofa10.2, using TopHat software. 65.20%–72.90% of total tags were aligned to the CDS region of the reference pig genome. The ratio of reads mapped in the 3′-UTR, intron and 5′-UTR regions was diminished (Table 1).

[Figure omitted. See PDF.]

Fig 2. The traits of backfat thickness and kidney fat in BH and BL.

The three traits are backfat thickness of live animals (a), backfat thickness of carcass (b) and weight of kidney fat (c). BH and BL are pigs with higher/lower thickness backfat of live using B-Ultrasound respectively. *: significantly different (p < 0.05)

https://doi.org/10.1371/journal.pone.0122396.g002

[Figure omitted. See PDF.]

Table 1. Summary of RNA-seq alignment.

https://doi.org/10.1371/journal.pone.0122396.t001

Use the Ensembl V68 as the reference genome annotation to classify the mapping tags into the different regions. Ratio of the tags mapping on the subregion of the gene was calculated as the tags on each region divided by the total tags on the whole genome.

To normalize gene expression across different samples, the method of Trimmed Mean of M-values (TMM) was used to quantify transcript expression levels. The total number of expressed genes in different animals was similar and ranged from 17,175 to 17,715 (S2 Table). Correlations between biological replicate samples had very high reproducibility. The deep sequencing data of total RNA have been submitted to NCBI Sequence Read Archive (SRA) under Accession no.SRP035333, Bioproject: PRJNA234335.

To validate the RNA-seq results, on the basis of differential expression, a total of six genes with functions related to fatty acid synthesis, including ACACA, LDHA, ELOVL6, CYP1A2, PDK1 and SCD, were selected and their expression assayed byRT-PCR. The commonly used reference gene, GAPDH, was used for validation. For all six genes, the RT-PCR fold-change ratios between BH and BL groups were consistent with the RNA-Seq data (S1 Fig).

Differential gene expression and biological function analysis

In our analysis, we treated the high/low backfat thickness samples as the biological samples as they showed similar phenotypes in each group (BH or BL). A total of 188 DEGs, including 32 up-regulated and 156 down-regulated genes, were selected (Fig 3). These DEGs were catalogued according expression in BH and BL adipose tissue (S3 Table). Moreover, we also analyzed the DEGs between pigs in each pair and then identified the DEGs that were common among the three different pairs (Fig 4). The DEG details, including gene name, expression level and statistical information, were shown in the S4 Table. Because of the sample limitations of the DEGs for each pair, we only show DEGs and their functional analysis that were obtained by treating three pairs as the three biological replicates.

[Figure omitted. See PDF.]

Fig 3. Gene expression level in BH and BL groups.

A. X-axis plots gene expression counts in the BH group after TMM quantification and the y-axis plots gene expression counts in the BL group after TMM quantification. The red points indicate significantly differentially expressed genes.

https://doi.org/10.1371/journal.pone.0122396.g003

[Figure omitted. See PDF.]

Fig 4. Venn diagrams showing differentially expressed genes (DEGs) among each of three pairs.

A. Up-regulated genes in BH compared with BL in three pairs. B. Down-regulated genes in BH compared with BL in three pairs.

https://doi.org/10.1371/journal.pone.0122396.g004

About 27 of 28 highly expressed genes in the BH group were mapped to GO and KEGG pathways in DAVID. The functional enrichment of pathways (Table 2) and biological functions (Fig 5) included fatty acid biosynthesis, the insulin signaling pathway, biosynthesis of unsaturated fatty acids, lipid synthesis, fatty acid metabolic process, response to hormone, lipase activity, and saccharide metabolic process. These results showed that the enzymes involved in fatty acid and unsaturated fatty acid synthesis, such as ACACA, FASN, SCD and ELOVL6, were highly expressed in the BH group. Meanwhile, several genes in the insulin signaling pathway were significantly increased in the BL group. From 122 down regulated genes, 119 were assigned to gene ontology categories, such as molecular metabolism, receptor signaling pathway, immune response and lipid metabolism. Biological processes (Fig 5) and pathways (Table 2) relating to the metabolism of several compounds, including retinol, caffeine, arachidonic acid, and drugs were enriched in cytochrome P450 (CYP) family genes. The details of gene ontology categories were provided in the S5 Table.

[Figure omitted. See PDF.]

Fig 5. Heatmap of differentially expressed genes (DEGs) in adipose tissue and GO analysis.

The heatmap was generated by hierarchical cluster analysis of DEGs. The processes listed opposite the upper half of the heatmap are gene ontology annotation processes of the up-regulated genes [biological process (BP)]. The lower processes are those (BP) of the down-regulated genes.

https://doi.org/10.1371/journal.pone.0122396.g005

[Figure omitted. See PDF.]

Table 2. Pathways (P <0.05) enriched in differentially expressed genes (DEGs) in adipose tissue(BH vs. BL).

https://doi.org/10.1371/journal.pone.0122396.t002

STRING protein interaction network analysis generated a network from the genes that were up-regulated in BH compared with BL that were related to de novo fatty acid synthesis, from providing the substrate to participating in synthesis. Another network is related to immune action through three systems: 1) direct involvement in the antiviral process, 2) regulation of interfere on antiviral effectors, 3) functions in ubiquitination by immune and proinflammatory responses (Fig 6). Most genes in the immune action network are down-regulated in the BH group compared with the BL group.

[Figure omitted. See PDF.]

Fig 6. STRING analysis shows that DEGs are involved in known and predicted protein-protein interactions.

STRING is used to analyze of DEGs in adipose tissue between BH and BL. The network nodes stand for those genes shown in S3 Table. Lines of different colors represent seven types of evidence used in predicting associations. Red line: fusion evidence; green line: neighborhood evidence; blue line: co-occurrence evidence; purple line: experimental evidence; yellow line: text mining evidence; light blue line: database evidence and black line: co-expression evidence.

https://doi.org/10.1371/journal.pone.0122396.g006

Genome sequencing, SNP/indel detection, and functional annotation of genomic variation

Each whole genomic resequencing pooled sample were harvested to about 25G clean data with average 10×for each pool. Complete data sets of whole genomic resequencing have been submitted to NCBI SRA under Accession no. SRP040044, Bioproject: PRJNA240950. Using GATK, we identified single nucleotide polymorphisms (SNPs) and indels in BH and BL groups, some of which were common (Fig 7). The number of homozygous SNPs in BH and BL groups, was smaller than that of heterozygous SNPs, with a ratio of 1:1.46 and 1:1.36 (homozygous: heterozygous), respectively, but among the indels, the result was the opposite (2.63:1 and 2.61:1). Most of the SNPs and indels are located between genes or within intergenic regions and introns.

[Figure omitted. See PDF.]

Fig 7. The venn diagram of SNPs and Indels in BH and BL.

The number of SNPs and Indels in genome located in BH and BL groups.

https://doi.org/10.1371/journal.pone.0122396.g007

Because the two groups have opposite backfat thickness phenotypes, we focused on the SNPs and indels that were different between groups. The SNPs that were only in BH or BL were annotated to 21,134 and 21,623 genes, respectively, using the NCBI Reference sequence database (S6 Table). Non-synonymous (NS) SNPs and indels within a coding DNA sequence (cIndel) may affect gene expression or function. We detected 5,099 NS SNPs in 2753 genes and 145 cIndels in 104 genes in the BH group. Similarly, we detected 6,010 NS SNPs in 3116 genes and 223 cIndels in 139 genes in the BL group. In total, there were 2818 and 3198 genes with NS SNPs/cIndels in the BH and BL groups, respectively. After conversion to homologous human genes, we conducted a DAVID functional annotation clustering analysis of genes containing variants (NS SNPs/cIndels) to identify molecular functions (MF), biological processes (BP), cellular components (CC) and KEGG pathways (S7 Table). The genes with variants (NS SNPs/cIndels) were found in numerous pathways, including synthesis, metabolism, oxidation reduction, enzymes activity, molecule transport and binding. The most significant pathway was olfactory transduction in both BH and BL, and is involved with 107 and 110 genes in BH and BL, respectively. Sixty and 67 genes were detected in the lipid related process such as lipid biosynthesis, lipid transport, lipid localization, lipid binding, fatty acid metabolic process, and regulation of lipid kinase activity. NS SNPs/cIndels genes were significantly enriched in the linoleic acid metabolism, fatty acid metabolism and starch and sucrose metabolism which are very important in the lipid metabolism (S7 Table).

Mutations in DEGs supported by genomic resequencing

A large number of SNPs and indels were detected in adipose tissue DEGs in the BH and BL groups (S8 Table). The genetic variations included frame shift, splice site, UTR and non-synonymous coding lesions that may affect gene expression, function and even phenotype. As mentioned above, through DEG functional enrichment analysis, several genes regulating lipid processes were found. For these genes, genetic variation between BH and BL is particularly important. Almost all indels were found in introns except for two frame-shifts: one in the ACACA gene and one in the CYS2 3′-UTR. It is noteworthy that a great number of indels were identified in ACACA and ELOVL6 genes. Many SNPs were detected in important regions of the key DEGs. For instance, SCD, a key gene in the biosynthesis of unsaturated fatty acids, had seven SNPs in the 3′-UTR. For FASN, a key gene in the biosynthesis of fatty acids, the genetic variation was more complex and included three 5′-UTR and nine non-synonymous coding SNPs. SNPs were similarly identified in other genes, such as CYS2, LGALS13, ME1, ACLY and CYP2B22. Genetic variations in introns and in the 5kb upstream and downstream regions of key genes were also abundant. We also matched our DEGs’ position to the QTLs related to fat traits in the pig QTLdb (http://www.animalgenome.org/cgi-bin/QTLdb/SS/index). There are overlapping regions in 13 DEGs and QTLs related to fat traits (S9 Table), including the important gene ELOVL6. Those mutations in both DEGs and QTL region related to fat deposition would be studied in the next step work.

Discussion

Lipid metabolism is one of the hot research topics. We compared with several previous transcriptome studies [11,22–24] to narrow down candidate genes which are critical important for lipid metabolism. Based on these comparisons, we found 53 genes showed the consistent expression patterns across different pig breeds (S10 Table).

Fatty acid synthesis

Previous studies have shown that almost all fatty acid synthesis occurs in adipose tissue [6]. The lipid synthetases, including ACACA, FASN, SCD and ELOVL6 are up-regulated in BH. The ACACA gene encodes acetyl-CoA carboxylase alpha, which catalyzes the carboxylation of acetyl-CoA to malonyl-CoA. The acetyl-CoA carboxylase catalyzed reaction is the rate-limiting step and the key regulator of de novo fatty acid synthesis [25]. In past reports, the regulation of the ACACA transcript level depends on tissue, breed, and phenotype (fat versus lean)[25]. Consistent with our results, a significant up-regulation of the ACACA gene and of other genes involved in de novo lipogenesis was observed in muscle samples from fatter Duroc pigs in comparison to leaner ones[26]. FASN is also rate-limiting in de novo fatty acid synthesis as its main function is to catalyze the synthesis of palmitate from acetyl-CoA and malonyl-CoA[27]. Palmitic acid is normally the end point of de novo fatty acids biogenesis. SCD and ELOVL6 are very important enzymes involved in fatty acid desaturation and elongation[28]. SCD is the rate-limiting enzyme in the biosynthesis of monounsaturated fatty acids [29]. In human, a high SCD activity is associated with a metabolic state favoring hepatic triglyceride accumulation and expansion of adipose triglyceride stores[30]. These genes are enriched in the biosynthesis of fatty acids and unsaturated fatty acids. ME1 takes part in the tricarboxylic acid cycle to supply NADPH and to transport acetyl-CoA from mitochondria to the cytosol for the biosynthesis of fatty acids[31]. ME1, whose functions in glucose metabolism contribute to the initial steps of lipogenesis, is worth noting[32]. Consistent to previous study[11], ME1 mRNA was more abundant in the fatty compared with the leaner group, which is consistent with the biological function of this gene; promoting adipose deposition[33,34].ACLY is the primary enzyme responsible in many tissues for the synthesis of cytosolic acetyl-CoA from citrate and CoA [35]. Leptin is a fat-derived cytokine whose levels directly correlate with fat mass and communicate the energy status of the organism to the brain[36]and inhibit ACACA and/or FASN[37]. In our study, mRNA level of Leptin is more abundant in BH group, similarity a full-sib study [11]. This negative feedback regulates body weight through inhibition of feed intake, and increasing fatty acid and glucose oxidation[38]. Leptin also increases glucose use, and this is consistent with the function of the insulin signaling pathway which is enriched with up-regulated genes expressed in adipose tissue. The insulin signaling pathway is an important regulator of lipid metabolism through increasing glucose uptake in fat, promoting the storage of substrates in fat by stimulating lipogenesis, and inhibiting lipolysis [39]. THRSP is abundantly expressed in lipogenic tissues and plays an important role in the biosynthesis of triglycerides with a medium-length fatty acid chain and in modulating lipogenesis[40]. In mammals, biotin serves as a covalently bound coenzyme for acetyl-CoA carboxylases (ACACA) and participates directly in fatty acid biosynthesis[41]. Holocarboxylase synthetase (HLCS) catalyzes the binding of biotin to acetyl-CoA carboxylases (ACACA) and, therefore, plays a pivotal role in biotin-dependent metabolic pathways[42]. The above mentioned genes, which are highly expressed in adipose tissue of BH compared with BL animals, take part in the de novo synthesis of fatty acid in three aspects: 1) lipid synthetase, 2)supply or transport of substrate, 3)direct or indirect regulatory function (Fig 8). Fatty acid synthesis in adipose tissue may, therefore, decide the backfat thickness. Much effort has been expended to find the effect of polymorphisms in the above genes on gene expression and fat-related traits in animals. Polymorphisms in the ELOVL6 gene are associated with ELOVL6 expression levels and fatty acid synthesis in pigs[43]. ELOVL6 is up-regulated in BH pigs and we identified indels and SNPs within the gene. For ACACA, two SNPs have been associated with fatness and performance traits in Polish breeds and breed-specific differences in the transcript level were observed (fat versus lean) [25].Furthermore, we identified SNPs and indels in ACACA, which is up-regulated in our results. Five SNPs are located in the 3′-UTR of ME1 and are associated with the ME1 mRNA levels in muscular and adipose tissue and backfat thickness[44]. In our study, one ME13′-UTR SNP was found between BH and BL groups and the expression level of ME1 was different between BH and BL. This indicates that gene variation in key DEGs is very important. Because our groups are from a local Chinese breed of pig, the gene variation found in previous studies was not found in our study. Abundant genetic variations found in whole-genome resequencing provide a platform for more research. For fat deposition, synthesis of lipid is one key factor; lipid mobilization and fatty acid oxidation is another important area.

[Figure omitted. See PDF.]

Fig 8. The function of DEGs in the de novo fatty acids synthesis.

Black arrow means enzymes take part in the de novo fatty acids synthesis; red arrow means inhibiting effect; green arrow means promoting effect; blue quadrangle means up-regulation expression genes between BH and BL in adipose tissue

https://doi.org/10.1371/journal.pone.0122396.g008

Metabolism in adipose tissue

The adipocyte is not a passive lipid storage depot but a dynamic cell that plays a fundamental role in energy balance and overall body homeostasis[45]. The fat cell functions as a sensor of lipid levels, transmitting information to a neural circuit affecting hunger and satiety [46]. The PPAR signaling pathways regulate energy balance, lipid metabolism, cellular differentiation, and proliferation through PPARs, as members of the nuclear receptor superfamily[47,48].Activation of the PPAR signaling pathway promotes fatty acid oxidation[49].Genes in the PPAR signaling pathway and genes involved in fatty acid metabolism are down-regulated in adipose tissue between high and low backfat animals in cattle[50]. CYP2A13, in the PPAR signaling pathway, is down-regulated in adipose tissues between Korean native (fat) and Yorkshire (lean) pigs[51]. Consistent with our results, the PPAR signaling pathway and CYP2A13 are down-regulated in adipose tissue between BH and BL. This indicates a higher fatty acid metabolism in adipose tissue in the BL group. Retinoids (vitamin A) are crucial for most forms of life[52], and many studies have identified connections between retinoid and lipid metabolism[53]. Retinoids regulate metabolism by activating specific nuclear receptors, including the retinoic acid receptor (RAR) and the retinoid Xreceptor (RXR), an obligate heterodimeric partner for other nuclear receptors, including PPARs, and help to coordinate energy balance[54].Retinoid metabolism was down-regulated in adipose tissue and up-regulated in liver (BH vs. BL) and may cause the difference in lipid metabolism between the two groups. Meanwhile, the P450 gene family is activated by RXR and PPARs[55]. The P450 genes play critical roles in catalyzing reactions in the metabolism of drugs, environmental pollutants and other xenobiotics and in the oxidation of unsaturated fatty acids to intracellular messengers[56]. Several P450 genes were found at different expression levels between Korean native (fat) and Yorkshire (lean) pigs[51]. In our results, some P450 genes behave similarly.

Immune response and metabolic regulation in adipose tissue

The immune response and metabolic regulation are highly integrated. Adipose tissue is not an immune organ, but provides a close relationship between the immune system and metabolism[57].As above mentioned, the DEGs and pathways identified in this study regulate metabolism, especially fatty acid metabolism. Genes related to the immune system, such as NFKBIA, IRF1, MX1 and OAS1 show a significantly higher level of expression in the adipose tissue of BL pigs (the leaner group). Adipocytokines are a key component not only of metabolic regulation, but also of the immune response. Leptin, an abundant adipocyte product, represents one of the best examples of an adipocytokine. Leptin is involved in the control of energy expenditure, lipid and carbohydrate metabolism, which are mentioned above. Leptin also regulates endocrine and immune functions. Genetic defects of leptin or leptin receptors leads to reduced macrophage phagocytosis and the expression of proinflammatory cytokines, both in vivo and in vitro, and exogenous leptin up-regulates these[58]. These results indicate that leptin is up-regulated in inflammatory immune responses. But in our results, the mRNA level of leptin and of immune response-related genes was the opposite. Another way that metabolic and immune processes are regulated by lipids involves transcription factors, like PPARs and LXRs. Ligands to all three family members suppress NFKB, and then induce the production of proinflammatory cytokines. But unliganded PPARβ has an inflammatory function. In our results, both the PPAR signing pathway and NFKBIA are down-regulated in adipose tissue[59]. Studies demonstrate that changing the fatty acid composition of immune cells affects phagocytosis and other functions[60]. Candidate genes affecting fatty acid composition, such as SCD and ELOVL6, show different mRNA levels in adipose tissues between BH and BL. This is another example of the immune response and metabolism having countless ties. In our RNA-seq results, DEGs in adipose tissue included genes related to the immune response and metabolism. The specific interactions among them are not clear, but these data provide candidates for future studies.

Conclusions

Here we showed whole genome expression differences in adipose tissues and genetic variation among pigs with opposite backfat thickness phenotypes. RNA-seq provided a high resolution map of transcriptional activities and around 80% of the total reads could be mapped to annotated references. Animals are classified into two groups according to their backfat thickness and 188 genes were found to be differentially expressed between groups. The identification of numerous DEGs confirms the key pathways and gene networks related to lipid synthesis and metabolism. The DEGs such ACACA, FASN, SCD and ELOVL6 in de novo fatty acid synthesis in adipose tissue may affect backfat thickness. We also sequenced the whole genome of BH and BL mix pools. By comparison with the reference genome sequence, we also identified 4,493,002 SNPs and 68,276 indels in BH and 5,027,840 SNPs and 97,017 indels respectively. Genetic variation of this study, especially in DEGs, will provide valuable information for function studies, as well as for marker development associated with traits related to lipid deposition. Further studies are required to investigate the roles of candidate genes in fat deposition to improve pig breeding programs.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Comparison of qPCR and RNA-Seq expression ratios (BH and BL groups) for the selected genes.

https://doi.org/10.1371/journal.pone.0122396.s001

(TIF)

S1 Table. qRT-PCR primers for six genes related to lipid metabolism.

https://doi.org/10.1371/journal.pone.0122396.s002

(XLS)

S2 Table. Gene expression count in the six samples.

https://doi.org/10.1371/journal.pone.0122396.s003

(XLS)

S3 Table. DEGs between BH and BL in adipose tissue.

https://doi.org/10.1371/journal.pone.0122396.s004

(XLS)

S4 Table. DEGs common among the three different pairs.

https://doi.org/10.1371/journal.pone.0122396.s005

(XLSX)

S5 Table. Details of gene ontology categories for differently expression genes.

https://doi.org/10.1371/journal.pone.0122396.s006

(XLSX)

S6 Table. SNPs and indels only in BH and BL annotated in genes.

https://doi.org/10.1371/journal.pone.0122396.s007

(XLS)

S7 Table. Genes with NS or cIndel annotated in DAVID.

https://doi.org/10.1371/journal.pone.0122396.s008

(XLSX)

S8 Table. SNP and indels by type in differently expression genes.

https://doi.org/10.1371/journal.pone.0122396.s009

(XLSX)

S9 Table. Information of overlapping between DEGs and QTLs related to fat traits.

https://doi.org/10.1371/journal.pone.0122396.s010

(DOCX)

S10 Table. Detail information of common DEGs compared with other studies.

https://doi.org/10.1371/journal.pone.0122396.s011

(DOCX)

Acknowledgments

Grateful acknowledgments are made to the reviewers for their constructive suggestions for our manuscript. Authors are also indebted to Weiliang Zhou from Tianjin Ninghe primary pig breeding farm for providing pigs for experiment and to molecular quantitative genetics team members from China Agricultural University for experiment.

Author Contributions

Conceived and designed the experiments: ZCH CDW KX. Performed the experiments: KX ZJW. Analyzed the data: ZCH FZ. Contributed reagents/materials/analysis tools: ZJW HJL YW KX LWZ SKC. Wrote the paper: KX ZCH.

Citation: Xing K, Zhu F, Zhai L, Liu H, Wang Y, Wang Z, et al. (2015) Integration of Transcriptome and Whole Genomic Resequencing Data to Identify Key Genes Affecting Swine Fat Deposition. PLoS ONE 10(4): e0122396. https://doi.org/10.1371/journal.pone.0122396

References

1. Ramayo-Caldas Y, Mach N, Esteve-Codina A, Corominas J, Castello A, Ballester M, et al. (2012) Liver transcriptome profile in pigs with extreme phenotypes of intramuscular fatty acid composition. BMC Genomics 13: 547. pmid:23051667

2. Tummaruk P, Tantasuparuk W, Techakumphu M, Kunavongkrit A (2009) The association between growth rate, body weight, backfat thickness and age at first observed oestrus in crossbred Landrace x Yorkshire gilts. Anim Reprod Sci 110: 108–122. pmid:18289804

3. Suzuki K, Inomata K, Katoh K, Kadowaki H, Shibata T (2009) Genetic correlations among carcass cross-sectional fat area ratios, production traits, intramuscular fat, and serum leptin concentration in Duroc pigs. Journal of animal science 87: 2209–2215. pmid:19329483

4. Houpt KA, Houpt TR, Pond WG (1979) The pig as a model for the study of obesity and of control of food intake: a review. The Yale journal of biology and medicine 52: 307. pmid:380187

5. Fam BC, Joannides CN, Andrikopoulos S (2012) The liver: Key in regulating appetite and body weight. Adipocyte 1: 259–264. pmid:23700543

6. O'Hea EK, Leveille GA (1969) Significance of adipose tissue and liver as sites of fatty acid synthesis in the pig and the efficiency of utilization of various substrates for lipogenesis. J Nutr 99: 338–344. pmid:5350989

7. Nguyen P, Leray V, Diez M, Serisier S, Bloc’h JL, Siliart B, et al. (2008) Liver lipid metabolism. Journal of animal physiology and animal nutrition 92: 272–283. pmid:18477307

8. Cao H, Gerhold K, Mayers JR, Wiest MM, Watkins SM, Hotamisligil GS (2008) Identification of a lipokine, a lipid hormone linking adipose tissue to systemic metabolism. Cell 134: 933–944. pmid:18805087

9. Li X, Yang H, Li G, Zhang G, Cheng J, Guan H, et al. (2012) Transcriptome profile analysis of porcine adipose tissue by high‐throughput sequencing. Animal genetics 43: 144–152. pmid:22404350

10. Wang T, Jiang A, Guo Y, Tan Y, Tang G, Mai M, et al. (2013) Deep sequencing of the transcriptome reveals inflammatory features of porcine visceral adipose tissue. International journal of biological sciences 9: 550. pmid:23781149

11. Chen C, Ai H, Ren J, Li W, Li P, Qiao R, et al. (2011) A global view of porcine transcriptome in three tissues from a full-sib pair with extreme phenotypes in growth and fat deposition by paired-end RNA sequencing. BMC Genomics 12: 448. pmid:21906321

12. Ramayo-Caldas Y, Mach N, Esteve-Codina A, Corominas J, Castelló A, Ballester M, et al. (2012) Liver transcriptome profile in pigs with extreme phenotypes of intramuscular fatty acid composition. BMC Genomics 13: 547. pmid:23051667

13. Stothard P, Choi J-W, Basu U, Sumner-Thomson JM, Meng Y, Liao X, et al. (2011) Whole genome resequencing of black Angus and Holstein cattle for SNP and CNV discovery. BMC genomics 12: 559. pmid:22085807

14. Fan W-L, Ng CS, Chen C-F, Lu M-YJ, Chen Y-H, Liu CJ, et al. (2013) Genome-wide patterns of genetic variation in two domestic chickens. Genome biology and evolution 5: 1376–1392. pmid:23814129

15. Xing K, Zhu F, Zhai L, Liu H, Wang Z, Hou Z, et al. (2014) The liver transcriptome of two full-sibling Songliao black pigs with extreme differences in backfat thickness. Journal of Animal Science and Biotechnology 5: 32. pmid:25053997

16. Tarazona S, Furió-Tarı P, Ferrer A, Conesa A (2013) NOISeq: Differential Expression in RNA-seq.

17. Tarazona S, Garcia-Alcalde F, Dopazo J, Ferrer A, Conesa A (2011) Differential expression in RNA-seq: a matter of depth. Genome Res 21: 2213–2223. pmid:21903743

18. Huang da W, Sherman BT, Lempicki RA (2009) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4: 44–57. pmid:19131956

19. Franceschini A, Szklarczyk D, Frankild S, Kuhn M, Simonovic M, Roth A, et al. (2013) STRING v9. 1: protein-protein interaction networks, with increased coverage and integration. Nucleic acids research 41: D808–D815. pmid:23203871

20. Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25: 1754–1760. pmid:19451168

21. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome research 20: 1297–1303. pmid:20644199

22. Sodhi SS, Park WC, Ghosh M, Kim JN, Sharma N, Shin KY, et al. (2014) Comparative transcriptomic analysis to identify differentially expressed genes in fat tissue of adult Berkshire and Jeju Native Pig using RNA-seq. Mol Biol Rep 41: 6305–6315. pmid:25008993

23. Moon JK, Kim KS, Kim JJ, Choi BH, Cho BW, Kim TH, et al. (2009) Differentially expressed transcripts in adipose tissue between Korean native pig and Yorkshire breeds. Anim Genet 40: 115–118. pmid:18945290

24. Li X, Yang H, Li G, Zhang G, Cheng J, Guan H, et al. (2012) Transcriptome profile analysis of porcine adipose tissue by high‐throughput sequencing. Anim Genet 43: 144–152. pmid:22404350

25. Stachowiak M, Nowacka-Woszuk J, Szydlowski M, Switonski M (2013) The ACACA and SREBF1 genes are promising markers for pig carcass and performance traits, but not for fatty acid content in the longissimus dorsi muscle and adipose tissue. Meat Sci 95: 64–71. pmid:23657179

26. Canovas A, Quintanilla R, Amills M, Pena RN (2010) Muscle transcriptomic profiles in pigs with divergent phenotypes for fatness traits. BMC Genomics 11: 372. pmid:20540717

27. Wakil SJ (1989) Fatty acid synthase, a proficient multifunctional enzyme. Biochemistry 28: 4523–4530. pmid:2669958

28. Jakobsson A, Westerberg R, Jacobsson A (2006) Fatty acid elongases in mammals: their regulation and roles in metabolism. Progress in lipid research 45: 237–249. pmid:16564093

29. Doran O, Moule SK, Teye GA, Whittington FM, Hallett KG, Wood JD (2006) A reduced protein diet induces stearoyl-CoA desaturase protein expression in pig muscle but not in subcutaneous adipose tissue: relationship with intramuscular lipid formation. Br J Nutr 95: 609–617. pmid:16512947

30. Flowers MT, Ntambi JM (2008) Role of stearoyl-coenzyme A desaturase in regulating lipid metabolism. Curr Opin Lipidol 19: 248–256. pmid:18460915

31. Wise EM Jr., Ball EG (1964) Malic Enzyme and Lipogenesis. Proc Natl Acad Sci U S A 52: 1255–1263. pmid:14231450

32. Corominas J, Ramayo-Caldas Y, Puig-Oliveras A, Estellé J, Castelló A, Alves E, et al. (2013) Analysis of porcine adipose tissue transcriptome reveals differences in de novo fatty acid synthesis in pigs with divergent muscle fatty acid composition. BMC genomics 14: 843. pmid:24289474

33. Zhou SL, Li MZ, Li QH, Guan JQ, Li XW (2012) Differential expression analysis of porcine MDH1, MDH2 and ME1 genes in adipose tissues. Genet Mol Res 11: 1254–1259. pmid:22614353

34. Schmid GM, Converset V, Walter N, Sennitt MV, Leung KY, Byers H, et al. (2004) Effect of high-fat diet on the expression of proteins in muscle, adipose tissues, and liver of C57BL/6 mice. Proteomics 4: 2270–2282. pmid:15274121

35. Ren Z-Q, Wang Y, Xu Y-J, Wang L-J, Lei M-G, Zuo B, et al. (2008) Identification of a differentially expressed gene, ACL, between Meishan× Large White and Large White× Meishan F1 hybrids and their parents. Genet Sel Evol 40: 625–637. pmid:18990355

36. Zhang Y, Proenca R, Maffei M, Barone M, Leopold L, Friedman JM (1994) Positional cloning of the mouse obese gene and its human homologue. Nature 372: 425–432. pmid:7984236

37. Lage R, Dieguez C, Vidal-Puig A, Lopez M (2008) AMPK: a metabolic gauge regulating whole-body energy homeostasis. Trends Mol Med 14: 539–549. pmid:18977694

38. Dieguez C, Vazquez MJ, Romero A, Lopez M, Nogueiras R (2011) Hypothalamic control of lipid metabolism: focus on leptin, ghrelin and melanocortins. Neuroendocrinology 94: 1–11. pmid:22143028

39. Saltiel AR, Kahn CR (2001) Insulin signalling and the regulation of glucose and lipid metabolism. Nature 414: 799–806. pmid:11742412

40. Colbert CL, Kim C-W, Moon Y-A, Henry L, Palnitkar M, McKean WB, et al. (2010) Crystal structure of Spot 14, a modulator of fatty acid synthesis. Proceedings of the National Academy of Sciences 107: 18820–18825. pmid:20952656

41. Wakil SJ, Gibson DM (1960) Studies on the mechanism of fatty acid synthesis VIII. The participation of protein-bound biotin in the biosynthesis of fatty acids. Biochimica et Biophysica Acta 41: 122–129.

42. Zempleni J, Kuroishi T (2012) Biotin. Adv Nutr 3: 213–214. pmid:22516729

43. Corominas J, Ramayo-Caldas Y, Puig-Oliveras A, Pérez-Montarelo D, Noguera JL, Folch JM, et al. (2013) Polymorphism in the ELOVL6 gene is associated with a major QTL effect on fatty acid composition in pigs. PLoS ONE 8: e53687. pmid:23341976

44. Vidal O, Varona L, Oliver M, Noguera J, Sanchez A, Amills M (2006) Malic enzyme 1 genotype is associated with backfat thickness and meat quality traits in pigs. Animal genetics 37: 28–32. pmid:16887000

45. Rosen ED, Spiegelman BM (2006) Adipocytes as regulators of energy balance and glucose homeostasis. Nature 444: 847–853. pmid:17167472

46. Kopecky J, Rossmeisl M, Flachs P, Brauner P, Sponarova J, Matejková O, et al. (2004) Energy metabolism of adipose tissue—physiological aspects and target in obesity treatment. Physiol Res 53 Suppl 1: S225–232. pmid:15119952

47. Lemberger T, Desvergne B, Wahli W (1996) Peroxisome proliferator-activated receptors: a nuclear receptor signaling pathway in lipid physiology. Annu Rev Cell Dev Biol 12: 335–363. pmid:8970730

48. Evans RM, Barish GD, Wang YX (2004) PPARs and the complex journey to obesity. Nat Med 10: 355–361. pmid:15057233

49. Roman J (2008) Peroxisome proliferator-activated receptor gamma and lung cancer biology: implications for therapy. J Investig Med 56: 528–533. pmid:18317436

50. Taniguchi M, Guan LL, Basarab JA, Dodson MV, Moore SS (2008) Comparative analysis on gene expression profiles in cattle subcutaneous fat tissues. Comparative Biochemistry and Physiology Part D: Genomics and Proteomics 3: 251–256. pmid:20494844

51. Choi K-M, Moon J-K, Choi S-H, Kim K-S, Choi Y-I, Kim J-J, et al. (2008) Differential expression of cytochrome P450 genes regulate the level of adipose arachidonic acid in Sus Scrofa. Asian-Aust J Anim Sci 21: 967–971.

52. Blomhoff R, Blomhoff HK (2006) Overview of retinoid metabolism and function. J Neurobiol 66: 606–630. pmid:16688755

53. Keller H, Dreyer C, Medin J, Mahfoudi A, Ozato K, Wahli W, et al. (1993) Fatty acids and retinoids control lipid metabolism through activation of peroxisome proliferator-activated receptor-retinoid X receptor heterodimers. Proc Natl Acad Sci U S A 90: 2160–2164. pmid:8384714

54. Ziouzenkova O, Plutzky J (2008) Retinoid metabolism and nuclear receptor responses: New insights into coordinated regulation of the PPAR-RXR complex. FEBS Lett 582: 32–38. pmid:18068127

55. Waxman DJ (1999) P450 gene induction by structurally diverse xenochemicals: central role of nuclear receptors CAR, PXR, and PPAR. Archives of Biochemistry and Biophysics 369: 11–23. pmid:10462436

56. Lewis DF (2004) 57 varieties: the human cytochromes P450. Pharmacogenomics 5: 305–318. pmid:15102545

57. Hotamisligil GS (2006) Inflammation and metabolic disorders. Nature 444: 860–867. pmid:17167474

58. Loffreda S, Yang S, Lin H, Karp C, Brengman M, Wang DJ, et al. (1998) Leptin regulates proinflammatory immune responses. The FASEB Journal 12: 57–65. pmid:9438411

59. Glass CK, Ogawa S (2006) Combinatorial roles of nuclear receptors in inflammation and immunity. Nat Rev Immunol 6: 44–55. pmid:16493426

60. Calder PC (2008) The relationship between the fatty acid composition of immune cells and their function. Prostaglandins, Leukotrienes and Essential Fatty Acids 79: 101–108. pmid:18951005