Integrated Analysis of Liver Transcriptome, miRNA, and Proteome of Chinese Indigenous Breed Ningxiang Pig in Three Developmental Stages Uncovers Significant miRNA–mRNA–Protein Networks in Lipid Metabolism

Liver is an important metabolic organ of mammals. During each transitional period of life, liver metabolism is programmed by a complex molecular regulatory system for multiple physiological functions, many pathways of which are regulated by hormones and cytokines, nuclear receptors, and transcription factors. To gain a comprehensive and unbiased molecular understanding of liver growth and development in Ningxiang pigs, we analyzed the mRNA, microRNA (miRNA), and proteomes of the livers of Ningxiang pigs during lactation, nursery, and fattening periods. A total of 22,411 genes (19,653 known mRNAs and 2758 novel mRNAs), 1122 miRNAs (384 known miRNAs and 738 novel miRNAs), and 1123 unique proteins with medium and high abundance were identified by high-throughput sequencing and mass spectrometry. We show that the differences in transcriptional, post-transcriptional, or protein levels were readily identified by comparing different time periods, providing evidence that functional changes that may occur during liver development are widespread. In addition, we found many overlapping differentially expressed genes (DEGs)/differentially expressed miRNAs (DEMs)/differentially expressed proteins (DEPs) related to glycolipid metabolism in any group comparison. These overlapping DEGs/DEMs/DGPs may play an important role in functional transformation during liver development. Short Time-series Expression Miner (STEM) analysis revealed multiple expression patterns of mRNA, miRNA, and protein in the liver. Furthermore, several diverse key Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, including immune defense, glycolipid metabolism, protein transport and uptake, and cell proliferation and development, were identified by combined analysis of DEGs and DGPs. A number of predicted miRNA–mRNA–protein pairs were found and validated by qRT-PCR and parallel reaction monitoring (PRM) assays. The results provide new and important information about the genetic breeding of Ningxiang pigs, which represents a foundation for further understanding the molecular regulatory mechanisms of dynamic development of liver tissue, functional transformation, and lipid metabolism.

Introduction

Metabolism of lipid in animals is highly dependent upon liver development and physiology (Liu et al., 2018). The type and content of lipids determine the quality of meat. Data from several mammalian studies suggest that microRNAs (miRNAs) play an important regulatory role in liver lipid metabolism. For example, miR-122, which is highly expressed in mammalian liver, specifically inhibitions the expression of genes involved in cholesterol biosynthesis by targeting the regulation of cationic amino acid transporter 1 (CPT1α), thereby regulating liver fatty acid (FA) and cholesterol synthesis (Iliopoulos et al., 2010). One miRNA can regulate the expression of several mRNAs, and the expression of one mRNA can be simultaneously regulated by several miRNAs. Li H. et al. (2016) found that fatty acid desaturase 1 (FADS1) was targeted by miR-365-3p, miR-218-5p, miR-181a-5p, miR-181b-5p, miR-29a-3p, and miR-23b-3p, while fatty acid desaturase 2 (FADS2) was targeted by miR-30c-1-3p, which was found in the study of differences in miRNA expression in liver of 20- and 30-week-old hens. In addition, it has been reported that fatty acid synthase (FASN) and stearoyl-COA desaturease (SCD) can be jointly targeted directly by miR-212-5p (Guo et al., 2017) and miR-27a (Zhang M. Y. et al., 2017), affecting the accumulation of triglyceride (TG) in hepatocytes.

In vertebrates, the transformation of liver metabolic program during each transitional period is a complex molecular regulatory system, whether in embryo, postnatal developmental stages and adulthood. For example, the yolk lipids of newly hatched chicks are rapidly depleted and metabolically converted to a carbohydrate-based energy source. Due to a poorly developed digestive system, they cannot process feed as efficiently as adult chickens (Uni et al., 2005). By using delayed feeding for 48 h to prevent liver metabolic transformation, found that delayed feeding resulted in increased miR-20b expression, while its target FADS1 expression was conversely decreased. It shows that the FADS1 gene is mediated by miRNA and participates in this transformation before and after transcription (Hicks et al., 2019). At the same time, the expression of miR-33 was significantly reduced, and miR-33, a well-known lipid metabolism regulator in vertebrates, has been proved to inhibit the oxidation of FAs in the liver (Shao et al., 2019). It is suggested that de novo synthesis of FAs can be almost inhibited during the hepatic metabolic transition, which may be caused by the decreased expression of SREBF1 and its downstream FADS, and miRNAs may be the key regulatory factors of the hepatic metabolic pathway during this metabolic transition. miRNAs may also play an important role in the prolongation of circulation responses of liver FAs. For example, Zhang M. et al. (2017) demonstrated that miR-218-5p can directly target ELOVL5 and indirectly promote the synthesis of LCPUFA in the liver, and this regulation may be dependent on miRNAs. Subsequently, it has been found that ELOVL2 and ultra-long chain PUFA concentrations were downregulated in the dedifferentiated primary human hepatocytes (PHHs). Furthermore, a variety of miRNAs were differentially expressed during the dedifferentiation of PHHs, such as up-regulated miR-27a and down-regulated miR-30, suggesting that these miRNAs are involved in the synthesis, accumulation, and secretion of PHH lipids (Kiamehr et al., 2019).

Ningxiang pig is one of the four most famous native pig breeds in China, with a domestication history of more than 5000 years. Ningxiang pigs, known for their meaty taste, have more than 5% intramuscular fat content, while compared with about 2% for imported commercial pigs (He et al., 2020). However, the molecular mechanism of Ningxiang pig liver regulating FA synthesis and metabolism and other meat quality related traits is poorly understood. FA synthesis and metabolism are affected by many factors such as feed, variety, and environment. In addition, there are many genes, transcripts, and proteins involved in each stage of tissue development of different pig breeds. The single factor molecular research based on the comparison of different pig breeds of the same age or weight cannot fully elucidate the molecular mechanisms and interactions in the organism. Therefore, the study of the dynamic expression of tissue development and the changes of transcription factors are important bases for understanding the transformation of tissue function during the development of pigs. However, due to the limitation of transcriptional regulation, transcriptional changes cannot accurately describe various changes in protein levels (Li et al., 2018). In this case, proteomics can help reveal changes in metabolic pathways in tissues, not only for biomarkers in breeding programs, but also for defining animal models (Bassols et al., 2014; Cui et al., 2016). Label-free liquid chromatography–mass spectrometry (LC–MS) has been increasingly used to quantify protein expression and compare different samples compared to traditional proteomic profiling methods (two-dimensional polyacrylamide gel electrophoresis and mass spectrometry). This approach may also offer the possibility of better studying low molecular weight or high molecular weight hydrophobic proteins, which are very challenging for analysis using traditional proteomic methods (Neilson et al., 2011). Bovo et al. (2018) used label-free LC–MS proteomics to compare the liver proteomic characteristics of Italian Duroc and Italian Great White pigs, and analyzed a total of 25 proteins that clearly distinguish Italian Duroc and Italian Great White pigs.

It has been reported that the scalability of protein studies techniques are not as good as nucleic acid (Manzoni et al., 2018). For this purpose, we combined analysis of mRNA-seq, miRNA-seq, and proteomics in liver from Ningxiang pigs at different ages (30, 90, and 210 days after birth), which is the first report on the integrated analysis of Ningxiang pig liver, also validated by qRT-PCR and parallel reaction monitoring (PRM) assays (Peterson et al., 2012). Our data suggest that the potential miRNA regulation of gene expression and the regulatory mechanism of protein expression serve as the basis for future research on pork quality traits in Ningxiang. The ultimate goal is to identify the key dynamic regulatory networks of liver on lipid metabolism and FA synthesis during the growth and development of Ningxiang pigs, and then to inform the biological mechanisms that may affect their performance.

Results

Morphological Characteristics of Liver of Ningxiang Pigs at Different Stages

In order to observe the morphological differences of the liver of Ningxiang pigs during the important period of postnatal growth and development. Liver samples from Ningxiang pigs at the same location were collected at 30, 90, and 210 days of age, respectively. These time points were determined according to the growth characteristics of Ningxiang pigs. Hematoxylin and eosin (H&E) staining showed a gradual increase in connective tissue in the liver at 90 and 210 days of age, and a gradual increase in hepatic lobular sepsis, compared with 30 days of age. Especially at 210 days of age, the hepatic lobules were significantly enlarged and the liver was mature (Figure 1A). These results indicate that a series of biological processes take place in the liver of Ningxiang pigs after birth. Therefore, we combined multiple omics sequencing analyses to systematically study the roles of miRNA, mRNA, and protein in the biological processes of different stages of liver development (Figure 1B).

FIGURE 1

Figure 1. Hematoxylin and eosin staining of spleens and overall study design. (A) H&E staining results of spleens in three developmental stages of Ningxiang pigs. (B) Schematic workflow of the experimental design of this study.

Summary of Transcriptome, miRNAome, and Proteome

All MS/MS spectra were processed using Proteome Discoverer software. Label free analysis of Ningxiang pig liver proteome showed 367,412 spectra in this database and resulted in 5624 identified peptides. In summary, a total of 22,411 unique genes (including 19,653 known mRNAs and 2758 novel mRNAs) and 1122 unique miRNAs (including 384 known miRNAs and 738 novel miRNAs) had been identified (Supplementary Files 1, 2). The quality control analysis of miRNAome can be found in Supplementary File 3. Peptide/protein identification and quantification of all proteins are shown in Supplementary File 4. Protein molecular weight distribution diagram, peptide length distribution diagram, specific peptide number distribution diagram, and protein sequences coverage are shown in Figure 2.

FIGURE 2

Figure 2. Identification and feature of proteins in Ningxiang pigs. (A) Distribution of protein molecular weight. (B) Length distribution of peptide. (C) Number distribution of peptide. (D) Distribution of protein’s sequnences coverage. Each sector represents the proportion of a coverage range. The number outside the sector represents the coverage range and the proportion of proteins distributed in this range.

Screening and Functional Analysis of Differentially Expressed Proteins

A total of 753 differentially expressed proteins (DEPs) were identified in the three groups (Figures 3A–C): 192 DEPs were identified at N90 vs. N30 d, of which 88 were upregulated and 104 were downregulated. N210 vs. N30 d had 268 DEPs, including 72 upregulated DEPs and 196 downregulated DEPs. N210 vs. N90 d had 293 DEPs, among which 75 were up-regulated and 215 were down-regulated. Based on the analysis of the number of DEPs in the three groups, it was found that the number of DEPs in the N90 vs. N30 d group was far less than that in the N210 vs. N90 d group and the N210 vs. N30 d group. The number of DEPs in N210 vs. N90 d and N210 vs. N30 d groups showed the same trend, and most DEPs were decreased at N210 and N90 d. In addition, we identified 29, 55, and 132 DEPs that were uniquely differentially expressed in the N90 vs. N30 d group, N210 vs. N30 d group, and N210 vs. N90 d group, respectively, and 17 DEPs were co-expressed in the three groups (Figure 3D). The protein expression of DEPS in liver tissues of Ningxiang pigs of different ages was analyzed by hierarchical clustering method (Figure 3E). Three biological replicas in each age group showed similar expression trends. In addition, N90 and N210 d DEPs showed similar expression patterns, which indicated that N30 to N90 d was the key stage for the liver growth and development of Ningxiang pigs. In order to further understand the DEPs in three key growth stages involved in regulating the liver function, only the adjacent time point to compare differences in protein Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, we found that compared with 30 days, the 90 days upregulated protein significantly enriched in the 63 regulation pathway, covers almost all of the liver function (Figure 4A). For example, Primary immunodeficiency, Autoimmune thyroid disease, Allograft rejection, and NF-kappa B signaling pathway related to immune defense are extremely significant enrichment. The metabolic pathways related to Glycine, serine, and threonine metabolism, Tryptophan metabolism, Lysine degradation, Metabolism of xenobiotics by cytochrome P450 and FA degradation were significantly enriched. However, compared with 90 days, the 210 days upregulated proteins were only significantly enriched in Protein and absorption (Figure 4B). At the same time, The downregulated DEPs in the N90 vs. N30 d group were enriched only in HIF-1 signaling pathway, Glucagon signaling pathway, ECM-receptor interaction, and Pentose phosphate pathway-related pathways that control cell proliferation and differentiation and maintain cell glucose homeostasis (Figure 4C). However, the number of downregulated pathways in N210 vs. N90 d group was much more than that in N90 vs. N30 d group, which were mainly involved in the control of RNA transcription and protein translation, which indicated that the development of hepatocytes slowed down (Figure 4D and Supplementary File 5). The results showed that the main stage of liver development was from 30 to 90 days of age, during which the liver function might have been gradually matured. Next, in the KEGG enrichment pathway of liver overlapping DEPs, only 14 pathways were significantly identified (P < 0.05, Figure 4E). Among these pathways, lipid metabolism and amino acid metabolism are more abundant, among which proteins Pig.02709.1, Pig.07780.1, Pig.03993.1, and MSTRG.6801_M.6654 are involved in the synthesis and elongation of FAs. Their protein-coding genes are Acyl-CoA synthetase long chain family member 1 (ACSL1), methylglutaconyl-CoA hydratase (AUH), and ethylmalonyl-CoA decarboxylase (ECHDC1), respectively. Although the number and function of overlapping DEPs protein enrichment pathways identified are far less than that of differentially expressed genes (DEGs) (Supplementary File 6), these overlapped co-expressed differential proteins play an irreplaceable role in FA synthesis and metabolic balance during liver development in Ningxiang pigs.

FIGURE 3

Figure 3. Overview of the differentially expressed proteins. (A–C) Differentially expressed protein between different growth stages in postnatal liver. (D) The Venn diagram of DEPs at three comparison groups. (E) Heatmap of proteins expression abundance at 30, 90, and 210 days.

FIGURE 4

Figure 4. Kyoto Encyclopedia of Genes and Genomes analysis of DEPs in the N90 vs. N30 d and N210 vs. N90 d groups. (A) Upregulated KEGG analysis of DE proteins in the N90 vs. N30 d group. (B) Upregulated KEGG analysis of DE proteins in the N210 vs. N90 d group. (C) Downregulated KEGG analysis of DE proteins in the N90 vs. N30 d group. (D) Downregulated KEGG analysis of DE proteins in the N210 vs. N90 d group. The top 20 enriched KEGG pathways ranked by P-values are shown. *P < 0.05, **P < 0.01, ***P < 0.001. (E) Chord diagram of overlapped DEPs enrichment.

Furthermore, the KEGG pathway enrichment analysis was performed for the combination of different genes and proteins in each group (Figure 5). We found that these differential genes and proteins are mainly concentrated in immune defense, antioxidant capacity, metabolism (amino acids, lipids, and carbohydrates), and protein processing. However, the number of differentially enriched genes in pathways between 90 and 210 days of age was far less than that of 90 vs. 30 d and 210 vs. 30 d. But the number of differentially enriched proteins did not change much, suggesting that gene regulation may be more active in the liver at 30 days of age.

FIGURE 5

Figure 5. The number of DEG and DEP enriched in KEGG pathway was analyzed by combined transcriptome and protein analysis.

Liver Proteome Combined Analysis With miRNA and Transcriptome

MicroRNAs are important epigenetic factors that bind to mRNA to inhibit post transcriptional gene expression in eukaryotes. By comparing any two groups, a total of 556 differentially expressed miRNAs (DEMs) were identified, of which 50, 48, and 15 DEMs were uniquely differentially expressed in N30 vs. N90 d group, N30 vs. N210 d group, and N90 vs. N210 d group, respectively. Eleven DEMs were co-expressed in the three groups (Figure 6A). We also found high expression of some miRNAs known to be involved in liver development and FA metabolism, such as ssc-miR-143-3p, ssc-miR-21-5p, and ssc-miR-148a-3p, etc. (Figure 6B).

FIGURE 6

Figure 6. Overview of possible negative miRNA–mRNA–protein interactions. (A) The Venn diagram of DEMs at three comparison groups. (B) Top 20 known miRNAs expressed in liver. (C) STEM analysis of miRNA, transcriptome, and proteome. And count the number and P-value of miRNA/mRNA/protein enriched in each profile. (D) The numbers of possible negative miRNA–mRNA/miRNA–mRNA–protein interactions based on STEM analysis expression patterns. (E) The information of seven potential miRNA–mRNA–protein interaction pairs based on STEM analysis. (F) Part of the miRNA–mRNA negative correlation network based on 2993 profile (III)–profile (IV) expression patterns. Red circles indicate miRNAs, Blue rectangle indicate coding genes.

Correlation analysis of miRNA, mRNA, and protein expression data can effectively identify and evaluate the function of miRNA. In order to explore the dynamic expression patterns of Transcriptomic, miRNA, and Proteomic, we adopted the Short Time-series Expression Miner (STEM) analysis method. We found that Transcriptomic, miRNA, and Proteomic had multiple significantly enriched pattern profiles in the three stages of liver development (Figure 6C, P-value < 0.05). Surprisingly, Time series analysis showed that Transcriptomic, miRNA, and Proteomic were dynamically expressed during liver development, with may be dynamically correlated expression patterns. This implies that there is a functional correlation or regulatory relationship between Transcriptomic, miRNA, and Proteomic. Take into account miRNAs negatively regulate the expression of their target mRNA through target RNA cleavage, it is generally believed that the expression pattern of miRNA target genes is negatively correlated with the expression pattern of miRNAs. In consideration of the three time points (30, 90, and 210 days), we conducted negative regulation analysis in expression patterns to dissect the possible miRNA–mRNA or miRNA–mRNA–protein interactions. We found that profile III and profile IV presented the opposite expression pattern, which may have a negative regulatory relationship. In order to predict the genes targeted by miRNAs, two softwares (miRanda 3.3a and RNAhybrida) were simultaneously used to identify miRNA binding sites. Finally, combine the data predicted by the two algorithms and calculate the overlaps. We found that there were 2993 negative miRNA–mRNA interaction pairs (Figure 6D and Supplementary File 7) in profile III (71)–profile IV (5968), involving 54 miRNAs (including 31 newly identified miRNAs) and 2138 mRNAs. And part of the miRNA–mRNA interactions network is shown in Figure 6F. Combined with STEM analysis of the proteome, we found a total of seven negative miRNA–mRNA–protein interaction pairs (Figure 6E).

The analysis of significantly enriched model profile VII revealed a connected protein network (Figure 7 and Supplementary File 8) divided in: (1) one big module composed by 18 nodes and 36 links, (2) a small component of three singletons. In this network most of the proteins interacted with multiple other partners. However, we found that the proteins uroporphyrinogen decarboxylase (UROD) and lactate dehydrogenase A (LDHA) have the highest degree of connectivity (seven edges), which may enable them to a role as “hub” proteins playing controllers inside biochemical pathways.

FIGURE 7

Figure 7. Protein–protein interaction network analysis based on the profile (VII) for STEM analysis.

Data Verification of DEMs, DEGs, and DEPs

Although these data in this study have been subjected to rigorous statistical and bioinformatics analysis, three biological duplicates were analyzed to eliminate most possible errors. However, we adopted qRT-PCR and PRM method analysis on selected DEMs/DEGs/DEPs to verify the mRNA-seq/miRNA-seq/label-free sequencing results (Figure 8). These miRNAs/mRNAs/proteins were artificially selected as representatives for their possible roles in lipid deposition and FA synthesis. These were found to be correlated with liver lipid metabolism [ssc-miR-7138-5p-PON2 (Fuhrman, 2012; Strakovsky et al., 2014) and Chr12_28849-UROD (Fader et al., 2017)] and liver glycogen catabolism [ssc-miR-6782-3p-FBN1 (Bhadel et al., 2020; Yuan et al., 2020)], disease, and inflammatory immunity [Chr03_9333-CHID1 (Kovaleva et al., 2021)] according to the possible negative miRNA–mRNA–protein interactions on the basis of the STEM analysis during liver development (30–90–210 days) expression patterns.

FIGURE 8

Figure 8. Relative miRNA/mRNA/protein expression of seven selected miRNA–mRNA–protein interaction pairs for comparison of three time periods (30, 90, and 210 days), in respect to miRNAome/transcriptome and qRT-PCR, or proteome and parallel reaction monitoring (PRM). The histograms represent the result of qRT-PCR/PRM verification. Data are represented as mean ± SEM, n = 3 per group. *P < 0.05, **P < 0.01.

The PRM method was used to analyze 15 target proteins from liver samples in three time periods. All the second-order mass spectrum matching figures of the target protein candidate peptides can be found in Supplementary Figure 1. See Supplementary File 9 for all the raw data of Panorama. The data analyses of PRM quantitative Skyline in the target peptide section can be found in Supplementary Figure 2. The results show that quantitative data of target peptide fragments could be obtained from all the samples (three replicates at 30, 90, and 210 days, respectively). The results of qRT-PCR/PRM revealed that the most selected miRNAs/mRNAs/proteins have similar expression tendencies with those from mRNA-seq/miRNA-seq/label-free data when compared with 30 days (P < 0.05).

In addition, 11 negative DE miRNA/DE mRNA involved in key biological processes such as growth and development, glucose metabolism, and FA synthesis, were selected from the 90 vs. 30 d, 90 vs. 210 d, and 210 vs. 30 d groups, respectively, and quantified by qRT-PCR to verify the transcriptome results (Supplementary Figure 3). For example, the synthesis of unsaturated FAs [Chr10_25757-FADS3 and ssc-miR-671-5p-ECHDC1 (Linster et al., 2011; Dewulf et al., 2019)], FA transport [ssc-miR-874-SLC27A3 and Chr03_9333-SLC27A1 (Gallardo et al., 2013)], and FA metabolism [ssc-miR-4332-UCP (Argentato et al., 2018) and ssc-miR-885-3p-CCN4 (Tacke et al., 2018)]. Meanwhile, PRM data were combined to verify the changes of key proteins involved in FA synthesis in the proteome, and serve as Supplementary Figure 4 in the discussion. Even if there are some quantitative differences between these analytical methods, the similarity of the expression tendencies between RNA-seq/miRNA-seq/label-free data and qRT-PCR/PRM indicates that these data is reproducible and reliable.

Discussion

Depending on species and age, the transformation of the hepatic metabolic program during each transitional period of life is a complex molecular regulatory system, many of which are regulated by nuclear receptors and transcription factors, as well as by other dietary and hormonal factors (e.g., carbohydrates, alcohol, insulin). Despite the complexity and degree of coordination of liver lipid mechanisms, any metabolite spillover or nuclear receptor changes will lead to impaired organ function and subsequent pathological changes (Nguyen et al., 2008; Hodson and Gunn, 2020). However, Ningxiang pigs, as a fatty breed, not only have the reputation of excellent meat quality, but also have higher levels of alanine aminotransferase (ALT), ascorbic aminotransferase (AST), and lactate dehydrogenase (LDH-L) in blood than European breeds (Zhao et al., 2012). The liver can be considered a key organ for translating growth rate, feed efficiency, and performance potentials of the animals by exploiting tissue specific metabolic functions. Therefore, the liver of Ningxiang pigs at different growth and development stages is a good model for the study of liver function transformation and lipid metabolism. But because of the influence of the long-term foreign pig breeds, Ningxiang pig population was reduced, especially pure Ningxiang pig were much rarer, and therefore limited understanding of the genetic background of Ningxiang pig. Our current results include detailed information on parallel mRNA, miRNA, and protein expression levels in the liver of pured Ningxiang pigs during the important growth and development periods of lactation (30 days of age), nursery (90 days of age), and fattening (210 days of age).

We identified the expression of 10, 218 DE genes in N90 vs. N30 d, N210 vs. N30 d, and N90 vs. N210 d groups from Ningxiang pig livers. Most DE genes in the three groups were enriched in cancer and diseases (Figure 4). Although “Cancer” does not represent physiological and biochemical reactions for pig, the presence of more differentiated genes from the side was associated with cell proliferation and development (Chen et al., 2017). Notably, the “Rap1signaling pathway,” “PI3K-Akt signaling pathway,” and “MAPK signaling pathway” may provide clues to transboundary regulatory genes involved in signal transduction, signaling molecules, and interactions of cell proliferation and development. Some genes associated with cancer, immune diseases, and cell development have also been found in other animal livers (Li Y. et al., 2016; Zhang et al., 2019). However, there was no complete correlation between the number of DEGs and proteins. For example, most DE proteins in N90 vs. N30 d and N210 vs. N30 d groups were enriched in glucose metabolism, amino acid metabolism, and FA metabolism pathways, which was consistent with the main functions of liver. In addition, through KEGG pathway analysis of liver overlapping proteins, 17 key proteins were identified, among which Pig.02709.1 and MSTRG.6801_M.6654 were involved in FA synthesis, and Pig.07780.1 and Pig.03993.1 were involved in FA elongation. These key overlapping DE genes/proteins may play an irreplaceable role in the functional transformation during the growth and development stage.

The development of mammalian tissues is regulated by various miRNAs. However, many of the miRNAs involved remain unknown. In this study, miRNAs in Ningxiang pig livers from weaning piglets to fattening stage were compared and analyzed by high-throughput sequencing. In total, 384 known miRNAs and 738 new miRNAs were identified and predicted in developing pig livers. To our knowledge, this is the first complete report to systematically identify miRNAs during liver development in Ningxiang pigs using RNA-seq data. Therefore, our findings will enrich the understanding of miRNAs related to FA metabolism, glucose metabolism, and amino acid metabolism in Ningxiang pig liver. The expression level of miRNAs showed obvious temporal dynamic characteristics. In this study, 11 DE miRNAs were coexpressed at three developmental stages. In addition, some miRNAs are expressed at higher levels during liver development. This study presented many known ssc-miR-143-3p, ssc-miR-21-5p, ssc-miR-148a-3p, and many new miRNAs (such as Chr07_18810, Chr07_18928, Chr03_8336, etc.), were abundant miRNAs expressed in the liver of Ningxiang pigs at different time points (Supplementary File 1). It has been reported that exogenous miRNA-143-3p can target the 3′UTR of ANGPTL 8 and regulate the level of high-density lipoprotein cholesterol in hepatocytes (DiStefano, 2019). In addition, miR-21-5p and miR-148a-3p target PDCD4 and ERB 3, respectively, and are involved in regulating liver cell viability and apoptosis (Xiong et al., 2020; Zhang et al., 2020). These identified miRNAs with high expression may play a crucial role in liver development.

As the most widely distributed mammalian FABP, liver fatty acid-binding protein (FABP 1) is highly expressed in the most active liver tissues of long-chain fatty acids (LCFAs) (Jiang et al., 2006). It has been reported that FABP 1 expression is increased during the proliferation and early differentiation of porcine mesenchymal vascular cells, and also regulates the growth and differentiation of embryonic stem cells (Schroeder et al., 2001). Meanwhile, FABP 1 is required for liver uptake and oxidation of LCFAs, and it interacts with PPARA, which regulates FA catabolism (Erol et al., 2004; Landrier et al., 2004). FABP 1 may be responsible for the transport and targeting of FAs to catabolic sites during proliferative interstitial – vascular cell and early adipocyte differentiation. In this study, FABP 1 was significantly higher at 30 days of age than at 90 and 210 days of age, suggesting that FABP 1 plays an important role in liver tissue development and FA metabolism of Ningxiang pigs during lactation. In addition, many other DEPs related to FA metabolism and fat deposition are down-regulated during liver tissue development, such as CD36, HSPA5, FBP, ACSL1, LDH, APMAP, ACADM, etc. FA translocase (CD36) is a single chain transmembrane glycoprotein with a special neck ring structure, which is expressed in a variety of tissues of mammals (Spitsberg et al., 1995; Febbraio et al., 2001). CD36 has been reported to play an important role in adipocyte development and fat deposition, and also act as a receptor for LCFAs, promoting the entry of FAs into cells, thereby promoting the utilization and storage of FAs (Silverstein and Febbraio, 2009; Christiaens et al., 2012). Medium-chain acyl-coenzyme A dehydrogenase (ACADM) is a member of the acyl-coenzyme A dehydrogenase (ACAD) family, which is mainly involved in the β-oxidation process of medium-chain FAs and amino acids (Kim and Miura, 2004). Lack of ACADM hinders the oxidation of FAs and leads to disease or death in mice. This indicates the importance of ACADM in FA metabolism (Smith et al., 2010). APMAP is a 46 kDa glycosylated type II transmembrane protein with an N-terminal anchor and a six-leaf structure of β-an extracellular domain with potential hydrolase activity and calcium binding (Albrektsen et al., 2001; Bogner-Strauss et al., 2010). APMAP is highly upregulated during adipogenic differentiation in various mouse and human cell lines, which is important for in vitro adipogenesis, and the absence of Apmap has a protective effect on diet-induced insulin resistance in mice and shows an overall improvement in glucose utilization in high-fat diets (Pessentheiner et al., 2017). It is speculated that these proteins regulate the synthesis and oxidation of FAs by down-regulating their expression during the growth of Ningxiang pig liver. In addition, many other proteins were significantly up-regulated, such as CD26, AUH, CDH1, EchDDC1, EchDDC2, HSPB1, and PON2, etc. Cadherin-1 (CDH1) is one of two adaptor proteins of the E3 ubiquitin ligase (APC/C) that controls mitosis and DNA replication. So far, the function of CDH1 in vivo has not been fully studied. However, it has been found that CDH1 may be related to the regulation of liver cell cycle. The deletion of CDH1 will upregulate Cyclin D1, lead to the prolongation of S phase and increase DNA replication (Cheng et al., 2017). Weight gain, increased fat mass, and hepatic fat deposition were also observed in CDH1-overexpressing mice. These up-regulated proteins may be closely related to the growth and development of liver cells (Baumeier et al., 2017).

At the same time, we found that the proteins AUH (Nakagawa et al., 1995), ECHDC1, and ACSL1 were differentially expressed in any two groups of liver comparison. ACSL1 is a member of the five known acyl-CoA synthases (ACSLs). Specific loss of ACSL1 in the liver reduces triglyceride synthesis and β-oxidation of LCFAs. And changes the composition of phospholipid FAs, while overexpression of ACSL1 can increase the content of TG in liver and blood, thus indirectly affecting meat quality (Parkes et al., 2006). It is well known that most FAs are straight-chain FAs synthesized by FASN using acetyl-CoA and malonyl CoA. However, the absence of ECHDC1 in adipocytes restricts the accumulation of cytoplasmic methyl and ethylmalonyl-coenzyme A, leading to a significant reduction in the synthesis of linear FAs, indicating that ECHDC1 may serve as a key modulator of methyl branched FA expression (Dewulf et al., 2019). These key overlapping DE proteins may play an irreplaceable role in the functional transformation during growth and development.

Time series analysis showed that protein-coding genes and miRNAs were expressed dynamically during liver development, with opposite expression patterns. These results indicate that this implies a functional correlation or regulatory relationship between Transcriptomic and miRNA. Since miRNA negatively regulates the expression of its target mRNA through cleavage of target RNA, which provides a basis for protein changes to some extent, the expression pattern of miRNA target genes is generally negatively correlated with the expression pattern of miRNA. For example, ssc-miR-4332-UCP3, Chr08_21275-ACACA, ssc-miR-885-3p-CCN4, ssc-miR-4332-ATP2A3, ssc-miR-874-LDLRAD4, ssc-miR-138-LDLRAD3, and other miRNA–mRNA pairs related to lipid anabolism are in line with the negative regulation law of STEM analysis, which can provide ideas for future research on liver function. From a biological point of view, genes/proteins with similar expression patterns may have common characteristics, such as being regulated by a certain gene at the same time, having similar biological functions, or having a common cellular origin. By STEM analysis, we also observed similar expression trends between mRNAs (profile IV) and proteins (profile VII). Then, for the functional connections established between these proteins with similar expression trends, we found that most of the proteins in this network have multiple other partners interacting with each other. In addition, the protein UROD presented the highest degree of connection (seven edges), which may assign to them a role as “hub” proteins playing a putative function of controllers inside biochemical pathways. It has been reported that UROD is a cytoplasmic enzyme, which is related to heme production. Heme is essential for all cells and involved in the oxidation of a variety of chemicals, including drugs, vitamins, FAs, and prostaglandins (Bonkowsky et al., 1979; Shimizu et al., 2019). Furthermore, proteins related to lipid metabolism, such as superoxide dismutase 1 (SOD1), LDHA, and paraoxonase 2 (PON2), were also included in this network. Combined with miRNAome data, we found seven negative miRNA–mRNA–protein interactions pairs (including protein UROD and PON2) from profile III–profile IV–profile VII, indicating that these proteins with similar expression profiles may be regulated by certain genes and participate in the same regulatory pathway. These interaction networks are expected to lay a foundation for the future study of the molecular regulation mechanism of lipid metabolism in Ningxiang pig liver at different periods.

Conclusion

In this study, RNA-seq, miRNA-seq, label-Free, and PRM were used to analyze the transcriptome and proteome of liver tissues at three key developmental time points (lactation period, nursery period, and fattening period) in Ningxiang pigs. The key dynamic regulatory network related to liver development and FA synthesis and metabolism of Ningxiang pigs was identified. Taken together, liver function is closely related to metabolic and digestive diseases of animals and meat quality traits of livestock. although this study cannot fully explain all the molecular mechanisms in regulation of lipid metabolism in pig liver. Research results from transcription, miRNA, and proteome analysis provide initial database in elucidating the genetic mechanisms underlying diseases and meat quality traits in livestock. It also provides a theoretical basis for future breeding applications.

Materials and Methods

Experimental Animals and Sampling

Nine purebred male Ningxiang pigs were selected for this experiment, including 30 days (30 days after birth, n = 3), 90 days (90 days after birth, n = 3), and 210 days (210 days after birth, n = 3). Three groups of pig are half-sibs, and the duplicate samples are full-sibs. According to standard procedures, three healthy Ningxiang pigs of similar weight were selected from three age groups. The pigs are provided by the Ningxiang Pig Farm of Hunan Dalong Animal Husbandry Technology Co., Ltd. Under standard environmental conditions, they were fed three times routinely with free drinking water. Samples of the left liver lobules were collected within 30 min after slaughter, frozen in liquid nitrogen, and stored at −80°C. These harvested samples were used to conduct mRNA-seq, miRNA-seq, proteomics, qRT-PCR, and PRM, respectively. All the above experiments including the mass spectrometry experiments had three biological replicates (N30d_1, N30d_2, N30d_3, N90d_1,N90d_2, N90d_3, N210d_1, N210d_2, and N210d_3).

Hematoxylin–Eosin Staining

The liver tissue is taken out from the fixative, and the dehydration machine is dehydrated with gradient alcohol. Melt the paraffin at 65° for 3 h. Embed the wax-soaked liver tissue in the embedding machine. The wax block is sliced in a paraffin microtome with a thickness of 4 μm. After roasting the wax, take it out at room temperature and store it for later use. Stain with hematoxylin for 3–5 min, wash with water. The sections were dehydrated with 85 and 95% gradient alcohol in sequence, and stained with eosin staining solution for 5 min. Seal the slices with neutral gum. Microscope inspection, Case Viewer 2.4 image analysis software collects and analyses pictures.

RNA Isolation, Library Preparation, and Transcriptome Sequencing

According to the TruSeq^TM strand total RNA kit from Illumina (San Diego, CA, United States), 5 μg total RNA was used to prepare an RNA-seq transcriptome strand library. After that, the ribosomal RNA (rRNA) was removed using the Ribo-Zero Magnetic kit. And fragmented by fragmentation buffer firstly. Second, random hexamer primers were used to synthesize first-strand complementary DNA (cDNA). Then the RNA template was removed and a replacement strand was synthesized, and dUTP was used instead of dTTP to generate ds cDNA. AMPure XP beads were used to separate ds cDNA from the second-strand reaction mix. An “A” nucleotide was added to the 3′ ends to prevent them from ligating to one another during the adapter ligation reaction. Finally, multiple indexing adaptors were ligated to the ends of the ds cDNA (Cui et al., 2018). After quantification by TBS380, Shanghai Majorbio Bio-Pharm Biotechnology Co., Ltd. (Shanghai, China) provides sequencing technical support.

Transcripts Assembly

The quality control of raw paired-end reads were performed through the default parameters of SeqPrep.¹ Then the HISAT2² software was used to compare the clean reads with the Ningxiang pig reference genome (accession number PPJNA531381) in a orientation mode. The mapped readings of each sample were assembled by StringTie³ in a reference-based method. After comparing with HISAT2, it is assembled by stringtie software. Compare with the known transcripts in the known genome annotation information (Gffcompare software) to obtain new transcripts, and select the classcode as “x,” “i,” “j,” “u,” “o” transcript. Transcripts with length ≥200 bp, Exon number ≥2, and ORF length ≤300 bp were selected as candidate lncRNAs for preliminary screening. Whether a transcript has coding potential is a key condition for judging whether a transcript is an lncRNA. The most widely used coding potential analysis method for lncRNA screening mainly includes four methods: CPC analysis, CNCI analysis, CPAT analysis, and Pfam protein domain analysis. The default screening conditions for these four methods are: CPC score <0.5; CNCI score <0; CPAT score <0.5; Pfam database has no annotations. The overlap analysis of lncRNA predicted by these four methods was carried out. Finally, the result of the new prediction of lncRNA takes the intersection of these softwares as the basis for subsequent analysis. Those with coding capabilities that are not predicted to be lncRNAs will be classified as new mRNAs.

Data Construction and Processing of miRNA Libraries

One microgram total RNA was extracted from the tissue as the starting template. The small RNA sequencing library was prepared using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, United States). The random primer was reversed into first cDNA (TruSeQ^TM Small RNA Sample Prep Kit). Library purification to remove adaptor dimer, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeat sequences (6% Novex TBE PAGE Gel, 1.0 mm, 10 well). TBS380 (Picogreen) quantitative, according to the proportion of data mixed on the computer. The constructed library was sequenced using Illumina Hiseq 2000/2500, and the length of sequenced reads was 1 × 50 bp with single end. Compare clean reads with miRNA precursor and mature body sequences in the miRBase database of selected wild boar species. Align the sRNA sequences that cannot be compared on Rfam and miRBase to the reference genome, intercept the surrounding sequences and use miRDeep2 software to predict the secondary structure. According to the prediction results, use Dicer restriction site information, energy value, and other characteristics to filter to identify new miRNAs.

The Prediction of Target Genes of miRNAs

In order to predict the genes targeted by miRNAs, we used two computational target prediction algorithms, miRanda 3.3a (score cutoff ≥160 and energy cutoff 20), and RNAHybrida (energy cutoff = 20) to jointly identify the binding sites of miRNAs. Finally, the data predicted by the two algorithms are combined and the overlapping target genes are calculated.

Analysis of Differentially Expressed Genes and Differentially Expressed miRNAs

Salmon (Patro et al., 2017) was used to perform expression level for unigenes by calculating Transcripts Per Million (TPM) (Mortazavi et al., 2008). DEMs and DEGs were identified by the DEseq2 with p-adjust < 0.05 and had fold-change ≥2 with N90 vs. N30 d, N210 vs. N30 d and N210 vs. N90 d. DEPs were identified by the t-test with P-value < 0.05 and had fold-change ≥1.2 or ≤0.83 with N90 vs. N30 d, N210 vs. N30 d, and N210 vs. N90 d. Gene function was annotated based on the following databases: NCBI non-redundant protein sequences (Nr,ftp⁴), Protein Family (Pfam⁵), String,⁶ and (KEGG⁷).

Time-Series Analysis

Time-series analysis was performed by STEM (Ernst and Bar-Joseph, 2006). The STEM analysis studies the dynamic behavior of gene expression and measures a series of processes that have a strong correlation with time points. The significantly enriched model profiles have a P-value of less than 0.05.

Prediction of Target Genes of miRNAs, and Network Analysis

To define all the possible negative relationships between miRNA and mRNA expression. RNAhybrid and miRanda were used to predict the miRNA-gene pairs. The correlation between mRNAs and miRNAs was evaluated using the Pearson correlation coefficient (PCC) from matched mRNA and miRNA expression profile data. The interaction network was built and visually displayed using Cytoscape software. Protein identification was employed against the above RNA-seq results for database searching, the proteins added miRNA–mRNA pairs directly to the miRNA–mRNA–protein (Morris et al., 2015).

Mass Spectrometry Analysis

All samples were analyzed using Q Exactive HF-X mass spectrometer (Thermo, United States) and UltiMate 3000 RSLC nanosystem (Thermo, United States). The peptide was injected onto a self-made C18 column (75 μm × 25mm, Thermo, United States) in buffer A (2% acetonitrile and 0.1% Formic acid) and separated with a linear gradient of buffer B (80% acetonitrile and 0.1% Formic acid) at a flow rate of 300 nl/min. Through nanoelectrospray ionization and Q-Exactive HF-X quadrupole orbit rap mass spectrometer interface, the integrated column heater is set at 50°C, and the electrospray voltage is 1.8 kV. The tandem mass spectra of the first 20 ions were scanned with a full range of 350–1300 m/z, the isolation width was set to 1.6 Da, the resolution was 70K, and the MS/MS scan resolution was 17.5K. In all cases, one microscan was recorded using dynamic exclusion of 18 s.

Sequence Database Selection and Searching

Proteome Discoverer TM Software 2.2 software was used for database search identification and quantitative analysis against the Ningxiang pig database (accession number PPJNA531381). Since the original mass spectrometer off-machine data was a raw file, the raw file is submitted to the Proteome Discoverer server, and the Ningxiang pig genome database is selected for database search. The parameters for protein searching were the default parameter. Trypsin digestion with up to two missed cleavages. The carbamoylmethylation of cysteine is used as a fixed modification, and the oxidation of methionine and the acetylation of the protein N-terminus were used as variable modifications. Based on the q-value, a 1% false discovery rate (FDR) was used to verify the peptide spectral matching.

Predicted miRNA–mRNA–Protein Validation by qRT-PCR and PRM

To validate the identified RNAs in Ningxiang pigs, total RNA was extracted using an Animal Total RNA Kit (Tiangen, China). Total RNA was extracted from the three stages (30, 90, and 210 days) livers used in the RNA-seq and reversed to cDNA with a Primescript RT Master Kit (Thermo Scientific, United States) with random primers in accordance with the manufacturer’s instructions. The PCR amplicon template sequence described above was then used to design the PCR primers using Primer 5 software. Primer sequences used in the qRT-PCR are listed in Supplementary File 10. All reactions were performed in triplicate for each sample. The relative expression level of the miRNAs and mRNAs was calculated via the comparative CT method (2^{–△△ CT}). U6/GAPDH was used as an internal control.

Mass spectrometric data was collected using the Triple TOF5600 and liquid mass spectrometry system (SCIEX). For PRM data collection, the optimized protein peptide list was added to the inclusion list. The peptides in the list were selected and separated one by one by mass spectrometry, and the secondary ion spectrum was collected. First, a part of all the peptide samples obtained after trypsin digestion were mixed for mass spectrometric DDA detection. The resulting MS–MS data were processed by Protein pilot^TM V4.5 search engine, then were processed with Skyline software and imported into the corresponding spectra library (MacLean et al., 2010). Screening the target protein peptide and establishing the method. The PRM detection method was designed by adding the target peptide m/z to the inclusion list to establish a mass spectrometry acquisition method. The mixed sample’s PRM data was collected by the PRM method to adjust and optimize the PRM data acquisition method and form the final PRM method for sample data acquisition. The optimized PRM mass spectrometry acquisition method was used to collect data for each sample, and PRM spectrum files was obtained. Through the extraction and analysis of these PRM map files, the quantitative information of protein can be obtained.

Statistical Analysis

An unpaired two-tailed Student’s t-test was used to test changes in candidate differential miRNAs and mRNAs levels between any two time points. ^∗P < 0.05, ^∗∗P < 0.01. Data were expressed as mean ± SEM.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by all animal experiments in this study were approved by the “Institutional Animal Care and Use Committee of Hunan Agricultural University” (Changsha, China).

Author Contributions

BL performed the data analysis and wrote the manuscript. QZ and JiH provided the site of experiment and some of test conditions support. YX and YG conducted part of animal experiments. YX completed the animal feeding experiment. YD, JY, KX, JuH, and HM designed this study and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by Genetic resource innovations and its application of indigenous pig breeds in Hunan province (2019NK2193); “double first Class” Construction Project of Hunan Agricultural University (kxk201801004); Strategic Priority Research Program of Chinese Academy of Sciences (XDA24030204); National Natural Science Foundation of China (U19A2037); Basic Research Project of Ningxiang Pig (15073); Hunan Natural Science Joint Fund (2018JJ4003), and Hunan Science and Technology Innovation Project (72010102-30218).

Conflict of Interest

QZ is employed by Dalong Livestock Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.709521/full#supplementary-material

Supplementary Figure 1 | All second-order mass spectrum matching figures of target protein candidate peptides.

Supplementary Figure 2 | Skyline analysis of target peptide PRM results.

Supplementary Figure 3 | Validation of the expression of DEGs/DEMs by qRT-PCR. Data from qRT-PCR are shown as column and Y-axis on the left, while the data from RNA-seq are shown as line and Y-axis on the right. Data are represented as mean ± SEM, n = 3 per group.

Supplementary Figure 4 | Validation of the expression of DEPs by parallel reaction monitoring (PRM).

Supplementary File 1 | List of expressed genes and miRNA during Ningxiang pig liver development. Gene ID, TPM value, gene name, GO/KO ID, position, and description of all acquired genes.

Supplementary File 2 | List of all miRNA sequences identified during Ningxiang pig liver development.

Supplementary File 3 | Quality control analysis of the miRNAome.

Supplementary File 4 | List of expressed proteins during Ningxiang pig liver development. Protein ID, abundance value, description of all acquired proteins.

Supplementary File 5 | Kyoto Encyclopedia of Genes and Genomes pathways for upregulated and downregulated DEPs from N90 vs. N30 d and N210 vs. N90 d. Pathways marked gray background represent significant enrichment, while other terms represent non-significant enrichment.

Supplementary File 6 | Kyoto Encyclopedia of Genes and Genomes pathways were identified by KEGG pathway analysis of hepatic overlapping differentially expressed (DE) mRNAs from three comparison groups.

Supplementary File 7 | The relationship of miRNAs and mRNAs regulatory networks based on profile III–profile IV.

Supplementary File 8 | List of protein interaction network node attribute.

Supplementary File 9 | Data analysis of PRM quantitative in the target peptide section. The yellow bottom represents the internal control protein alpha-actinin-1.

Supplementary File 10 | Primers of quantitative real-time RT-PCR.

Abbreviations

ncRNA, non-coding RNA; miRNA, microRNA; N30 d, Ningxiang pig 30 days after birth; N90 d, Ningxiang pig 90 days after birth; N210 d, Ningxiang pig 210 days after birth; vs, versus; DEGs, differentially expressed mRNAs; DEMs, differentially expressed miRNAs; DEPs, differentially expressed proteins; FDR, false discovery rate; TPM, Transcripts Per Kilobase of exon model per Million mapped reads; KEGG, Kyoto Encyclopedia of Genes and Genomes; STEM, Short Time-series Expression Miner; RPM, reads per million mapped reads; RNA-seq, RNA sequencing; CIRI 2, CircRNA identifier 2; qRT-PCR, quantitative real-time RT-PCR; PRM, parallel reaction monitor.

Footnotes

1. ^ https://github.com/OpenGene/fastp
2. ^ https://ccb.jhu.edu/software/hisat2/index.shtml
3. ^ https://ccb.jhu.edu/software/stringtie/index.shtml
4. ^ http://ftp.ncbi.nlm.nih.gov/blast/db/
5. ^ http://pfam.xfam.org/
6. ^ https://string-db.org
7. ^ http://www.genome.jp/kegg/

References

Albrektsen, T., Richter, H. E., Clausen, J. T., and Fleckner, J. (2001). Identification of a novel integral plasma membrane protein induced during adipocyte differentiation. Biochem. J. 359(Pt 2), 393–402. doi: 10.1042/0264-6021:3590393

CrossRef Full Text | Google Scholar

Argentato, P. P., de Cássia César, H., Estadella, D., and Pisani, L. P. (2018). Programming mediated by fatty acids affects uncoupling protein 1 (UCP-1) in brown adipose tissue. Br. J. Nutr. 120, 619–627. doi: 10.1017/s0007114518001629

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassols, A., Costa, C., Eckersall, P. D., Osada, J., Sabria, J., and Tibau, J. (2014). The pig as an animal model for human pathologies: a proteomics perspective. Proteomics Clin. Appl. 8, 715–731. doi: 10.1002/prca.201300099

PubMed Abstract | CrossRef Full Text | Google Scholar

Baumeier, C., Schlüter, L., Saussenthaler, S., Laeger, T., Rödiger, M., Alaze, S. A., et al. (2017). Elevated hepatic DPP4 activity promotes insulin resistance and non-alcoholic fatty liver disease. Mol. Metab. 6, 1254–1263. doi: 10.1016/j.molmet.2017.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhadel, P., Shrestha, S., Sapkota, B., Li, J. Y., and Tao, H. (2020). Asprosin and type 2 diabetes mellitus: a novel potential therapeutic implication. J. Biol. Regul. Homeost. Agents 34, doi: 10.23812/19-244-e [Epub ahead of print].

CrossRef Full Text | Google Scholar

Bogner-Strauss, J. G., Prokesch, A., Sanchez-Cabo, F., Rieder, D., Hackl, H., Duszka, K., et al. (2010). Reconstruction of gene association network reveals a transmembrane protein required for adipogenesis and targeted by PPARγ. Cell. Mol. Life Sci. 67, 4049–4064. doi: 10.1007/s00018-010-0424-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bonkowsky, H. L., Sinclair, P. R., and Sinclair, J. F. (1979). Hepatic heme metabolism and its control. Yale J. Biol. Med. 52, 13–37.

Google Scholar

Bovo, S., Di Luca, A., Galimberti, G., Dall’Olio, S., and Fontanesi, L. (2018). A comparative analysis of label-free liquid chromatography-mass spectrometry liver proteomic profiles highlights metabolic differences between pig breeds. PLoS One 13:e0199649. doi: 10.1371/journal.pone.0199649

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Zhang, Y. H., Lu, G., Huang, T., and Cai, Y. D. (2017). Analysis of cancer-related lncRNAs using gene ontology and KEGG pathways. Artif. Intell. Med. 76, 27–36. doi: 10.1016/j.artmed.2017.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, R., Liang, X., Zhao, Q., Lian, Z., Tang, L., Qiu, C., et al. (2017). APC(Cdh1) controls cell cycle entry during liver regeneration. Exp. Cell Res. 354, 78–84. doi: 10.1016/j.yexcr.2017.03.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Christiaens, V., Van Hul, M., Lijnen, H. R., and Scroyen, I. (2012). CD36 promotes adipocyte differentiation and adipogenesis. Biochim. Biophys. Acta 1820, 949–956. doi: 10.1016/j.bbagen.2012.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, J., Xu, P., Meng, J., Li, J., Jiang, N., and Luan, Y. (2018). Transcriptome signatures of tomato leaf induced by Phytophthora infestans and functional identification of transcription factor SpWRKY3. Theor. Appl. Genet. 131, 787–800. doi: 10.1007/s00122-017-3035-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, Y. J., Hao, Y., Li, J. L., Bao, W. G., Li, G., Gao, Y. L., et al. (2016). Chronic heat stress induces immune response, oxidative stress response, and apoptosis of finishing pig liver: a proteomic approach. Int. J. Mol. Sci. 17:393. doi: 10.3390/Ijms17050393

PubMed Abstract | CrossRef Full Text | Google Scholar

Dewulf, J. P., Gerin, I., Rider, M. H., Veiga-da-Cunha, M., Van Schaftingen, E., and Bommer, G. T. (2019). The synthesis of branched-chain fatty acids is limited by enzymatic decarboxylation of ethyl- and methylmalonyl-CoA. Biochem. J. 476, 2427–2447. doi: 10.1042/bcj20190500

PubMed Abstract | CrossRef Full Text | Google Scholar

DiStefano, J. K. (2019). Angiopoietin-like 8 (ANGPTL8) expression is regulated by miR-143-3p in human hepatocytes. Gene 681, 1–6. doi: 10.1016/j.gene.2018.09.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Ernst, J., and Bar-Joseph, Z. (2006). STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics 7:191. doi: 10.1186/1471-2105-7-191

PubMed Abstract | CrossRef Full Text | Google Scholar

Erol, E., Kumar, L. S., Cline, G. W., Shulman, G. I., Kelly, D. P., and Binas, B. (2004). Liver fatty acid binding protein is required for high rates of hepatic fatty acid oxidation but not for the action of PPARalpha in fasting mice. FASEB J. 18, 347–349. doi: 10.1096/fj.03-0330fje

PubMed Abstract | CrossRef Full Text | Google Scholar

Fader, K. A., Nault, R., Kirby, M. P., Markous, G., Matthews, J., and Zacharewski, T. R. (2017). Convergence of hepcidin deficiency, systemic iron overloading, heme accumulation, and REV-ERBα/β activation in aryl hydrocarbon receptor-elicited hepatotoxicity. Toxicol. Appl. Pharmacol. 321, 1–17. doi: 10.1016/j.taap.2017.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Febbraio, M., Hajjar, D. P., and Silverstein, R. L. (2001). CD36: a class B scavenger receptor involved in angiogenesis, atherosclerosis, inflammation, and lipid metabolism. J. Clin. Invest. 108, 785–791. doi: 10.1172/jci14006

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuhrman, B. (2012). Regulation of hepatic paraoxonase-1 expression. J. Lipids 2012:684010. doi: 10.1155/2012/684010

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallardo, D., Amills, M., Quintanilla, R., and Pena, R. N. (2013). Mapping and tissue mRNA expression analysis of the pig solute carrier 27A (SLC27A) multigene family. Gene 515, 220–223. doi: 10.1016/j.gene.2012.11.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y. J., Yu, J. J., Wang, C. X., Li, K., Liu, B., Du, Y., et al. (2017). miR-212-5p suppresses lipid accumulation by targeting FAS and SCD1. J. Mol. Endocrinol. 59, 205–217. doi: 10.1530/Jme-16-0179

PubMed Abstract | CrossRef Full Text | Google Scholar

He, J., Wu, X. L., Zeng, Q., Li, H., Ma, H., Jiang, J., et al. (2020). Genomic mating as sustainable breeding for Chinese indigenous Ningxiang pigs. PLoS One 15:e0236629. doi: 10.1371/journal.pone.0236629

PubMed Abstract | CrossRef Full Text | Google Scholar

Hicks, J. A., Porter, T. E., Sunny, N. E., and Liu, H. C. (2019). Delayed feeding alters transcriptional and post-transcriptional regulation of hepatic metabolic pathways in peri-hatch broiler Chicks. Genes 10:272. doi: 10.3390/Genes10040272

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodson, L., and Gunn, P. J. (2020). The regulation of hepatic fatty acid synthesis and partitioning: the effect of nutritional state (vol 15, pg 689, 2019). Nat. Rev. Endocrinol. 16, 340–340. doi: 10.1038/s41574-020-0345-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Iliopoulos, D., Drosatos, K., Hiyama, Y., Goldberg, I. J., and Zannis, V. I. (2010). MicroRNA-370 controls the expression of microRNA-122 and Cpt1alpha and affects lipid metabolism. J. Lipid Res. 51, 1513–1523. doi: 10.1194/jlr.M004812

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y. Z., Li, X. W., and Yang, G. X. (2006). Sequence characterization, tissue-specific expression and polymorphism of the porcine (Sus scrofa) liver-type fatty acid binding protein gene. Yi Chuan Xue Bao 33, 598–606. doi: 10.1016/s0379-4172(06)60089-9

CrossRef Full Text | Google Scholar

Kiamehr, M., Heiskanen, L., Laufer, T., Dusterloh, A., Kahraman, M., Kakela, R., et al. (2019). Dedifferentiation of primary hepatocytes is accompanied with reorganization of lipid metabolism indicated by altered molecular lipid and miRNA profiles. Int. J. Mol. Sci. 20:2910. doi: 10.3390/Ijms20122910

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. J., and Miura, R. (2004). Acyl-CoA dehydrogenases and acyl-CoA oxidases. Structural basis for mechanistic similarities and differences. Eur. J. Biochem. 271, 483–493. doi: 10.1046/j.1432-1033.2003.03948.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kovaleva, O. V., Rashidova, M. A., Samoilova, D. V., Podlesnaya, P. A., Tabiev, R. M., Mochalnikova, V. V., et al. (2021). CHID1 is a novel prognostic marker of non-small cell lung cancer. Int. J. Mol. Sci. 22:450. doi: 10.3390/ijms22010450

PubMed Abstract | CrossRef Full Text | Google Scholar

Landrier, J. F., Thomas, C., Grober, J., Duez, H., Percevault, F., Souidi, M., et al. (2004). Statin induction of liver fatty acid-binding protein (L-FABP) gene expression is peroxisome proliferator-activated receptor-alpha-dependent. J. Biol. Chem. 279, 45512–45518. doi: 10.1074/jbc.M407461200

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Ma, Z., Jia, L. J., Li, Y. M., Xu, C. L., Wang, T. A., et al. (2016). Systematic analysis of the regulatory functions of microRNAs in chicken hepatic lipid metabolism. Sci. Rep. 6:31766. doi: 10.1038/Sre31766

CrossRef Full Text | Google Scholar

Li, J., Levitan, B., Gomez-Jimenez, S., and Kültz, D. (2018). Development of a gill assay library for ecological proteomics of threespine sticklebacks (Gasterosteus aculeatus). Mol. Cell Proteomics 17, 2146–2163. doi: 10.1074/mcp.RA118.000973

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Li, X. C., Sun, W. K., Cheng, C., Chen, Y. H., Zeng, K., et al. (2016). Comparison of liver microRNA transcriptomes of Tibetan and Yorkshire pigs by deep sequencing. Gene 577, 244–250. doi: 10.1016/j.gene.2015.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Linster, C. L., Noël, G., Stroobant, V., Vertommen, D., Vincent, M. F., Bommer, G. T., et al. (2011). Ethylmalonyl-CoA decarboxylase, a new enzyme involved in metabolite proofreading. J. Biol. Chem. 286, 42992–43003. doi: 10.1074/jbc.M111.281527

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C. D., Shen, L. Y., Du, J. J., Wu, X. Q., Luo, J., Pu, Q., et al. (2018). The effect of lipid metabolism-related genes on intramuscular fat content and fatty acid composition in multiple muscles. Anim. Prod. Sci. 58, 2003–2010. doi: 10.1071/AN16292

CrossRef Full Text | Google Scholar

MacLean, B., Tomazela, D. M., Shulman, N., Chambers, M., Finney, G. L., Frewen, B., et al. (2010). Skyline: an open source document editor for creating and analyzing targeted proteomics experiments. Bioinformatics 26, 966–968. doi: 10.1093/bioinformatics/btq054

PubMed Abstract | CrossRef Full Text | Google Scholar

Manzoni, C., Kia, D. A., Vandrovcova, J., Hardy, J., Wood, N. W., Lewis, P. A., et al. (2018). Genome, transcriptome and proteome: the rise of omics data and their integration in biomedical sciences. Brief. Bioinform. 19, 286–302. doi: 10.1093/bib/bbw114

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, J. H., Vijay, D., Federowicz, S., Pico, A. R., and Ferrin, T. E. (2015). CyAnimator: simple animations of cytoscape networks. F1000Res 4:482. doi: 10.12688/f1000research.6852.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakagawa, J., Waldner, H., Meyer-Monard, S., Hofsteenge, J., Jenö, P., and Moroni, C. (1995). AUH, a gene encoding an AU-specific RNA binding protein with intrinsic enoyl-CoA hydratase activity. Proc. Natl. Acad. Sci. U. S. A. 92, 2051–2055. doi: 10.1073/pnas.92.6.2051

PubMed Abstract | CrossRef Full Text | Google Scholar

Neilson, K. A., Ali, N. A., Muralidharan, S., Mirzaei, M., Mariani, M., Assadourian, G., et al. (2011). Less label, more free: approaches in label-free quantitative mass spectrometry. Proteomics 11, 535–553. doi: 10.1002/pmic.201000553

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, P., Leray, V., Diez, M., Serisier, S., Le Bloc’h, J., Siliart, B., et al. (2008). Liver lipid metabolism. J. Anim. Physiol. Anim. Nutr. (Berl.). 92, 272–283. doi: 10.1111/j.1439-0396.2007.00752.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Parkes, H. A., Preston, E., Wilks, D., Ballesteros, M., Carpenter, L., Wood, L., et al. (2006). Overexpression of acyl-CoA synthetase-1 increases lipid deposition in hepatic (HepG2) cells and rodent liver in vivo. Am. J. Physiol. Endocrinol. Metab. 291, E737–E744. doi: 10.1152/ajpendo.00112.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419. doi: 10.1038/nmeth.4197

PubMed Abstract | CrossRef Full Text | Google Scholar

Pessentheiner, A. R., Huber, K., Pelzmann, H. J., Prokesch, A., Radner, F. P. W., Wolinski, H., et al. (2017). APMAP interacts with lysyl oxidase-like proteins, and disruption of Apmap leads to beneficial visceral adipose tissue expansion. FASEB J. 31, 4088–4103. doi: 10.1096/fj.201601337R

PubMed Abstract | CrossRef Full Text | Google Scholar

Peterson, A. C., Russell, J. D., Bailey, D. J., Westphall, M. S., and Coon, J. J. (2012). Parallel reaction monitoring for high resolution and high mass accuracy quantitative, targeted proteomics. Mol. Cell. Proteomics 11, 1475–1488. doi: 10.1074/mcp.O112.020131

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroeder, F., Atshaves, B. P., Starodub, O., Boedeker, A. L., Smith, R. R. III, Roths, J. B., et al. (2001). Expression of liver fatty acid binding protein alters growth and differentiation of embryonic stem cells. Mol. Cell. Biochem. 219, 127–138. doi: 10.1023/a:1010851130136

CrossRef Full Text | Google Scholar

Shao, F., Wang, X., Yu, J., Shen, K., Qi, C., and Gu, Z. (2019). Expression of miR-33 from an SREBP2 intron inhibits the expression of the fatty acid oxidation-regulatory genes CROT and HADHB in chicken liver. Br. Poult. Sci. 60, 115–124. doi: 10.1080/00071668.2018.1564242

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimizu, T., Lengalova, A., Martínek, V., and Martínková, M. (2019). Heme: emergent roles of heme in signal transduction, functional regulation and as catalytic centres. Chem. Soc. Rev. 48, 5624–5657. doi: 10.1039/c9cs00268e

PubMed Abstract | CrossRef Full Text | Google Scholar

Silverstein, R. L., and Febbraio, M. (2009). CD36, a scavenger receptor involved in immunity, metabolism, angiogenesis, and behavior. Sci. Signal. 2:re3. doi: 10.1126/scisignal.272re3

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, E. H., Thomas, C., McHugh, D., Gavrilov, D., Raymond, K., Rinaldo, P., et al. (2010). Allelic diversity in MCAD deficiency: the biochemical classification of 54 variants identified during 5 years of ACADM sequencing. Mol. Genet. Metab. 100, 241–250. doi: 10.1016/j.ymgme.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Spitsberg, V. L., Matitashvili, E., and Gorewit, R. C. (1995). Association and coexpression of fatty-acid-binding protein and glycoprotein CD36 in the bovine mammary gland. Eur. J. Biochem. 230, 872–878. doi: 10.1111/j.1432-1033.1995.tb20630.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Strakovsky, R. S., Zhang, X., Zhou, D., and Pan, Y. X. (2014). The regulation of hepatic Pon1 by a maternal high-fat diet is gender specific and may occur through promoter histone modifications in neonatal rats. J. Nutr. Biochem. 25, 170–176. doi: 10.1016/j.jnutbio.2013.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Tacke, C., Aleksandrova, K., Rehfeldt, M., Murahovschi, V., Markova, M., Kemper, M., et al. (2018). Assessment of circulating Wnt1 inducible signalling pathway protein 1 (WISP-1)/CCN4 as a novel biomarker of obesity. J. Cell Commun. Signal. 12, 539–548. doi: 10.1007/s12079-017-0427-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Uni, Z., Ferket, P. R., Tako, E., and Kedar, O. (2005). In ovo feeding improves energy status of late-term chicken embryos. Poult. Sci. 84, 764–770. doi: 10.1093/ps/84.5.764

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, J., Ni, J. B., Chen, C. Y., and Wang, K. Z. (2020). miR-148a-3p regulates alcoholic liver fibrosis through targeting ERBB3. Int. J. Mol. Med. 46, 1003–1012. doi: 10.3892/ijmm.2020.4655

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, M. Y., Li, W. D., Zhu, Y., Yu, B. Y., and Wu, J. (2020). Asprosin: a novel player in metabolic diseases. Front. Endocrinol. 11:64. doi: 10.3389/fendo.2020.00064

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, G., Li, J., Zhang, J., Liang, X., Zhang, X., Wang, T., et al. (2019). Integrated analysis of transcriptomic, mirna and proteomic changes of a novel hybrid yellow catfish uncovers key roles for miRNAs in Heterosis. Mol. Cell. Proteomics 18, 1437–1453. doi: 10.1074/mcp.RA118.001297

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Li, C. C., Li, F., Li, H., Liu, X. J., Loor, J. J., et al. (2017). Estrogen promotes hepatic synthesis of long-chain polyunsaturated fatty acids by regulating ELOVL5 at Post-transcriptional level in laying hens. Int. J. Mol. Sci. 18:1405. doi: 10.3390/Ijms18071405

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M. Y., Sun, W. L., Zhou, M. H., and Tang, Y. (2017). MicroRNA-27a regulates hepatic lipid metabolism and alleviates NAFLD via repressing FAS and SCD1. Sci. Rep. 7:14493.

Google Scholar

Zhang, Y. J., Xiao, Y. Y., Ma, Y., Liang, N. J., Liang, Y. H., Lu, C., et al. (2020). ROS-mediated miR-21-5p regulates the proliferation and apoptosis of Cr(VI)-exposed L02 hepatocytes via targeting PDCD4. Ecotoxicol. Environ. Saf. 191:110160.

Google Scholar

Zhao, S. Q., Wang, R. Q., Tang, Z. L., and Li, K. (2012). Analysis of the blood biochemistry parameters of three Chinese indigenous pig breeds. China Anim. Husband. Vet. Med. 39, 96–100. doi: 10.3969/j.issn.1671-7236.2007.04.022

CrossRef Full Text | Google Scholar