Research | Open | Published:
A genome-wide landscape of mRNAs, lncRNAs, and circRNAs during subcutaneous adipogenesis in pigs
Journal of Animal Science and Biotechnologyvolume 9, Article number: 76 (2018)
Preadipocyte differentiation plays a critical role in subcutaneous fat deposition in pigs. However, the roles of different RNAs, such as messenger RNAs (mRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) in the differentiation process of subcutaneous preadipocytes, are still largely unclear. In the present study, a transcriptome analysis, including the analysis of mRNAs, lncRNAs, and circRNAs, during different differentiation stages, namely, day 0 (D0), day 2 (D2), day 4 (D4), and day 8 (D8), of subcutaneous preadipocytes from Chinese Erhualian pigs was performed.
A total of 1554, 470, 1344, 1777, and 676 differentially expressed (DE) mRNAs, 112, 58, 95, 136, and 93 DE lncRNAs, and 902, 787, 710, 932, and 850 DE circRNAs were identified between D2 and D0, between D4 and D2, between D8 and D4, between D4 and D0, and between D8 and D0, respectively. Furthermore, functional enrichment analysis showed that the common DE mRNAs during the entire differentiation process were mainly involved in lipid metabolic and cell differentiation processes. Additionally, co-expression network analysis identified the potential lncRNAs related to adipogenesis, e.g., MSTRG.131380 and MSTRG.62128.
Our study provides new insights of the expression changes of RNAs during adipogenic differentiation, which might contribute to the phenotype of subcutaneous adipogenesis. These results greatly improve our understanding of the molecular mechanisms regulating subcutaneous fat deposition in pigs.
Being a crucial part of body energy metabolism and the endocrine system, subcutaneous fat plays an important role in the growth and meat quality of pigs. On one hand, excessive subcutaneous fat deposition in pigs, especially in obese pig breeds, greatly decreases the growth performance and production efficiency, which results in high costs of feeding and production. On the other hand, a sufficient layer of subcutaneous fat is necessary to obtain high quality processed products such as dry-cured hams [1, 2]. Additionally, the thickness of subcutaneous fat has been reported to positively correlate with intramuscular fat (IMF) in the longissimus dorsi and gluteus medius muscles . Therefore, understanding the mechanism of subcutaneous fat formation would be greatly beneficial to improve the production efficiency and meat quality. Because the increase in the size of porcine adipocytes is closely related with adipose tissue expansion, controlling subcutaneous preadipocyte differentiation could be considered a good strategy for regulating subcutaneous fat management.
In recent years, two main types of noncoding RNAs (ncRNAs), including long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), have been regarded as key regulators, because they play a vital role in regulating gene expressions at the transcriptional and post-transcriptional levels [4, 5]. LncRNAs are one class of ncRNAs that are longer than 200 nucleotides in length, but contain no open reading frames ; they have been revealed to affect preadipocyte differentiation by regulating adipogenic-related key genes such as PPARG and CEBPA in humans and mice [7,8,9], demonstrating that lncRNAs might have an essential role in adipogenesis. CircRNAs are a unique class of ncRNAs with a covalently closed continuous loop without 5′ caps and 3′ tails; they have been shown to be widely expressed in animal tissues and cells . More interestingly, certain circRNAs have tissue-specific and stage-specific expression patterns , indicating that circRNAs would be a specific type of regulator in cell or tissue development processes. Additionally, accumulating evidences indicate that circRNAs have an important role in mammalian cell differentiation [12, 13]. However, thus far, the functions of lncRNAs and circRNAs in porcine subcutaneous preadipocyte differentiation are still largely unknown.
Previously, several studies discovered the messenger RNA (mRNA) or microRNA (miRNA) expression profiles of subcutaneous tissues between high and low backfat pigs [14, 15]. Furthermore, the changes of mRNAs in subcutaneous preadipocytes during adipogenic differentiation were studied in pigs and mice [16, 17]. Moreover, lncRNAs related to castration-induced subcutaneous fat changes were identified in Huainan male pigs . However, the expression profiling of RNAs during porcine subcutaneous preadipocyte differentiation has not yet been well studied. For example, the mRNA transcriptome analysis of porcine subcutaneous preadipocytes during their differentiation is only performed at the early and middle stages of differentiation in Large White pigs . Furthermore, little is known about the expression characteristics of lncRNAs and circRNAs during adipogenic differentiation. Accordingly, it is necessary to further analyze the expression patterns of RNAs, including coding and noncoding RNAs, during porcine subcutaneous preadipocyte differentiation. Additionally, as a typical indigenous pig breed with plenty of subcutaneous fat, the Chinese Erhualian pig is a good model for studying subcutaneous fat formation. As such, the expression characters of mRNAs, lncRNAs, and circRNAs during different differentiation stages (day 0 (D0), day 2 (D2), day 4 (D4), and day 8 (D8)) of subcutaneous preadipocytes in Erhualian pigs were investigated using RNA sequencing (RNA-Seq) technology. Our results demonstrate the genome-wide changes of molecular events during adipogenic differentiation, thus giving us newer insights regarding subcutaneous fat management of pigs.
Three five-day-old Chinese Erhualian piglets were purchased from Changzhou Erhualian Pig Production Cooperation (Changzhou, Jiangsu, China). The experimental procedures were approved by the Institutional Animal Care and Use Committee of Nanjing Agricultural University.
Preadipocyte culture and differentiation
Newly isolated subcutaneous adipose tissue from each piglet was washed thrice with phosphate-buffered saline (PBS). Then, the tissue was minced and digested with 1 mg/mL collagenase type I (Invitrogen, Carlsbad, CA, USA) at 37 °C for 60 min, followed by filtration through a 200 μm nylon mesh for removing the undigested fractions. The filtrated solution was centrifuged at 1,000 r/min for 10 min to collect the preadipocytes, and then, the cells were cultured in Dulbecco’s modified Eagle’s medium/Ham’s F-12 (DMEM-F12) growth medium containing 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin at 37 °C with an atmosphere of 5% CO2. After the preadipocytes reached confluence (D0), the DMI hormone cocktail (1 μmol/L dexamethasone, 0.5 mmol/L 3-isobutyl-1-methylxanthine, and 5 μg/mL insulin) was added to the growth medium to induce the cell differentiation for 2 d. Next, the cells were subjected to maintenance medium (growth medium supplemented with 5 μg/mL insulin) for an additional 2 d. After that, the growth medium was changed every alternate day until adipocyte maturation.
Oil Red O staining and triglyceride assay
After removing the culture medium, the adipocytes were washed thrice with PBS and fixed in 10% formaldehyde for 15 min. Next, the cells were washed thrice with PBS, and then stained with Oil Red O for 20 min. Finally, the cells were washed thrice with PBS and photographed using an inverted microscope (Leica, Wetzlar, Germany). The absorbance values of Oil Red O-stained cells were measured at the wavelength of 510 nm to quantify the lipid accumulation. Meanwhile, triglyceride contents were determined using a commercial triglyceride assay kit (Applygen, Beijing, China), according to the manufacturer’s protocol.
RNA extraction, library preparation, and sequencing
Total RNA was isolated at each time point (D0, D2, D4, and D8) using the Trizol reagent (Invitrogen, Carlsbad, CA, USA). The qualities and quantities of the RNA were measured using Bioanalyzer 2100 (Agilent Technologies, Santa Clare, CA, USA) and 1% agarose gel electrophoresis, which showed that the RNA integrity number (RIN) values of all samples were 10. Ribosomal RNA from each sample was removed using the Ribo-Zero™ GoldKits (Epicentre, Madison, WI, USA). Equal amounts of total RNA of the same stage of differentiation from three Erhualian piglets were pooled into one sample. Then, cDNA libraries were prepared using a NEB Next Ultra Directional RNA LibraryPrep Kit (NEB, Ispawich, MA, USA), according the manufacturer’s instructions and sequenced using the Illumina HiSeq X Ten system (Illumina, San Diego, CA, USA).
Identification of lncRNA and circRNA
Low-quality and adaptor-polluted reads were firstly removed from the raw data. The clean reads from each sample were aligned onto the pig reference genome (Sus scrofa 11.1) using the HiSAT2 v2.0.5 program . In addition, because the sequences of circRNA cannot be directly aligned to the reference genome, the slicing alignment was mapped to the genome for obtaining the circRNA using the Burrows-Wheeler Aligner-maximal exact match (BWA-MEM) algorithm . The candidate lncRNAs were identified using the following criteria: 1) transcripts were filtered by removing those shorter than 200 nt, those with less than two exons, and those with a read coverage less than five in all samples, to avoid unreliable transcripts or those with inconsistent lengths, and sequences consisting of the known mRNAs and other non-coding RNAs (ribosomal RNA (rRNA), transfer RNA (tRNA), small nucleolar RNA (snoRNA), and small nuclear RNA (snRNA)); 2) putative lncRNAs for their protein-coding ability were determined using four approaches, including coding-non-coding index (CNCI), coding potential calculator (CPC), protein folding domain database (PFAM), and coding potential assessing tool (CPAT). Finally, the remaining transcripts were defined as novel lncRNAs. Furthermore, the candidate circRNAs were recognized using the CIRI (circRNA identifier) algorithm . In brief, paired chiastic clipping, paired-end mapping, and GT-AG splicing signals were discovered via scanning the above obtained slicing alignments. Next, the alignment files were scanned again using a dynamic programming algorithm for detecting additional junction reads and eliminating false positive circRNA candidates. The final circRNAs were obtained by retaining sequences with ≥2 junction reads.
Analysis of differentially expressed (DE) genes
The fragments per kilobase of transcript per million reads mapped (FPKM) value was used to estimate the expression levels of mRNAs and lncRNAs, while the spliced reads per billion mapping (SRPBM) value was utilized to determine the amount of circRNAs . Genes with an FPKM or SRPBM value of < 1 in no less than 50% of samples were defined as unreliably expressed genes, while those with an FPKM or SRPBM value of ≥1 in more than 50% of samples were considered as reliably expressed genes. Differentially expressed (DE) genes including mRNAs, lncRNAs, and circRNAs were analyzed using DEGseq v1.18.0 , which defined DE genes as reliably expressed genes with |log2_ratio| ≥ 1 and false discovery rate (FDR) < 0.05 between any two groups. Meanwhile, genes differentially expressed in three comparisons (D2 versus D0, D4 versus D0, and D8 versus D0) were defined as common DE genes.
Gene ontology (GO) analysis
The function of DE lncRNAs was predicted by the GO analysis of their cis- and trans-target mRNAs, which were screened based on their genomic positional relation 50 kilobase pairs (kb) upstream and 50 kb downstream, for cis-target mRNAs and based on the Pearson correlation coefficient of lncRNA-RNA pairs being ≥0.9, for trans-target mRNAs. The function of DE circRNAs was revealed via GO analysis of their parental genes. In brief, all genes were firstly mapped to GO terms using the Gene Ontology database (http://www.geneontology.org), and then, the functional enrichment analysis was performed using Omicshare tools (www.omicshare.com/tools). Terms with P-values less than 0.01 and more than 3 matching genes were identified as significant or enriched terms.
Co-expression network construction
The co-expression network of common DE lncRNAs with their cis- and trans-target common DE mRNAs was constructed using the Cytoscape software (v3.2.1) to explore the function of key lncRNAs .
Quantitative real-time RT-PCR
Primers used in quantitative real-time RT-PCR (qRT-PCR) were listed in Additional file 1. Three RNAs samples per differentiation stage were reverse transcribed using a PrimeScript™ Master Mix (Takara, Dalian, China), according to the manufacturer’s instructions. Next, qRT-PCR was performed using SYBR Premix Ex Taq™ (Takara, Dalian, China) on the StepOnePlus™ Real-time PCR System (Applied Biosystems, Foster City, CA, USA). The reaction conditions were as follows: denaturation for 30 s at 95 °C, followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. Meanwhile, RPLP0 (ribosomal protein lateral stalk subunit P0), PPIA (peptidylprolyl isomerase A), and HPRT1 (hypoxanthine phosphoribosyltransferase 1) were used as the normalization controls, and the experiments were performed in triplicate. The 2-ΔΔCT method was used to calculate the relative gene expression levels.
Statistical analysis of the data from triglyceride and qRT-PCR assay was performed using the SPSS software (version 20.0). Statistical comparison among the groups was analyzed using one-way analysis of variance (ANOVA), followed by Tukey’s multiple comparison test. P-values less than 0.05 were considered significant, while P-values less than 0.01 or 0.001 were considered highly significant.
The sequencing data have been submitted to the Gene Expression Omnibus (GEO) database (accession number GSE114583).
Phenotypic changes during preadipocyte differentiation
Compared to the cell shapes in the initial phase (D0), the preadipocytes gradually turned from fibrous into a spherical shape at D2. Subsequently, lipid droplets were visibly observed at D4 and gradually increased until D8 (Fig. 1a). In accordance with the results of the adipocyte shapes, the triglyceride levels progressively accumulated, accompanying the increase of the differentiation process (Fig. 1b-c). These data indicate that the subcutaneous preadipocyte differentiation process is normal and can be further analyzed.
Characters of RNA-seq
After quality control, a total of 155981654, 137762268, 168562368, and 119115476 clean reads with greater than 94.80% of Q30 were obtained in D0, D2, D4, and D8, respectively (Table 1). Among them, a total of 97.49%, 97.60%, 97.66%, and 97.88% reads from D0, D2, D4, and D8, respectively, were mapped to the pig reference genome (Sus scrofa 11.1) (Table 1). Additionally, BWA-MEM was used for sequence split comparison to accurately identify circRNAs; a mapping rate of more than 99.91% from each sample was discovered (Table 1).
Gene expression profiles during preadipocyte differentiation
A total of 16192 mRNAs, 5615 lncRNAs, and 8623 circRNAs were obtained from four stages (D0, D2, D4, and D8) of differentiation. LincRNAs and classic circRNAs accounted for the maximum proportion of novel lncRNAs and circRNAs, respectively (Fig. 2a-b). Based on an FPKM or SRPBM value of ≥1 in more than 2 groups, 9191 mRNAs, 430 lncRNAs (including 208 annotated lncRNAs and 222 novel lncRNAs), and 2172 circRNAs were identified as reliably expressed genes (Fig. 2c and Additional file 2). Furthermore, five comparison groups (D2 vs D0, D4 vs D2, D8 vs D4, D4 vs D0, and D8 vs D0) were set to discover the differentially expressed (DE) genes during the differentiation. A total of 2758 mRNAs, 230 lncRNAs, and 1760 circRNAs were defined as DE genes among these comparison groups (Fig. 3a-c). Of these, 1554, 470, 1344, 1777, and 676 DE mRNAs, 112, 58, 95, 136, and 93 DE lncRNAs, and 902, 787, 710, 932, and 850 DE circRNAs were identified between D2 and D0, between D4 and D2, between D8 and D4, between D4 and D0, and between D8 and D0, respectively (Fig. 3d and Additional files 3, 4, 5, 6 and 7).
Gene ontology (GO) analysis of DE genes between D2 and D0
Between D2 and D0, significantly upregulated GO terms were mainly involved in: 1) lipid metabolism-related processes, response to oxygen-containing compound, and organic acid metabolic process for DE mRNAs, 2) carboxylic acid metabolic process, regulation of developmental growth, and small molecule catabolic process for DE lncRNAs, and 3) chromatin modification, histone methylation, and regulation of cell communication for DE circRNAs (Fig. 4a and Additional file 8). Additionally, significantly downregulated GO terms between D2 and D0 were mainly related to: 1) cytoskeleton organization and cell cycle-related processes for DE mRNAs, 2) positive regulation of developmental growth and cell cycle-related processes for DE lncRNAs, and 3) vesicle-mediated transport, developmental cell growth, and regulation of cell morphogenesis for DE circRNAs (Fig. 4a and Additional file 8).
Gene ontology (GO) analysis of DE genes between D4 and D2
Between D4 and D2, significantly upregulated GO terms were mainly involved in: 1) cell cycle-related processes and DNA packaging for DE mRNAs, 2) cytokine-mediated signaling pathway, RNA secondary structure unwinding, and cellular response to type I interferon for DE lncRNAs, and 3) regulation of cell morphogenesis, microtubule bundle formation, and chromatin modification for DE circRNAs (Fig. 4b and Additional file 9). Additionally, significantly downregulated GO terms between D4 and D2 were mainly related to: 1) collagen fibril organization, extracellular matrix organization, and positive regulation of multicellular organismal process for DE mRNAs, 2) carboxylic acid metabolic process, single-organism catabolic process, and regulation of hormone levels for DE lncRNAs, and 3) chromatin modification, histone demethylation, and regulation of cell communication for DE circRNAs (Fig. 4b and Additional file 9).
Gene ontology (GO) analysis of DE genes between D8 and D4
Between D8 and D4, significantly upregulated GO terms were mainly involved in: 1) actin cytoskeleton organization-related processes and biological adhesion for DE mRNAs, 2) regulation of protein processing, negative regulation of cell proliferation, and system development for DE lncRNAs, and 3) regulation of posttranscriptional gene silencing, protein modification process, and protein localization to Golgi apparatus for DE circRNAs (Fig. 4c and Additional file 10). Additionally, significantly downregulated GO terms between D8 and D4 were mainly related to: 1) regulation of cell proliferation, system development, and single-multicellular organism process for DE mRNAs, 2) gene expression, RNA metabolic process, and macromolecule metabolic process for DE lncRNAs, and 3) positive regulation of gene expression, regulation of nucleobase-containing compound metabolic process, and regulation of nitrogen compound metabolic process for DE circRNAs (Fig. 4c and Additional file 10).
Gene ontology (GO) analysis of DE genes between D4 and D0
Between D4 and D0, significantly upregulated GO terms were mainly involved in: 1) lipid metabolism-related processes and monocarboxylic acid metabolic process for DE mRNAs, 2) regulation of hormone levels, protein processing, and cytokine-mediated signaling pathway for DE lncRNAs, and 3) chromatin modification, organelle organization, and cell part morphogenesis for DE circRNAs (Fig. 4d and Additional file 11). Furthmore, significantly downregulated GO terms were mainly enriched in: 1) cytoskeleton organization-related processes and multicellular organismal process for DE mRNAs, 2) negative regulation of growth, lipid catabolic process, and regulation of osteoclast differentiation for DE lncRNAs, and 3) negative regulation of macromolecule metabolic process, negative regulation of gene expression, and cellular protein modification process for DE circRNAs (Fig. 4d and Additional file 11).
Gene ontology (GO) analysis of DE genes between D8 and D0
Between D8 and D0, significantly upregulated GO terms were mainly involved in: 1) lipid metabolism-related processes and small molecule metabolic process for DE mRNAs and lncRNAs, and 2) chromatin modification, clathrin coat assembly, and chromosome organization for DE circRNAs (Fig. 4e and Additional file 12). Furthmore, significantly downregulated GO terms were mainly related to: 1) multicellular organismal process, system development, and cell adhesion for DE mRNAs, 2) lipid metabolic process, steroid metabolic process, and carbohydrate derivative metabolic process for DE lncRNAs, and 3) regulation of cellular process, regulation of nucleobase-containing compound metabolic process, and regulation of biological process for DE circRNAs (Fig. 4e and Additional file 12).
Function analysis of common DE genes during preadipocyte differentiation
Compared to the initial phase of differentiation (D0), 394 mRNAs, 37 lncRNAs, and 297 circRNAs were identified as common DE genes during the entire differentiation process (Fig. 5a-f and Additional file 13). GO analysis indicated that the terms were mainly enriched in: 1) cell adhesion, cell differentiation, and lipid metabolic process for DE mRNAs, 2) organic acid metabolic process, lipid metabolism-related processes, and single-organism process for DE lncRNAs, and 3) histone methylation, protein alkylation, and regulation of cellular response to growth factor stimulus for DE circRNAs (Fig. 5g and Additional file 14).
Construction of the lncRNA-mRNA co-expression network
To identify the key lncRNAs related to the regulation of lipid metabolic and cell differentiation processes, 132 common DE mRNAs associated with these two processes and 36 common DE lncRNAs targeting them were chosen to build the mRNA-lncRNA co-expression network. The results demonstrated that the co-expression network comprised 1328 connections and each lncRNA might correlate with multiple mRNAs (Fig. 6 and Additional file 15). More importantly, a total of 20 lncRNAs were found to be co-expressed with PPARG, a key adipocyte differentiation marker, while four lncRNAs including MSTRG.6375, SMIM4, MSTRG.88924, and MSTRG.65804 were shown to be co-expressed with other lipid metabolism-related markers such as APOE, LIPE, and ADIPOQ (Fig. 6 and Additional file 15), indicating that these lncRNAs might play an important role in regulating adipogenesis.
Validation of DE genes by qRT-PCR
A total of eighteen genes, including twelve mRNAs (CCND1, CCNA2, ARPC2, ARPC3, PPARG, CEBPA, PLIN1, SCD, ELOVL6, FABP4, IL11, and ATF3) related to cell cycle, actin cytoskeleton, cell differentiation, and lipid metabolism and six random lncRNAs (MSTRG.131380, MSTRG.146410, LOC100513133, MSTRG.28, MSTRG.42239, and MSTRG.62128), were selected for qRT-PCR verification. After comparison with the RNA-Seq data, similar expression trends for qRT-PCR were discovered, showing the strong consistency between qRT-PCR and RNA-Seq data (Fig. 7a-b).
As a key physiological process of normal body fat storage, preadipocyte differentiation provides a great opportunity for resolving the formation of fat deposition. In the present study, we observed that the shape of subcutaneous adipocytes changed from the shuttle-like form into the circlular form at the early stage of differentiation (D2), compared to the initial phase (D0) (Fig. 1a). Fat droplets were produced at D4, and then, clusters of large lipid droplets were formed at D8 (Fig. 1a). Consistent with the morphological changes of the adipocytes, triglyceride contents were observed to gradually increase, accompanying the increase in the differentiation time, strongly supporting the fact that preadipocyte differentiation is a complex process including both adipocyte growth and lipid deposits (Fig. 1b-c). Meanwhile, we identified more than 470 DE mRNAs, 58 DE lncRNAs, and 710 DE circRNAs among the different differentiation stages (Fig. 3d and Additional files 3, 4, 5, 6 and 7), which were involved in multiple biological processes including lipid metabolic and cell differentiation processes. As such, our data provide a comprehensive view of understanding the transcriptional regulation mechanism during the differentiation of porcine subcutaneous preadipocytes.
Generally speaking, preadipocyte differentiation mainly consists of three important stages including growth arrest, mitotic clonal expansion, and late events and terminal differentiation [25,26,27]. In the early stage of differentiation (D2), we found that upregulated mRNAs were mainly enriched in lipid metabolic process, while downregulated mRNAs were closely involved in cell cycle-related processes (Fig. 4a and Additional file 8). Many key markers, e.g., LIPE, APOE, PLIN1, DGAT2, ADIPOQ, and LPL for lipid metabolism and CDK1, CCND1, and E2F1 for cell cycle (Additional file 3), were identified in this study. Meanwhile, triglyceride contents at D2 were observed to increase but did not reach significant levels compared with those at D0, suggesting that upregulated expression levels of these lipid-related markers at the early stage of differentiation does not significantly alter triglyceride phenotype (Fig. 1b-c). However, further studies are needed to confirm this speculation. As expected, two critical well-documented markers involved in preadipocyte differentiation, PPARG and CEBPA, were shown to highly significant increase at D2 (Additional file 3). In the past few decades, PPARG and CEBPA have been deeply investigated for their determinant role in initiating and regulating preadipocyte differentiation [28, 29]. Besides, these two markers have been found to be involved in growth arrest of differentiation [30, 31]. Before preadipocyte differentiation, growth arrest is a necessary process for blocking the cell in the G1 phase [32, 33]; restricted cell proliferation has been found to appear at D2 of differentiation in human mesenchymal stem cells [32, 34]. In the G1/S phase, CDK1, CCND1, and E2F1 have been shown to play a key role in regulating the cell cycle . Furthermore, the expression of CDK1 and CCND1 decreased at the early stage of differentiation in human and porcine preadipocytes [34, 36]. Additionally, PPARG has been reported to control the cell cycle by decreasing the activity of E2F and CCND1 and upregulating the expression of cyclin-dependent kinase inhibitors [30, 37, 38], while CEBPA can repress the expression of E2F1, which results in the impairment of the ability of cell differentiation, and simultaneously suppresses cell proliferation . Consistent with these data, our results found that CDK1, CCND1, and E2F1 levels decreased significantly at D2, accompanying the increase in the levels of PPARG and CEBPA, demonstrating that PPARG and CEBPA might impact growth arrest at D2, via the downregulation of cell cycle-related markers.
Conversely, compared with the early stage of differentiation (D2), cell cycle-related marker levels were also found to increase at the middle stage of differentiation (D4) (Fig. 4b and Additional file 9). This is not surprising because preadipocytes reenter the cell cycle with at least one round of mitotic clonal expansion for increasing the proportion of adipocytes after the G1 phase growth-arrested; this is a synchronous process for adipogenesis . Moreover, mitotic clonal expansion is mainly related to the composition of fat inducers, especially insulin . Consistent with these results, our data supported that mitotic clonal expansion might be induced by insulin, which resulted in increased levels of cell cycle-related markers, e.g., CDK1, CCNA2, CCNE2, and E2F1 (Additional file 4). As mentioned above, CDK1 and E2F1 are key markers for controlling the G1/S phase of the cell cycle. Furthermore, CCNA2 is an important gene of the S phase of the cell cycle in combination with CDK1, while CCNE2 is a critical factor of the G1/S phase of the cell cycle in combination with E2F1 . Hence, the upregulation of these genes might promote cell cycle process from D2 to D4. Interestingly, our data indicated that extracellular matrix-related genes such as COL14A1 and MFAP5 were downregulated at D4, compared to D2 (Additional files 4 and 9). COL14A1 is a gene encoding fibril-associated collagen; it has been shown to have an antiproliferative role in reducing de novo DNA synthesis in 3T3-L1 preadipocytes . MFAP5 encodes a microfibril-associated glycoprotein and its expression levels were reduced during human preadipocyte differentiation . Here, the downregulation of these extracellular matrix-related genes at D4 might be associated with the changes of the adipocyte morphology, but additional investigation is required to decipher the role of the extracellular matrix in porcine adipogenesis.
From D4 to D8, lipid droplets grow larger and eventually form the mature adipocytes. Interestingly, the levels of many important markers related to actin cytoskeleton remodeling, e.g., ARPC2, ARPC3, and DSTN, were found to significantly increase at the later stage of differentiation (Additional files 5 and 10). As two major components of the actin cytoskeleton, ARPC2 and ARPC3 play an important role in adherens junction and intracellular motility of lipid vesicles [42, 43]. The knockdown of the ARP2/3 complex severely disrupted adipocyte differentiation . DSTN is another actin-depolymerizing factor, and its knockdown inhibited adipocyte differentiation of human stromal stem cells . Because the actin cytoskeleton is closely related to lipid droplet formation of the adipocytes , upregulation of these actin cytoskeleton-related genes at D8 might contribute to the formation of mature adipocytes at the later stage of differentiation. Meanwhile, we observed that the levels of other well-investigated markers associated with lipid lipolysis, e.g., ABHD5, LIPE, PNPLA2, and ACOX1, were significantly reduced at D8, compared to D4 (Additional file 5). Previously, ABHD5, LIPE, and PNPLA2 were confirmed to be the master regulators of tricylglycerol hydrolysis [47,48,49], while ACOX1 is the first enzyme involved in the fatty acid beta-oxidation process . Combined with the increase in triglyceride levels, these data support the claim that the changes of triglyceride phenotype might result from the downregulation of lipid lipolysis-related markers.
More importantly, compared to the initial phase of differentiation (D0), we discovered that the common DE mRNAs during the entire differentiation process were mainly involved in lipid metabolic and cell differentiation processes (Additional file 14). The levels of many well-known key markers, e.g., LIPE, PLIN1, DGAT2, PNPLA2, LPL, and SCD for lipid metabolic process, PPARG for cell differentiation, and APOE and ADIPOQ for both these processes, changed significantly during the entire differentiation process (Additional file 13), suggesting that these markers might play critical roles in phenotypic changes during the differentiation. In addition, the co-expressed network of the common DE lncRNAs and their target common DE mRNAs revealed that 36 lncRNAs targeted 132 lipogenesis-related mRNAs, indicating that these lncRNAs might participate in adipogenesis by positively or negatively regulating their target mRNAs (Fig. 6 and Additional file 15). For example, the co-expressed network showed that multiple lncRNAs could interact with PPARG, which is a decisive marker in adipogenesis (Fig. 6 and Additional file 15). Consequently, our results provide new evidence that lncRNAs are involved in porcine preadipocyte differentiation.
In summary, a genome-wide view of the expression profiling of mRNAs, lncRNAs, and circRNAs during porcine subcutaneous preadipocyte differentiation was investigated. Moreover, a large number of DE genes, which could contribute to the phenotypic changes in adipocytes at different stages of differentiation, were further identified. Our study provides a comprehensive basis for the expression levels of various RNAs during adipocyte differentiation, thus giving us a novel clue for understanding the mechanism of the molecular regulation of subcutaneous adipogenesis in pigs.
Analysis of variance
Burrows-Wheeler Aligner-maximal exact match
Coding potential assessing tool
Coding potential calculator
Dexamethasone, 3-isobutyl-1-methylxanthine, and insulin
Fetal bovine serum
False discovery rate
Fragments per kilobase of transcript per million reads mapped
Gene expression omnibus
Long non-coding RNA
Protein folding domain database
Quantitative real-time RT-PCR
RNA integrity number
Small nucleolar RNA
Small nuclear RNA
Spliced reads per billion mapping
Bosi P, Russo V. The production of the heavy pig for high quality processed products. Ital J Anim Sci. 2004;3:309–21.
Candek-Potokar M, Skrlep M. Factors in pig production that impact the quality of dry-cured ham: a review. Animal. 2012;6:327–38.
Ros-Freixedes R, Reixach J, Bosch L, Tor M, Estany J. Genetic correlations of intramuscular fat content and fatty acid composition among muscles and with subcutaneous fat in Duroc pigs. J Anim Sci. 2014;92:5417–25.
Gao J, Xu W, Wang J, Wang K, Li P. The role and molecular mechanism of non-coding RNAs in pathological cardiac remodeling. Int J Mol Sci. 2017;18:608.
Jiang Q, Wang J, Wu X, Ma R, Zhang T, Jin S, et al. LncRNA2Target: a database for differentially expressed genes after lncRNA knockdown or overexpression. Nucleic Acids Res. 2015;43:D193–6.
Kapranov P, Cheng J, Dike S, Nix DA, Duttagupta R, Willingham AT, et al. RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science. 2007;316:1484–8.
Chen J, Liu Y, Lu S, Yin L, Zong C, Cui S, et al. The role and possible mechanism of lncRNA U90926 in modulating 3T3-L1 preadipocyte differentiation. Int J Obes. 2016;41:299–308.
Liu W, Ma C, Yang B, Yin C, Zhang B, Xiao Y. LncRNA Gm15290 sponges miR-27b to promote PPARγ-induced fat deposition and contribute to body weight gain in mice. Biochem Biophys Res Commun. 2017;493:1168.
Nuermaimaiti N, Liu J, Liang X, Jiao Y, Zhang D, Liu L, et al. Effect of lncRNA HOXA11-AS1 on adipocyte differentiation in human adipose-derived stem cells. Biochem Biophys Res Commun. 2018;495:1878–84.
Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, et al. Circular RNA: a new star of noncoding RNAs. Cancer Lett. 2015;365:141–8.
Ebbesen KK, Kjems J, Hansen TB. Circular RNAs: identification, biogenesis and function. Biochim Biophys Acta. 2015;1859:163.
Legnini I, Di Timoteo G, Rossi F, Morlando M, Briganti F, Sthandier O, et al. Circ-ZNF609 is a circular RNA that can be translated and functions in myogenesis. Mol Cell. 2017;66:22–37.
Li H, Yang J, Wei X, Song C, Dong D, Huang Y, et al. CircFUT10 reduces proliferation and facilitates differentiation of myoblasts by sponging miR-133a. J Cell Physiol. 2017;233:4643–51.
Zambonelli P, Gaffo E, Zappaterra M, Bortoluzzi S, Davoli R. Transcriptional profiling of subcutaneous adipose tissue in Italian large white pigs divergent for backfat thickness. Anim Genet. 2016;47:306–23.
Davoli R, Gaffo E, Zappaterra M, Bortoluzzi S, Zambonelli P. Identification of differentially expressed small RNAs and prediction of target genes in Italian Large White pigs with divergent backfat deposition. Anim Genet. 2018; Epub ahead of print.
Jiang S, Wei H, Song T, Yang Y, Peng J, Jiang S. Transcriptome comparison between porcine subcutaneous and intramuscular stromal vascular cells during adipogenic differentiation. PLoS One. 2013;8:e77094.
Yi F, Yang F, Liu X, Chen H, Ji T, Jiang L, et al. RNA-seq identified a super-long intergenic transcript functioning in adipogenesis. RNA Biol. 2013;10:991–1001.
Wang J, Hua L, Chen J, Zhang J, Bai X, Gao B, et al. Identification and characterization of long non-coding RNAs in subcutaneous adipose tissue from castrated and intact full-sib pair Huainan male pigs. BMC Genomics. 2017;18:542.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.
Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:1303.3997. 2013. p. 1-3.
Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 2015;16:4.
Li Y, Zheng Q, Bao C, Li S, Guo W, Zhao J, et al. Circular RNA is enriched and stable in exosomes: a promising biomarker for cancer diagnosis. Cell Res. 2015;25:981–4.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26:136–8.
Morris JH, Vijay D, Federowicz S, Pico AR, Ferrin TE. CyAnimator: simple animations of Cytoscape networks. F1000Res. 2015;4:482.
Nakajima I, Muroya S, Chikuni K. Growth arrest by octanoate is required for porcine preadipocyte differentiation. Biochem Biophys Res Commun. 2003;309:702–8.
Tang QQ, Otto TC, Lane MD. Mitotic clonal expansion: a synchronous process required for adipogenesis. Proc Natl Acad Sci U S A. 2003;100:44–9.
Gregoire FM, Smas CM, Sul HS. Understanding adipocyte differentiation. Physiol Rev. 1998;78:783–809.
Tontonoz P, Hu E, Spiegelman BM. Stimulation of adipogenesis in fibroblasts by PPAR gamma 2, a lipid-activated transcription factor. Cell. 1994;79:1147–56.
Umek RM, Friedman AD, McKnight SL. CCAAT-enhancer binding protein: a component of a differentiation switch. Science. 1991;251:288–92.
Altiok S, Xu M, Spiegelman BM. PPARgamma induces cell cycle withdrawal: inhibition of E2F/DP DNA-binding activity via down-regulation of PP2A. Genes Dev. 1997;11:1987–98.
Porse BT, Pedersen TA, Xu X, Lindberg B, Wewer UM, Friis-Hansen L, et al. E2F repression by C/EBPalpha is required for adipogenesis and granulopoiesis in vivo. Cell. 2001;107:247–58.
Krawisz BR, Scott RE. Coupling of proadipocyte growth arrest and differentiation. I. Induction by heparinized medium containing human plasma. J Cell Biol. 1982;94:394–9.
Ailhaud G, Dani C, Amri EZ, Djian P, Vannier C, Doglio A, et al. Coupling growth arrest and adipocyte differentiation. Environ Health Perspect. 1989;80:17–23.
Marquez M, Alencastro F, Madrigal A, Jiminez J, Blanco G, Gureghian A, et al. The role of cellular proliferation in adipogenic differentiation of human mesenchymal stem cells. Stem Cells Dev. 2017;26:1578–95.
Lim S, Kaldis P. Cdks, cyclins and CKIs: roles beyond cell cycle regulation. Development. 2013;140:3079–93.
Kim HS, Hausman GJ, Hausman DB, Martin RJ, Dean RG. The expression of cyclin D1 during adipogenesis in pig primary stromal-vascular cultures. Obes Res. 2001;9:572–8.
Wang C, Fu M, D'Amico M, Albanese C, Zhou JN, Brownlee M, et al. Inhibition of cellular proliferation through IkappaB kinase-independent and peroxisome proliferator-activated receptor gamma-dependent repression of cyclin D1. Mol Cell Biol. 2001;21:3057–70.
Morrison RF, Farmer SR. Role of PPARgamma in regulating a cascade expression of cyclin-dependent kinase inhibitors, p18(INK4c) and p21(Waf1/Cip1), during adipogenesis. J Biol Chem. 1999;274:17088–97.
Qiu Z, Wei Y, Chen N, Jiang M, Wu J, Liao K. DNA synthesis and mitotic clonal expansion is not a required step for 3T3-L1 preadipocyte differentiation into adipocytes. J Biol Chem. 2001;276:11988–95.
Ruehl M, Erben U, Schuppan D, Wagner C, Zeller A, Freise C, et al. The elongated first fibronectin type III domain of collagen XIV is an inducer of quiescence and differentiation in fibroblasts and preadipocytes. J Biol Chem. 2005;280:38537–43.
Vaittinen M, Kolehmainen M, Rydén M, Eskelinen M, Wabitsch M, Pihlajamäki J, et al. MFAP5 is related to obesity-associated adipose tissue and extracellular matrix remodeling and inflammation. Obesity (Silver Spring). 2015;23:1371–8.
Verma S, Shewan AM, Scott JA, Helwani FM, den Elzen NR, Miki H, et al. Arp2/3 activity is necessary for efficient formation of E-cadherin adhesive contacts. J Biol Chem. 2004;279:34062–70.
Zhou K, Sumigray KD, Lechler T. The Arp2/3 complex has essential roles in vesicle trafficking and transcytosis in the mammalian small intestine. Mol Biol Cell. 2015;26:1995–2004.
Yang W, Thein S, Lim CY, Ericksen RE, Sugii S, Xu F, et al. Arp2/3 complex regulates adipogenesis by controlling cortical actin remodelling. Biochem J. 2014;464:179–92.
Chen L, Hu H, Qiu W, Shi K, Kassem M. Actin depolymerization enhances adipogenic differentiation in human stromal stem cells. Stem Cell Res. 2018;29:76–83.
Padilla-Benavides T, Velez-Delvalle C, Marsch-Moreno M, Castro-Muñozledo F, Kuri-Harcuch W. Lipogenic enzymes complexes and cytoplasmic lipid droplet formation during adipogenesis. J Cell Biochem. 2016;117:2315–26.
Haemmerle G, Zimmermann R, Zechner R. Letting lipids go: hormone-sensitive lipase. Curr Opin Lipidol. 2003;14:289–97.
Lass A, Zimmermann R, Haemmerle G, Riederer M, Schoiswohl G, Schweiger M, et al. Adipose triglyceride lipase-mediated lipolysis of cellular fat stores is activated by CGI-58 and defective in Chanarin-Dorfman syndrome. Cell Metab. 2006;3:309–19.
Smirnova E, Goldberg EB, Makarova KS, Lin L, Brown WJ, Jackson CL. ATGL has a key role in lipid droplet/adiposome degradation in mammalian cells. EMBO Rep. 2006;7:106–13.
Varanasi U, Chu R, Chu S, Espinosa R, LeBeau MM, Reddy JK. Isolation of the human peroxisomal acyl-CoA oxidase gene: organization, promoter analysis, and chromosomal localization. Proc Natl Acad Sci U S A. 1994;91:3107–11.
This work was supported by the National Transgenic Project of China (2018ZX0800928B), the National Key Technology Research and Development Program of the Ministry of Science and Technology of China (2015BAD03B01), and the National Natural Science Foundation of China (31501930).
Availability of data and materials
The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.
All animal experiments in this study were approved by the Animal Ethics Committee at Nanjing Agricultural University (Nanjing, China).
Consent for publication
The authors declare that they have no competing interests.
Primers of quantitative real-time RT-PCR. (XLSX 71 kb)
List of expressed genes during porcine preadipocyte differentiation. Gene ID, FPKM or SRPBM value, gene name, biotype, position and description of all acquired genes with FPKM or SRPBM ≥1 in more than 50% of samples during porcine preadipocyte differentiation. (XLSX 2094 kb)
List of differentially expressed (DE) genes for D2 vs D0. Gene ID, fold change, gene name, biotype, and position of differentially expressed mRNAs, lncRNAs, and circRNAs between D2 versus D0. (XLSX 345 kb)
List of differentially expressed (DE) genes for D4 vs D2. Gene ID, fold change, gene name, biotype, and position of differentially expressed mRNAs, lncRNAs, and circRNAs between D4 versus D2. (XLSX 167 kb)
List of differentially expressed (DE) genes for D8 vs D4. Gene ID, fold change, gene name, biotype, and position of differentially expressed mRNAs, lncRNAs, and circRNAs between D8 versus D4. (XLSX 295 kb)
List of differentially expressed (DE) genes for D4 vs D0. Gene ID, fold change, gene name, biotype, and position of differentially expressed mRNAs, lncRNAs, and circRNAs between D4 versus D0. (XLSX 385 kb)
List of differentially expressed (DE) genes for D8 vs D0. Gene ID, fold change, gene name, biotype, and position of differentially expressed mRNAs, lncRNAs, and circRNAs between D8 versus D0. (XLSX 206 kb)
GO terms for upregulated and downregulated mRNAs, lncRNAs, and circRNAs from D2 vs D0. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 4497 kb)
GO terms for upregulated and downregulated mRNAs, lncRNAs, and circRNAs from D4 vs D2. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 3928 kb)
GO terms for upregulated and downregulated mRNAs, lncRNAs, and circRNAs from D8 vs D4. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 4067 kb)
GO terms for upregulated and downregulated mRNAs, lncRNAs, and circRNAs from D4 vs D0. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 4373 kb)
GO terms for upregulated and downregulated mRNAs, lncRNAs, and circRNAs from D8 vs D0. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 3796 kb)
List of common differentially expressed (DE) genes for three differentiation stages (D2, D4, and D8) vs D0. Gene ID, fold change, gene name, biotype, and position of common differentially expressed (DE) mRNAs, lncRNAs, and circRNAs. (XLSX 5248 kb)
GO terms for common differentially expressed (DE) mRNAs, lncRNAs, and circRNAs from three differentiation stages (D2, D4, and D8) versus D0. Terms marked gray background represent significant enrichment, while other terms represent non-significant enrichment. (XLSX 1763 kb)
The relationship of common DE lncRNAs and their targeted common DE mRNAs involved in lipid metabolic and cell differentiation processes. (XLSX 34 kb)