INTRODUCTION
Adipose tissue in pigs is associated with important traits of carcass and meat quality [
1], and backfat deposition greatly influences porcine growth performance, meat production and final farming profit [
2]. Importantly, the pig is emerging as an attractive biomedical model for studying obesity and related diseases in human because of the similar physiology, anatomy, and metabolic features [
3,
4]. Therefore, understanding the molecular mechanism of backfat deposition in pigs is not only conducive to the progress of genetic breeding, but also expands the comprehension of obesity related metabolic diseases in human.
microRNAs (miRNAs) are a kind of small non-coding RNAs of about 19 to 23 nucleotides in length that negatively regulates gene expression [
5]. lncRNAs, a class of non-protein coding RNAs of >200 nucleotides in length, are poorly conserved and expressed in a cell-, tissue-, and stage-specific manner, and lncRNAs can be classified into sense, antisense, intergenic, and bidirectional lncRNAs according to their genomic position [
6]. Accumulating evidence indicates that miRNAs play various crucial roles in the processes of adipocyte differentiation, adipogenesis and lipid metabolism [
7]. In recent years, studies on identification and functional analysis of porcine lncRNAs were progressively performed, and several differentially expressed lncRNAs and potential regulatory lncRNAs were obtained either in adipose tissue or in the process of adipogenesis from two different pig breeds [
8–
13].
Daweizi (DWZ), a famous indigenous pig breed in China, is characterized by tender meat, slow growth rate and high fat percentage [
14]. In contrast, Yorkshire, a worldwide well-known pig breed, exhibits rapid growth rate with low fat percentage under the intensive selection. In view of the distinctive difference in term of fat content, the two pig breeds can be regarded as the appropriate objects to study the molecular mechanism of fat deposition.
However, the molecular mechanism of the obvious difference in backfat deposition between DWZ and Yorkshire pigs has not yet been studied. In the present investigation, the expression profiles of miRNAs, lncRNAs, and mRNAs were compared in backfat tissue between DWZ and Yorkshire pigs by RNA sequencing (RNA-seq) technology. Then, the differentially expressed miRNAs (DE miRNAs), differentially expressed lncRNAs (DE lncRNAs) and differentially expressed mRNAs (DE mRNAs) associated with porcine adipogenesis were identified, and the functional enrichments and lncRNA-miRNA-mRNA interaction network were further analyzed. The present investigation provides more insights into the mechanism of backfat deposition from the point of view of miRNAs, lncRNAs, and mRNAs.
MATERIALS AND METHODS
Experimental animal and tissue collection
Three healthy castrated male DWZ pigs (180-day-old, average slaughter weight 73.83±1.88 kg) and three healthy castrated male Yorkshire pigs (180-day-old, average slaughter weight 121.30±1.33 kg) were raised under the same conditions at Tianfu Ecological Agricultural Limited Company in Hunan province, China. After slaughtered in a commercial slaughter plant, the left side carcass was hung upside down, and the midline backfat thickness at the position of the thickest point in shoulder, last rib and lumbosacral junction was separately measured, and then the average value was calculated. Subsequently, backfat tissues collected at the position of the thickest point in shoulder were divided into two parts. One part was fixed in 4% paraformaldehyde for histological analysis, and the other part was immediately frozen in liquid nitrogen and stored in a −80°C refrigerator. All the experiments in this study were reviewed and approved by the Institutional Animal Care and Use Committee of Hunan Institute of Animal and Veterinary Science (Approval number: 20200110).
Histological analysis of backfat tissue
The paraformaldehyde-fixed backfat samples were embedded in paraffin. The serial tissue sections were cut using cryostat (Leica RM2235; Leica company, Wetzlar, Germany) and were then stained with hematoxylin/eosin. The sections were viewed at 200× magnification using microscope (Leica DM3000, Germany), and five areas were randomly selected in each sample for measuring the diameter of adipocyte.
RNA extraction and library construction
Total RNA from backfat tissue was isolated using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instruction. The quantity and integrity of total RNA were assessed by Agilent 2100 bioanalyzer (Santa Clara, CA, USA), which exhibited that the RNA integrity number was more than 8.0 for each sample. The main processes of small RNA library construction include: i) Total RNA was purified by electrophoretic separation on a 15% urea denaturing polyacrylamide gel electrophoresis (PAGE) gel, and small RNA regions corresponding to the 18 to 30 nt bands were excised and recovered. ii) Small RNAs were ligated to adenylated 3′-adapters annealed to unique molecular identifiers (UMI), followed by ligation of 5′-adapters. iii) The adapter-ligated small RNAs were transcribed into cDNA by SuperScript II Reverse Transcriptase (Invitrogen, USA) and then several rounds of polymerase chain reaction (PCR) amplification were performed to enrich cDNA fragments. iv) PCR products with 110 to 130 bp in length were acquired on PAGE gel and were then purified by QIAquick Gel Extraction Kit (QIAGEN, Düsseldorf, Germany). Finally, these PCR products were circularized and then sequenced using the DNBseq platform (BGI-Shenzhen, China). The brief procedures for construction of mRNAs and lncRNAs libraries were as follows: i) Ribosomal RNA (rRNA) was removed using MGIEasy rRNA kit (BGI, China), and the remained RNA was fragmented into small pieces using divalent cations under elevated temperature. ii) The cleaved RNA fragments were transcribed into first strand cDNA using reverse transcriptase and random primers, followed by second strand cDNA synthesis using DNA polymerase I and RNase H with dUTP instead of dTTP. iii) These cDNA fragments were ligated with an ‘A’ base and sequencing adapter. iv) The products were enriched with PCR and then heat denatured and circularized by the splint oligo sequence. v) The single strand circle DNA was formatted as the final library which was sequenced on DNBseq platform (BGI-Shenzhen, China).
Sequencing data analysis
For analysis of small RNA sequencing data, the raw reads containing polyA, shorter than 18 nt in length, with low-quality reads (the number of bases with quality score less than 10 is ≤4 and with quality score less than 13 is ≤6), with 5′-adapter contamination, without 3′-adapter or inserted fragments were filtered to obtain clean reads by SOAPnuke (v1.5.0). Subsequently, the clean reads were mapped to porcine reference genome (GCF_000003025.6_Sscrofa11.1) and other small RNA databases including miRNA from miRbase (v22.1), and siRNA, piRNA, and snoRNA from NCBI GenBank by Bowtie2. Particularly, cmsearch was performed for mapping Rfam database (
http://rfam.xfam.org/). To make each unique sRNA have a unique annotation, the priority order of miRbase>pirnabank>snoRNA>Rfam>other sRNA was followed to traverse the annotation. The unannotated sequences were used to predict novel miRNA candidates by miRDeep2 program (v0.1.3) based on the secondary structure. The novel miRNAs were aligned to mature miRNAs from other mammals in the miRbase by BLAST.
For analysis of mRNA and lncRNA sequencing data, clean reads were obtained from raw reads by removing adapter pollution, low-quality reads (more than 20% bases with quality score less than 15) and reads whose unknown base ratio was greater than 5% using SOAPnuke. These clean reads were aligned to the porcine reference genome using HISAT2 (v2.0.4), and Bowtie2 (v2.2.5) was applied to align the clean reads to the reference sequence. Then the mapped reads were assembled using StringTie. The transcripts were then screened by Pfam database and three softwares including coding potential calculator (score <0), txCdsPredict (score <500), and coding-non-coding index (score <0). The transcripts that unmatched Pfam database and passed through at least two of the three softwares were considered to be lncRNAs.
Identification of differentially expressed transcripts
After mapping clean reads to the porcine reference sequence by Bowtie2, the expressions of small RNAs were calculated by counting absolute numbers of molecules using UMI, and then DE miRNAs analysis was performed using DEGseq. The fragments per kilobase of transcript per million (FPKM) reads mapped value was used to estimate the expressions of lncRNAs and mRNAs by RSEM software (v1.2.12), and analysis of DE mRNAs and DE lncRNAs were examined by DEseq2. Q value (adjusted p value) ≤0.05 and |log2(DWZ/Yorkshire)| ≥1 were set as thresholds for significant differential expression. The heatmap of DE transcripts was drawn by pheatmap (v1.0.8) with default parameter.
Target genes prediction and functional analysis of DE miRNAs and DE lncRNAs
The predicted target genes of miRNAs were implemented by using RNAhybrid, miRanda and TargetScan softwares. lncRNAs can regulate target genes by acting
in cis and
in trans. If lncRNAs and genes exhibit similar expression patterns, their biological functions may be highly correlated. Accordingly, the targets genes of DE lncRNAs were predicted as follows. Based on the Spearman and Pearson correlation coefficients of lncRNA-mRNA pair being ≥0.6, mRNAs located within 50 kb upstream and 50 kb downstream of DE lncRNAs were selected for
cis-acting regulation, and RNAplex was utilized to predict the combination of lncRNA and mRNA for
trans-acting regulation with binding energy <–20. To explore the potential biological functions of DE transcripts, GO and KEGG analysis were carried out for the DE mRNAs and targets of DE miRNAs and DE lncRNAs based on GO (
http://www.geneontology.org/ ) and KEGG (
https://www.kegg.jp/ ) databases, and then the functional enrichment analysis was performed by Phyper based on Hypergeometric test. The significant levels of terms and pathways were calculated with a rigorous threshold (Q value ≤0.05) by Bonferroni.
Construction of PPI and lncRNA-miRNA-mRNA network
The STRING was applied to construct protein-protein interaction (PPI) network for DE mRNAs, which was visualized by Cytoscape software. As mentioned above, the target genes of DE miRNAs and DE lncRNAs were predicted, and the interactions of lncRNA-miRNA, lncRNA-mRNA, and miRNA-mRNA were obtained. Consequently, the potential regulatory network of lncRNA-miRNA-mRNA was also visualized using Cytoscape.
Quantitative real-time polymerase chain reaction analysis
The validation of miRNA expression levels was detected by stem-loop QPCR method [
15]. cDNA was synthesized using RevertAid first strand cDNA synthesis kit (K1622, Fermentas) according to the manufacturer’s instructions. QPCR analysis was performed using SYBR Green Supermix (Biomed, Beijing, China) on CFX96 machine (Bio-Rad, Hercules, CA, USA). Porcine glyceraldehyde-3-phosphate dehydrogenase (
GAPDH) was used as endogenous control for lncRNAs and mRNAs, and U6 snRNA was used for miRNAs. Each QPCR reaction was performed in triplicate, and the relative expression level of transcripts was calculated using the 2
−ΔΔCt method. The sequences of QPCR primers for mRNAs, miRNAs and lncRNAs are listed in
Supplementary Table S1.
Statistical analysis
Statistical analysis of the data from adipocyte diameter and QPCR assay was performed by one-way analysis of variance program with SPSS 20.0 software. Mean values and standard error were presented, and p value of less than 0.05 was deemed statistically significant (* p≤0.05 and ** p≤0.01).
DISCUSSION
Fat deposition is closely linked with the genetic background and is influenced by various transcription factors, crucial genes and signaling pathways. To better understand the regulatory network controlling backfat deposition, expression profiles of miRNAs, lncRNAs, and mRNAs in backfat tissue between Chinese indigenous (obese-type) DWZ and foreign lean-type Yorkshire pig breeds were analyzed for the first time by RNA-seq technology. A total of 1,625 DE transcripts were identified, including 72 DE miRNAs, 183 DE lncRNAs and 1,370 DE mRNAs. And the expression levels of several DE transcripts were validated by QPCR analysis, indicating that data of RNA-seq was reliable with strong consistency.
Several of highly abundant miRNAs identified in this study are closely related to adipogenesis. miR-26a, miR-206, miR-378, miR-499, and miR-191 have been reported to participate in adipocyte differentiation and adipogenesis [
17–
21]. Meanwhile, some of the DE miRNAs are indispensable for adipogenesis, lipid metabolism and adipose development. For instances, miR-375 exerts a positive regulatory effect on adipocyte differentiation. miR 204-5p favors the adipogenic differentiation and lipid synthesis, of human adipose-derived mesenchymal stem cells by suppressing Wnt/β-catenin signaling, of 3T3–L1 cells by targeting KLF transcription factor 3 (
KLF3). Recently, transcriptome and miRNAome analysis have been used to identify miRNA expression profiles and characterize the possible regulatory relationships for backfat deposition in pigs [
1,
22–
26]. In these researches and our study, miR-26a, miR-206, miR-99a-5p, miR-16, let-7f-5p, miR-1, and miR-24-3p were the top 20 most abundant miRNAs, additionally, miR-375, miR-204–5p, miR-205, miR-192, miR-194a-5p, miR-29b, miR-29c, miR-708–5p, miR-127, miR-369, miR-493–5p, miR-323, miR-432–5p, and miR-493–5p were the identified DE miRNAs. However, some of the DE miRNAs in our study were different from that in previous research, which may be caused by the differences in pig breeds, days of age, or threshold of the fold change.
To identify the differential regulatory networks in backfat tissue between the two pig breeds, GO and KEGG enrichment analysis were fulfilled. For DE mRNAs, several signaling pathways associated with lipid metabolism and adipogenesis were discovered, such as PPAR, p53, AMPK, and PI3K-Akt signaling pathways [
8,
27–
29]. Several hub proteins were identified in PPI network of DE mRNAs.
ACOX1,
ALDOB, protein kinase AMP-activated catalytic subunit alpha 2 (
PRKAA2),
JAK3, and protein tyrosine phosphatase non-receptor type 11 (
PTPN11) have been confirmed to be closely related to β-oxidation of fatty acid, lipogenesis, cholesterol metabolism, and lipodystrophy [
30–
35]. Moreover, several genes, including
PPARδ,
PPARα, fatty acid binding protein 3 (
FABP3), acyl-CoA synthetase long chain family member 4 (
ACSL4), B cell lymphoma 2 (
BCL2),
FOXO3 and TBC1 domain family member 1 (
TBC1D1), are closely associated with fatty acid metabolism, adipogenic differentiation and cell cycle progression.
PPARδ and
PPARα are vital regulators in cell differentiation, lipid accumulation, fatty acid oxidation and lipid catabolism.
FABP3 is positively correlated with the backfat thickness in beef cattle, and it regulates transport of fatty acids and lipid deposition. Furthermore,
FABP3 is notably downregulated in subcutaneous adipose tissue from Yorkshire pigs than that in Chinese indigenous Laiwu pigs [
8], which is similar to our study.
ACSL4 converts free long-chain fatty acids into fatty acyl-CoA esters which are the key intermediates in the synthesis of complex lipids. Previous study demonstrated that the expression of
ACSL4 gradually increases during adipogenesis of porcine primary intramuscular preadipocytes and
ACSL4 knockdown decreases lipid accumulation. However, another research found that the expression of
ACSL4 was gradually decreased during different differentiation stages of subcutaneous preadipocytes in Erhualian pigs [
10], and our result also exhibited that the expression of
ACSL4 in backfat tissue is significantly lower in DWZ pigs than that in Yorkshire pigs. In addition,
AKT2 knockdown remarkably reduces preadipocyte proliferation, adipogenic differentiation and fat mass in human, but the present result showed that
AKT2 expression was drastically lower in DWZ pigs than that in Yorkshire, which is contrary to the result reported in a former article about Bamei and Large White pigs. The discrepancy may be because of the different experimental methods and cell sources.
BCL2 mediates obedience to proliferation and resistance to apoptosis in adipocytes. In addition,
FOXO3 and
TBC1D1 have been shown to be involved in lipid accumulation, fatty acid oxidation and obesity development. In our study, lower expressions of
ACOX1,
AKT2,
PRKAA2,
ACSL4,
FOXO3 and higher expressions of
ALDOB,
JAK3,
PTPN11,
PPARδ,
PPARα,
FABP3,
BCL2,
TBC1D1 were observed in DWZ pigs. Taken together, these regulatory relationships may partly illuminate the mechanism of porcine backfat deposition.
The biological function of lncRNAs is generally mediated by their targets. Therefore, the targets of DE lncRNAs were predicted and then underwent enrichment analysis. Several pathways were enriched, among which AMPK, apelin and insulin signaling pathways are closely related to adipocyte differentiation and lipid metabolism. AMPK signaling pathway is believed to act as a key master switch that modulates cholesterol synthesis and lipid metabolism [
36,
37]. The apelin signaling pathway inhibits adipogenesis and lipolysis through distinct molecular pathways. Insulin is a key regulator and activates the transcription and proteolytic maturation of sterol regulatory element binding transcription factor 1c (SREBP1c), and then SREBP1c induces the expression of a family of genes involves in fatty acid synthesis, such as acetyl-CoA carboxylase (
ACC) and fatty acid synthase (
FAS).
The present study indicated that BGIR9823_87741, BGIR 9823_78329, BGIR9823_79919, and BGIR9823_85300 targets
FABP3,
BCL2,
FOXO3, and
PPARα, respectively. Meanwhile, it is noteworthy that XR_002341917.1, BGIR9823_80765, BGIR9823_79960, BGIR9823_93430, BGIR9823_84443, BGIR9823_80428, BGIR9823_85841 and XR_002340248.1 potentially regulates H2A histone family member Y (
H2AFY), neuronal regeneration related protein (
NREP), carboxypeptidase A4 (
CPA4), optic atrophy 1 (
OPA1), FKBP prolyl isomerase 5 (
FKBP5), monoacylglycerol O-acyltransferase 1 (
MGAT1), angiopoietin like 2 (
ANGPTL2) and protein kinase AMP-activated non-catalytic subunit gamma 2 (
PRKAG2), respectively. In our study,
H2AFY,
NREP, and
CPA4 are highly expressed in Yorkshire pig. Former research showed that deletion of
H2AFY results in lipid accumulation in murine liver [
38]. Treatment of HepG2 cells with
NREP knockdown exhibits greater lipid droplet accumulation and increases triglyceride and cholesterol content through TGFβR/PI3K/AKT signaling pathway.
CPA4 knockdown enhances differentiation of human preadipocytes and
CPA4 expression in subcutaneous adipose tissue negatively correlates with indices of insulin sensitivity. These results demonstrated that
H2AFY,
NREP, and
CPA4 might inhibit fat deposition in Yorkshire pigs. Moreover,
FABP5 and
MGAT1 were markedly upregulated in backfat tissue in DWZ pigs.
FKBP5 expression in human subcutaneous adipose tissue tends to be increased in type II diabetes subjects and is associated with genes involved in lipid metabolism and adipogenesis. Another research suggested that mice lacking
FKBP5 gene had reduce body weight and were resistant to diet-induced obesity, and knockdown of
FKBP5 in 3T3-L1 cells had a strong anti-adipogenic impact. It has been found that
MGAT1 encodes the enzyme that catalyzes monoacylglycerol and fatty acyl-CoA to form diacylglycerol, which is indispensable for triacylglycerol synthesis, and
MGAT1 knockdown exhibits a significant reduction in lipid accumulation.
ANGPTL2, predominantly secreted by adipose tissue, enhances fatty acid synthesis and lipid accumulation in mice, and
ANGPTL2 knockdown inhibits adipogenic differentiation of 3T3-L1 cells. On above basis,
FABP5,
MGAT1, and
ANGPTL2 might enhance adipogenic differentiation and lipid formation in DWZ pigs. Accordingly, the DE lncRNAs target these genes might play crucial roles in adipogenesis and lipid metabolism in porcine backfat tissue, and further studies are needed to be fulfilled to validate this speculation.
The target relationships between miRNAs and mRNAs were further analyzed.
ANXA3, the predicted target of miR-122-3p which plays important roles in cholesterol synthesis and lipogenesis [
39,
40], is downregulated at an early phase of adipocyte differentiation in 3T3–L1 cells, and suppression of
ANXA3 causes elevation of the
PPARγ2 mRNA level and lipid droplet accumulation. In addition, miR-583-3p, miR-161-3p, miR-671-3p, and miR-817-3p targets activating transcription factor 3 (
ATF3), triacylglycerol synthase 1 (
TGS1), protein kinase C and casein kinase substrate in neurons 2 (
PACSIN2), and ataxin 2 (
ATXN2), respectively, and these genes have been verified to be associated with fat weight, triglyceride synthase and lipid droplet formation. Furthermore,
AQP7 is the predicted target of miR-194a-5p, and previous report showed that the body weight and fat mass increases significantly in
AQP7 knockout mice, and adipocytes are large and exhibits accumulation of triglyceride by elevating adipose glycerol kinase activity [
41]. These results indicated that the BGIR9823_87926/miR-194a-5p/
AQP7 network may affect the process of adipocyte differentiation or adipogenesis (
Figure 8), but their functions and interactions needed further validation.