RNA-seq analysis of subcutaneous adipocytes from Jiaxing black pig and Large white pig
Primary cultures of SC adipocytes were isolated from 3-day-old Jiaxing Black (JX) and Large White (LW) pigs, and subjected eight days of differentiation. The differentiated adipocytes were harvested and subjected to RNA-seq analysis in three biological replicates. The Illumina HiSeq platform provided an average of 15.2 GB of clean reads for each sample. The percentage of clean reads among the raw data in each library ranged from 91.444–95.103%. For each sample, 90.22%, 87.20%, 86.05%, 84.79%, 85.59% and 85.76% unique reads could be mapped to the current version of the pig genome (Sscrofa 11.1), representing 23290, 22269, 22597, 22201, 21809 and 21981 genes, respectively (Supplementary Table 1). Gene numbers within a certain range of FPKM values were analyzed, and each sample gave similar results (Supplementary Fig. 1A). The distribution of log10 transformed FPKM values indicated the expression levels of transcripts. In each breed, the mRNA expression levels were relatively higher than those of lncRNAs, as expected, while both mRNAs and lncRNAs showed similar distribution in both breeds (Supplementary Figs. 1B and 1C). Taken together, both the biological replicates and sequencing data indicated sufficiently good quality for further analysis.
Differentially expressed lncRNAs and genes in subcutaneous adipocytes from the two breeds
To further understand the differences of SC adipocytes between the two breeds, comparative transcriptome analysis was conducted. A total of 4553 differentially expressed RNAs were identified, including 1258 up- and 3295 downregulated RNAs in SC adipocytes between the two breeds (Fig. 1A). Among these differentially expressed RNAs, 3371 were mRNAs (DEGs) and 1182 were lncRNAs (DELs). Among the 3371 DEGs, 797 were up- and 2574 were downregulated. Among the 1182 DELs, 461 were up- and 721 were downregulated (Fig. 1B). Among these differentially expressed mRNAs and lncRNAs, 2464 differentially expressed genes (DEGs) and 429 DELs had been previously annotated, and 907 DEGs and 753 DELs were novel (Fig. 1C).
Functional enrichment analysis of differentially expressed genes
The potential functions and signaling pathways of the differentially expressed genes (DEGs) were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. GO analysis based on biological process was conducted and the top 20 most highly enriched categories were listed (Fig. 2A). The results showed that DEGs related to cell proliferation, differentiation, migration, adhesion, apoptosis and cell-cell signaling were significantly enriched. Four processes related to immune response, namely “response to lipopolysaccharide”, “inflammatory response”, “immune response”, and “defense response to virus”, were also detected. Two terms closely associated with lipid metabolism were also identified. These were “calcium ion binding” and “positive regulation of ERK1 and ERK2 cascade”. In addition, the enrichment based on KEGG analysis was performed and the top 20 pathways are presented in Fig. 2B. Among these results, several immune-related pathways also appeared, such as “Chemokine signaling pathway”, “Systemic lupus erythematosus”, “Tuberculosis”, “Cytokine-cytokine receptor interaction”, “Leukocyte transendothelial migration”, “Phagosome” and “Rheumatoid arthritis”. Moreover, “Ras signaling pathway”, “PI3K-Akt signaling pathway”, “Osteoclast differentiation”, and “MAPK signaling pathway” were significantly enriched, all of which are highly associated with adipocyte differentiation and lipid accumulation. The results of enrichment analysis illustrated the regulatory differences of SC fat deposition between JX and LW pigs.
Protein-protein interaction network analysis
In unsupervised hierarchical clustering analysis, heat maps were generated using the differentially expressed DEGs, and they clearly self-segregated into different clusters for JX and LW pigs. These results reflected the distinct mRNA expression profiles of the two breeds (Fig. 3A). The protein-protein network was constructed using the top 20 DEGs ranked using the Maximal Clique Centrality (MCC) topological algorithm. As can be seen in the network shown in Fig. 3B, six DEGs were closely related to the distinctions of SC fat differentiation in both breeds. Compared to LW pigs, three DEGs (MMP9, MKI67 and VCL) were up- and three (SPTAN1, TLR2 and KIT) were downregulated in differentiated subcutaneous adipocytes of JX pigs (Fig. 3C).
Functional enrichment analysis of lncRNAs based on target genes
Based on RNA-seq data, the potential target genes of DELs were predicted to explore their potential functions. In GO enrichment analysis, the two immunity-related terms “positive regulation of I-κB kinase/NFκB signaling” and “innate immune response”, were detected. The categories cell proliferation, apoptosis and cell adhesion, which were detected in the DGE analysis, were also identified here (Fig. 4A). In KEGG enrichment analysis, four pathways related to adipocyte differentiation and lipid accumulation were identified, including “PI3K-Akt signaling pathway”, “cGMP-PKG signaling pathway”, “MAPK signaling pathway” and “Calcium signaling pathway” (Fig. 4B). In addition, the protein-protein network was constructed to offer new insights into the related biological processes. Compared to expression levels in LW pigs, a quarter of the target genes were up- and three quarters were downregulated in JX pigs (Fig. 4C). Among the interacting genes, Rab7a, WDR12, LPAR1, TBX5, Dicer1, WEE1, CDC25B, WDR12, CAPZB, UVRAG, and ST6Gal-I, are known to participate in the regulation of cell proliferation or differentiation. Furthermore, Dicer1, TBX5, TMEM173, PRPF8, AKAP9, ZBP1, UVRAG, ST6Gal-I, LRRFIP2, and HDAC10 are known to mediate the immune response. Moreover, LPAR1, LPAR1and AKAP9 were also reported to be associated with the PI3K/AKT signaling pathway or MAPK signaling pathway. Significantly, Nfe2l1and PLAG1 were reported to have an impact on the plasticity of adipose tissue. Thus, the DELs might play an essential role in the distinct adipogenesis of JX pigs.
Validation of the differentially expressed genes and lncRNAs
To validate the reliability of the RNA-seq results, 12 differentially expressed genes (DEGs) and 10 lncRNAs (DELs) were randomly chosen for quantitative PCR (qPCR) verification (Fig. 5A-5D). Compared with the RNA-seq data, 10 DEGs and 8 DELs gave consistent results, while two DEGs (MGP and RESTI) and two DELs (aXR_002337668.1 and LTCONS_00084076) showed statistically different results by qPCR analysis. Overall, 81.8% of the results were in agreement between the two techniques.
Verification of the pathway analysis
Because the MAPK signaling pathway was identified in the functional enrichment analysis of both DEGs and DELs, its activity in SC fat tissues of the two breeds was examined. Expression levels and phosphorylation of two kinases in the MAPK pathway, ERK1/2 and p38, were determined by western blot analysis. ERK1/2 showed no differences in total protein abundance or phosphorylation between the two breeds. However, while the total protein abundance of p38 was similar, the abundance of phosphorylated p38 showed a remarkably large difference between the two breeds. While p38 phosphorylation was practically absent in LW pig samples, it presented a heavy band in the JX pig samples (Fig. 6). The difference of p38 phosphorylation indicated that the activity of the MAPK pathway in SC fat tissue varied significantly between the two breeds, which confirmed the results of pathway enrichment analysis.