Abstract
Background/Aims: Intrauterine growth restriction (IUGR) is a risk factor for adult metabolic syndrome, but how this disease is regulated by lncRNAs and circRNAs remains elusive. Methods: Here, we employed adult IUGR and normal pigs as models to evaluate the expression of various global lncRNAs and circRNAs in pig livers using RNA-seq. Results: In total, we obtained 1,162 million raw reads of approximately 104.54 Gb high quality data. After a strict five-step filtering process, 3,368 lncRNAs were identified, including 300 differentially expressed lncRNAs (p < 0.05) in the IUGR group relative to the control group. The cis-regulatory analysis identified target genes that were enriched in specific GO terms and pathways (p < 0.05), including amino acid metabolism, oxidoreductase activity, PPAR signaling pathway, and insulin signaling pathway. These are closely related to the observed phenotypes of increased gluconeogenesis and impaired mitochondrial oxidative phosphorylation in adulthood of the IUGR group. Additionally, we also identified 403 circRNAs, of which 44 were differentially expressed (p < 0.05). Interestingly, our results identified ATF4-miR214-circRNA7964 and TCF7-miR22-3p-circRNA16347 as two competing endogenous networks, which were closely associated with the observed increase in hepatic gluconeogenesis in the IUGR group. Conclusion: Together, this study reveals a multitude of candidate lncRNAs and circRNAs involved in the development of IUGR pigs, which could facilitate further researches on the molecular mechanisms of metabolic syndrome.
Introduction
Intrauterine growth restriction (IUGR) is usually attributed to utero-placental insufficiency. IUGR is a common complication of pregnancy, whichis associated with the future development of metabolic diseases in adulthood, such as type 2 diabetes and cardiovascular disease [1, 2]. It is widely known that the liver is a central organ of metabolic homeostasis, and also a major organ for synthesis, metabolism, storage, and redistribution of carbohydrates, proteins, and lipids [3]. Therefore, the liver increasingly became the primary tissue that was used to study the metabolic syndrome induced by IUGR. For example, fetuses exposed to IUGR have increased fetal hepatic gluconeogenic capacity and reduced hepatic mRNA translation initiation [4]. Peterside et al. revealed that IUGR in rats impaired hepatic mitochondrial oxidative phosphorylation and increased hepatic glucose production by suppressing pyruvate oxidation [5]. Additionally, Yan et al. reported that IUGR in pigs increased fatty acid flux to the liver and reduced lipolysis and fatty acid oxidation [6]. In recent decades, many experimental animal models and epidemiological data have revealed that IUGR is accompanied by the metabolic syndrome. However, there is a lack of evidence in whole-genome differential expression profiles to elucidate the genetic mechanism of metabolic syndrome induced by IUGR, especially in long non-coding RNA and circRNA transcriptomes.
Genome-wide transcriptional studies revealed that large regions of eukaryotic genomes transcribe into non-coding RNAs (ncRNAs), such as microRNA (miRNA) [7], PIWI-interacting RNA (piRNA) [8], small interfering RNA (siRNA) [9], long non-coding RNA (lncRNA) [10], and circular RNAs (circRNAs) [11]. Generally, linear ncRNAs longer than 200 nucleotides are considered as lncRNAs [12], and circRNAs are characterized by a covalently closed-loop structure generated from non-sequential back-spliced of pre-mRNA transcripts [13]. Recently, an increasing number of lncRNAs and circRNAs have been discovered in mammals, including Homo sapiens [14, 15], Mus musculus [16, 17], Bos taurus [18, 19], and Sus scrofa [20, 21]. Accordingly, lncRNAs and circRNAs play crucial roles in multiple biological processes, including transcriptional regulation [22, 23], cell differentiation and development [12, 24], and also in some diseases [25, 26]. However, little is known about the biological functions and molecular mechanisms of lncRNAs and circRNAs in the development of metabolic syndrome in mammals with IUGR.
The pig is an ideal biomedical model for metabolic studies because it’s similar to humans with regard to physiology and metabolism [27]. Therefore, in this study, the pig was used as a model to identify differential expression patterns of whole-genome lncRNAs and circRNAs in the livers between adult IUGR and normal pigs. To our knowledge, this is the first study to elucidate the underlying molecular mechanisms of lncRNAs and circRNAs responding to the metabolic diseases arising from IUGR.
Materials and Methods
Ethics Statement
All animal experimental procedures and sample collections were performed in accordance with the guidelines of the Institutional Animal Care and Use Committee of the College of Animal Science and Technology of Sichuan Agricultural University, Sichuan, China, under permit No. DKY-B20131403 (Ministry of Science and Technology, China, revised in June 2004).
Animal and tissue preparation
A total of 8 pairs of normal and IUGR PIC breed male piglets chosen from 8 sows were used in this study. One IUGR piglet and one normal piglet were selected from each litter. All piglets were living in the same environment and fed with milk formula ad libitum before 28 days of age. The formula composition and nutrient levels of the milk formula fed to the pigs are provided in Supplementary Table S1 (For all supplemental material see www.karger.com/10.1159/000494794/). From 28 days of age until 150 days of age, the pigs were fed twice daily with formulas that meet the National Research Council (NRC 1998) recommendations for different growth phases and were given water ad libitum (Supplementary Table S2). Approximately 24 h before they were euthanized, the feed was removed, but access to water ad libitum was maintained. The pigs were electrically stunned, exsanguinated, scalded, and rinsed. Liver samples were obtained immediately after exsanguination and rapidly frozen in liquid nitrogen for sequencing analyses and phenotype measurements.
Measurement of metabolite phenotypes
The concentrations of hepatic glycogen (No. A043) and triglyceride (No. A110-2) were determined by commercial kits (Nanjing Jiancheng Institute of Bioengineering, Nanjing, Jiangsu, China), according to the manufacturer’s instructions for the automatic biochemical analyzer (Model 7020, Hitachi, Tokyo, Japan). The activity of ATPase (No. A070-2) was also determined using commercial kits (Nanjing Jiancheng Institute of Bioengineering, Nanjing, Jiangsu, China) andcommercial RIA kits (R&D Systems Europe Ltd., Abingdon, UK).
Intravenous glucose tolerance test
According to our previous study, pigs were administered an intravenous glucose tolerance test (i.v.GTT) following an overnight fast at 149 days old. Briefly, dextrose (500 g/L) was infused continuously through ear venipuncture over 6 minutes at an infusion dose of 0.5 g/kg of body weight and an infusion rate of 10 g glucose/minute. Before dextrose infusion, blood samples were obtained at -6, -4, -2, and 0 min, and post-administration blood samples were obtained at 5, 10, 15, 30, 60, 90, and 120 min. Blood glucose was measured using an Ascensia Elite glucometer (Bayer Healthcare Company, Leverkusen, Germany).
Total RNA and DNA extraction
Total RNA was extracted from the liver tissue using TRIzol (Invitrogen, CA, USA). A portion of the total RNA was reverse-transcribed to cDNA for qPCR analysis using a PrimeScriptTM RT Reagent Kit with gDNA Eraser (TAKARA, Dalian, China), which had oligo(dT) and random hexamer primers. Another portion of the total RNA was further purified with an RNeasy column for RNA-Seq (Qiagen, USA), according to the manufacturer’s protocol. The RNA integrity and concentration were assessed using Bio-analyzer 2100 (Agilent Technologies) and NanoDrop (Thermo Technologies), respectively. DNA was isolated from the liver tissue to measure mtDNA copy numbers using the DNeasy Blood & Tissue Kit (Qiagen, USA).
Measurement of mitochondrial DNA (mtDNA) copy numbers
The relative mtDNA copy numbers were determined by qPCR. The ratios of mitochondrial genes (ATP6 and COX2) to nuclear DNA single copy gene (GCG) within the same sample were used to calculate the mtDNA contents (primer sequences are listed in Supplementary Table S3). All reactions were performed in triplicate, and the 2∆Ct method was used to determine the relative mtDNA copy numbers.
Construction of cDNA libraries and high-throughput sequencing
We randomly selected samples of three pigs from each group for deep sequencing analysis. Only RNA samples that had an RNA Integrity Number (RIN) score > 8 were used for sequencing. And then, ribosomal RNA of RNA samples was removed using Ribo-ZeroTM Magnetic Kit (Human/Mouse/Rat) (Epicentre). The strand-specific sequencing libraries were constructed following a previously described protocol [28]. Paired-end sequencing (2 × 90 bp) was performed on an Illumina Hiseq 2000 sequencer. Three biological repeats were used for construction of libraries. All high-throughput sequencing data have been deposited in NCBI’s Gene Expression Omnibus under GEO Series accession numbers GSE81766.
Reads mapping and transcriptome assembly
The raw data of directional paired-end reads were quality-checked with FastQC, and the adapter contamination and low-quality tags in the raw data were discarded. Ribosomal RNA data were also removed from the remaining data by alignment. Thus, all downstream analyses were performed with high-quality data only. Thereafter, the clean reads from the six cDNA libraries were merged and mapped to the Sus scrofa genome sequence (Sscrofa11.1) using the spliced read aligner TopHat (v2.0.14). To construct transcriptomes, the mapped reads were assembled de novo using Cufflinks (v2.1.1), which was also used to calculate FPKM (fragment per kilobase of exon model per million mapped reads) scores for transcripts in each library. All transcripts were required to have exons greater than 1 and longer than 200 bp in length.
Identification of lncRNAs and circRNAs
To annotate the assembled transcripts, we used the Cuffcompare program from the Cufflinks package. The known protein-coding transcripts were identified according to the annotation of Sus scrofa genome sequence (Sscrofa11.1). The remaining unaligned transcripts were used to identify potential lncRNAs. First, we excluded transcripts that shorter than 200 nt in length, containing less than 2 exons, and containing less than three reads. Next, the coding potential of the remaining transcripts was calculated using the Coding Potential Calculator (CPC) and Coding-Non-Coding Index (CNCI). A transcript with a CPC value lower than -1 and a CNCI value lower than 0 was taken to be an lncRNA. To identify circRNAs, the software find_circ was used to extend the anchor sequences, and the back-spliced reads containing at least two supporting reads were considered to be circRNAs [29]. Briefly, the pipeline was used to find potential circRNAs from the unaligned back-spliced junction reads, and all identified circRNAs that were expressed in at least two samples were retained for further analysis. Differentially expressed genes (DEGs), including those of putative lncRNAs and circRNAs, were quantified as fragments per kilobase of exon model per million fragments mapped (FPKM) using the Cuffdiff program from the Cufflinks package. Gene expression differences were evaluated using an adjusted p value < 0.05 and |log2(fold change)| ≥ 1 as the threshold.
Functional enrichment analysis
LncRNAs have been implicated in regulation of the expression of neighboring genes (cis-acting regulation) [30, 31]. Thus, we searched for coding genes 10 kb upstream and downstream of all identified lncRNAs and then predicted their functional roles in cis-regulation. Thereafter, the names of the neighboring genes were used to form a gene list, which wasused for input into DAVID software for GO and KEGG enrichment analysis. In all tests, the P-values were calculated using Benjamini corrected modified Fisher’s exact test. Only P-values < 0.05 were considered significant.
Quantitative RT-PCR
The expression levels of selected miRNAs, mRNAs, lncRNAs and circRNAs were quantified using qRT-PCR. mRNAs, lncRNAs and circRNAs were performed with reverse transcription using oligo (dT) random hexamer primers provided in the PrimeScript RT Master Mix kit (TaKaRa, Dalian, China), and miRNAs were performed with reverse transcription using PrimeScriptTM miRNA RT-PCR Kit (TaKaRa, Dalian, China), following the manufacturer’s recommendations. The qRT-PCR reactions were performed using the SYBR Green Real-time PCR Master Mix (TaKaRa, Dalian, China) on a CF96 Real-Time PCR Detection System (Bio-Rad, Richmond, CA, USA). The reaction volume contained 10 μl SYBR Premix Ex TaqTM II, 1 μl of 10 μM forward and reverse primers, 2 μl template cDNA, and dH2O to a final volume of 20 μl. The PCR primer sequences used are presented in Table S3. ACTB, TBP, and TOP2B genes were used as internal control genes for normalization. All measurements contained negative controls that containing no cDNA template. Each RNA sample was analyzed in triplicate. The 2-ΔΔCt method was used to determine the relative abundance of each transcript, and data were expressed as least square mean ± standard error of the mean (SEM). All primers used are listed in Supplementary Table S3.
Statistical analyses
Data were analyzed using one-way ANOVA to test the homogeneity of variances followed by a Student’s t test in SPSS (21.0 version). Statistical significance was determined as p < 0.05.
Results
Differential phenotypes between IUGR and normal pigs
Consistent with our previous study [32], IUGR pigs had a birth weight 36% lighter than their control littermates (Fig. 1A, p < 0.01). Under the same conditions, the difference of body weight was more robust between IUGR and control pigs over the course of experiment (Fig. 1A). Interestingly, we also found that adult IUGR pigs had reduced hepatic glycogen and increased hepatic triacylglycerol when compared to normal pigs (Fig. 1B-1C, p < 0.05). Additionally, an intravenous glucose tolerance test (i.v.GTT) indicates that IUGR impairs the capacity of glucose tolerance in adult IUGR pigs (Fig. 1D). Moreover, mitochondria, which is essential to energy metabolism, was also linked with IUGR [36]. In the current study, mitochondrial DNA content and ATP activity were both significantly decreased in adult IUGR pigs (p < 0.05), as shown in Fig. 1E and 1F. Overall, these results indicate that IUGR affects the hepatic metabolic pathway of adult pigs.
Phenotypical differences between adult normal and IUGR pigs. (A) The body weights of normal and IUGR pigs (n = 8). The hepatic contents of glycogen (B) and triacylglycerol (C) in adult normal and IUGR pigs (n = 8). (D) Plasma glucose concentrations after an intravenous glucose tolerance test (i.v.GTT) at 149 days old. Time indicates minutes relative to completion of dextrose infusion (n = 8). (E) The copy numbers of hepatic mitochondrial DNA in adult normal and IUGR pigs (n=8). (F) The ATP contents in livers of adult normal and IUGR pigs (n = 8). Data are means ± SEM. Statistical significance was calculated by Student’s t-test. Significant differences levels: * p< 0.05, ** p< 0.01.
Phenotypical differences between adult normal and IUGR pigs. (A) The body weights of normal and IUGR pigs (n = 8). The hepatic contents of glycogen (B) and triacylglycerol (C) in adult normal and IUGR pigs (n = 8). (D) Plasma glucose concentrations after an intravenous glucose tolerance test (i.v.GTT) at 149 days old. Time indicates minutes relative to completion of dextrose infusion (n = 8). (E) The copy numbers of hepatic mitochondrial DNA in adult normal and IUGR pigs (n=8). (F) The ATP contents in livers of adult normal and IUGR pigs (n = 8). Data are means ± SEM. Statistical significance was calculated by Student’s t-test. Significant differences levels: * p< 0.05, ** p< 0.01.
Identification of lncRNAs in pig livers
To identify lncRNAs expressed in IUGR and normal porcine livers, 6 cDNA libraries from two groups were constructed and sequenced using the Illumina HiSeq 2000 platform. A total of 1, 161, 593, 872 raw reads (104.54 Gb) were generated in all libraries (Supplementary Table S4). After discarding adaptor sequences and low-quality reads, approximately 85.84% - 91.41% of all reads were located within exons (Supplementary Fig. S1), and 64.30% - 68.68% of all reads were uniquely mapped to the UCSC pig reference genome (Sus scrofa 11.1). Additionally, we obtained 1.85%-3.21% of back-spliced junction reads for further circRNAs identification (Supplementary Table S4). Subsequently, we developed the rigorous criteria filtering pipeline to discard transcripts that did not have all characteristics of lncRNAs (Fig. 2). A total of 3368 lncRNAs from the pig livers were identified and subjected to further analyses, including 3285 intergenic lncRNAs and 83 anti-sense lncRNAs (Supplementary Table S5).
To explore the differential features between lncRNAs and mRNAs, we estimated the average expression level of the two groups as log10 (FPKM+1). We found that the average lncRNAs expression level was lower than the mRNAs average expression level (Fig. 3A). Protein-coding transcripts had an average length of 2234 bp and 9.6 exons, which was longer than the lncRNA genes that averaged 684 bp in length and 2.5 exons (Fig. 3B-3C). However, the exon size in the lncRNA genes was greater than the exon size in the protein-coding genes (Supplementary Fig. S2). Furthermore, the average length of the predicted ORFs in the lncRNAs was approximately 56.05 bp. In contrast, the average length of mRNA ORFs was 429.45 bp (Fig. 3D).
Comparisons of genomic architecture between mRNAs and lncRNAs. (A) Box plot indicating expression levesl of total mRNAs and lncRNAs by log10 (FPKM + 1). (B) Distribution of transcript lengths of mRNAs and lncRNAs in the livers. The x axis represents the lengths of transcripts, and the y axis represents composing proportions. (C) Distribution of the number of exons in mRNAs and lncRNAs. Single-exon lncRNAs were filtered out from the pig genome. (D) Distribution of the lengths of open reading frames (ORFs) in mRNAs and lncRNAs. The ORF was identified using Estscan.
Comparisons of genomic architecture between mRNAs and lncRNAs. (A) Box plot indicating expression levesl of total mRNAs and lncRNAs by log10 (FPKM + 1). (B) Distribution of transcript lengths of mRNAs and lncRNAs in the livers. The x axis represents the lengths of transcripts, and the y axis represents composing proportions. (C) Distribution of the number of exons in mRNAs and lncRNAs. Single-exon lncRNAs were filtered out from the pig genome. (D) Distribution of the lengths of open reading frames (ORFs) in mRNAs and lncRNAs. The ORF was identified using Estscan.
Identification of differentially expressed lncRNAs between IUGR and normal pigs
To ensure the reproducibility and reliability of RNA-seq libraries in this study, all transcripts in each library were clustered using principal components analysis (PCA). As shown in Fig. 4A, there was high uniformity between samples obtained from two groups. These results confirm the high reliability of lncRNAs identified in this study. Moreover, the expression distribution of lncRNAs from the six libraries was drawn along 20 chromosomes, and they exhibited no obvious location preference (Fig. 4B). These results are consistent with a previous report [33]. The expression levels of lncRNA transcripts were estimated by FPKM using Cuffdiff, which was also used to identify the global lncRNA changes. As shown in Fig. 4C, we have identified 300 lncRNA transcripts that are differentially expressed between IUGR and normal pigs (Supplementary Table S6). Additionally, the heat map of the hierarchical clustering analysis indicated that the differentially expressed lncRNAs were also highly reproducible (Fig. 4C). Of these 300 lncRNAs, 63 lncRNAs were highly expressed in IUGR pigs, and 237 lncRNAs were highly expressed in normal pigs. All these differentially expressed lncRNAs had more than a 2-fold change (Fig. 4D, p < 0.05).
Identification the differentially expressed lncRNAs between normal and IUGR pigs. (A) Principal components analysis (PCA) of FPKM from each sample. (B) Circos plot showing the distribution of lncRNAs in different chromosomes. The different samples read from the inside circle to the outer circle: normal rep1, normal rep2, normal rep3, IUGR rep1, normal rep2, normal rep3. (C) The clustered heat map of 300 differentially expressed lncRNAs for the comparisons between normal and IUGR pigs. (D) The distribution of differentially expressed lncRNAs between normal and IUGR pigs.
Identification the differentially expressed lncRNAs between normal and IUGR pigs. (A) Principal components analysis (PCA) of FPKM from each sample. (B) Circos plot showing the distribution of lncRNAs in different chromosomes. The different samples read from the inside circle to the outer circle: normal rep1, normal rep2, normal rep3, IUGR rep1, normal rep2, normal rep3. (C) The clustered heat map of 300 differentially expressed lncRNAs for the comparisons between normal and IUGR pigs. (D) The distribution of differentially expressed lncRNAs between normal and IUGR pigs.
Enrichment analysis of nearest-neighbor genes of differentially expressed lncRNAs (cis-regulation)
Extensive research has shown that lncRNAs may affect the gene expression of their chromosomal neighborhood 10 kb upstream and downstream in the cis-regulatory [30, 31]. To investigate the possible functions of lncRNAs, we analyzed gene pairs formed by lncRNAs and their neighboring genes. We identified 300 differentially expressed lncRNAs that were transcribed to (< 10 kb) 118 differentially expressed protein-coding genes (Supplementary Table S7). Gene Ontology (GO) analysis of the cis lncRNA target genes was performed to explore their functions. We found 262 GO terms that were significantly enriched (p < 0.05); the top 10 enriched terms included endopeptidase and amino acid metabolism and oxidoreductase activity (Fig. 5 and Supplementary Table S7). Interestingly, we also found amino acid metabolism enrichment in adulthood of IUGR pigs. This agrees with the observed gluconeogenesis phenotypes of reduced hepatic glycogen and impaired glucose tolerance capacity in IUGR pigs (Fig. 1B, 1D). Furthermore, we found genes, including phosphoenolpyruvate carboxykinase-1 (PCK1) and glucose 6-phosphatase (G6PC), which were annotated with the hepatic gluconeogenic related GO term: d-amino-acid oxidase activity (GO: 0003884). Besiedes, we also found a significantly enriched GO term related to OXPHOS: oxidoreductase activity (GO: 0016641). This is consistent with the observed reduction of mitochondrial DNA content and ATP activity in adult IUGR pigs (Fig. 1E-1F). Moreover, we predicted potential targets of lncRNAs in cis-regulatory relationships using pathway analyses. It revealed that the top 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were related to the insulin signaling pathway, PPAR signaling pathway, and pyruvate metabolism (Fig. 5 and Supplementary Table S8).
GO and KEGG annotations for neighboring gene functions of predicated lncRNAs (cis-regulation). The EASE score, which indicated the significance of the comparison, was calculated by Benjamini-corrected modified Fisher’s exact test.
GO and KEGG annotations for neighboring gene functions of predicated lncRNAs (cis-regulation). The EASE score, which indicated the significance of the comparison, was calculated by Benjamini-corrected modified Fisher’s exact test.
Identification of differentially expressed circRNAs in pig livers
Recent researchers have shown increasing interest in the functions of circRNAs, which could negatively regulate the activity of miRNAs by acting as miRNA sponges and competing endogenous RNA (ceRNA) [34, 35]. To explore the mechanism of circRNAs in regulating IUGR development, we identified potential circRNAs from the unaligned back-spliced reads. Here, we identified 403 circRNAs, which showed no noticeable preference for genomic loci (Supplementary Table S9). Additionally, there was no significant difference in exon numbers, and most circRNAs contained two to several exons with an average of 3.6 exons (Supplementary Table S9). This is comparable to the characteristics of circRNAs found in a cattle model [36]. Based on the analysis of circRNA expression, we found that the average expression and expression density of circRNAs were not significantly different between IUGR and normal pigs (Supplementary Fig. S3-S4). Thus, there were only 44 differently expressed circRNAs identified (p < 0.05) (Supplementary Table S10). As shown in the hierarchical clustering analysis of differentially expressed circRNAs (Fig. 6A), the three libraries of each group were clearly assigned to one cluster. The high reproducibility of circRNAs suggests that circRNA formation is regulated through specific pathways and does not merely represent random transcriptional noise. As recent studies have shown, circRNA expression does not always correlate with the expression of the linear transcript from which it is derived [37, 38]. Thus, we only analyzed the ceRNA function of circRNAs and constructed the mRNA-miRNA-circRNA network by bioinformatics prediction.
Identification of differentially expressed circRNAs between normal and IUGR pigs. (A) Hierarchical clustering analyses for differentially expressed circRNAs between normal and IUGR pigs. (B) Bioinformatics predicted the ATF4-miR214-circRNA7964 competing endogenous regulation network, and qPCR was used to detect expression levels of ATF4, miR-214, and circRNA7964 in normal and IUGR pigs. (C) Bioinformatics predicted the TCF7-miR22-3P-circRNA16347 competing endogenous regulation network, and qPCR was used to detect expression levels of TCF7, miR-22-3P, and circRNA16347 in normal and IUGR pigs. Data are means ± SEM. Statistical significance was calculated by Student’s t-test. Significant difference levels: * p< 0.05, ** p< 0.01.
Identification of differentially expressed circRNAs between normal and IUGR pigs. (A) Hierarchical clustering analyses for differentially expressed circRNAs between normal and IUGR pigs. (B) Bioinformatics predicted the ATF4-miR214-circRNA7964 competing endogenous regulation network, and qPCR was used to detect expression levels of ATF4, miR-214, and circRNA7964 in normal and IUGR pigs. (C) Bioinformatics predicted the TCF7-miR22-3P-circRNA16347 competing endogenous regulation network, and qPCR was used to detect expression levels of TCF7, miR-22-3P, and circRNA16347 in normal and IUGR pigs. Data are means ± SEM. Statistical significance was calculated by Student’s t-test. Significant difference levels: * p< 0.05, ** p< 0.01.
Compared to the normal pigs, there were 32 circRNAs that were significantly downregulated in IUGR pigs (Supplementary Table S9). However, the top 10 differentially expressed circRNAs between adult IUGR and normal pigs presented in Table 1 revealed that only one circRNA was significantly downregulated in IUGR pigs. Based on our bioinformatics analysis, we found that two of them, circRNA7964 and circRNA16347, shared binding sites with the seed sequences of miR-214 and miR-22-3p, respectively. Here, we found that circRNA7964 was highly expressed in IUGR pigs. CircRNA7964 is a sponge of miR-214 and decreases its expression in IUGR pigs, resulting in an increase in the expression of its target gene, ATF4 (Fig. 6B). Interestingly, previous studies have shown that ATF4 related to positively regulating hepatic gluconeogenesis [39]. This agrees with the increased hepatic gluconeogenesis capacity observed in IUGR pigs. These results suggest a potential role of the ATF4-miR214-circRNA7964 network in regulating hepatic gluconeogenesis. Besides, Wnt-responsive transcription factor (TCF7) was proved to accelerate the metabolism of gluconeogenesis [40]. As shown in Fig. 6C, the expression pattern of the TCF7-miR22-3P-circRNA16347 network in IUGR and normal pigs was consistent with the observed phenotype of increased hepatic gluconeogenesis in IUGR pigs.
Validation of differentially expressed lncRNAs and circRNAs
To validate the expression levels of lncRNAs and circRNAs in this study, we randomly selected three differentially expressed lncRNAs and circRNAs for detection of their expression patterns in IUGR and normal pigs by qPCR. As shown in Fig. 7A-7B, the qPCR results confirmed that the expression patterns of the three lncRNAs and circRNAs were consistent with their expression levels calculated from the RNA-seq data. In addition, we designed convergent primers and discrete primers to validate the bioinformatics methods used to identify circRNAs. As shown in Fig. 7C, only the discrete primers could amplify PCR products from circRNAs transcripts. For example, the junction regions (381bp) of circRNA4816 was amplified using discrete qPCR primers, but no PCR products were amplified from linear transcripts, such as ACTB. Taken together, our results show that our RNA-seq data and the bioinformatics analysis pipeline provided reliable information about the relative abundance of lncRNAs and circRNAs in vivo.
Validation of RNA-seq data by qPCR. (A) The expression patterns of 3 randomly selected lncRNAs in normal and IUGR pigs validated by qPCR. (B) The expression patterns of 3 randomly selected circRNAs in normal and IUGR pigs validated by qPCR. (C) The schematic view illustrating the design of primers for circRNAs used in qPCR.
Validation of RNA-seq data by qPCR. (A) The expression patterns of 3 randomly selected lncRNAs in normal and IUGR pigs validated by qPCR. (B) The expression patterns of 3 randomly selected circRNAs in normal and IUGR pigs validated by qPCR. (C) The schematic view illustrating the design of primers for circRNAs used in qPCR.
Discussion
Phenotypes results demonstrated that the body weight of IUGR pigs higher than normal pigs over the course of experiment. This agrees with previous findings that adult IUGR pigs had lower body weights than their normal counterparts, which were accompanied by reduced daily food intakes [6]. This results also suggests that the influence of metabolic levels in neonates with IUGR will persistent into the adulthhood [41]. Previous studies have reported that IUGR mammals have increased hepatic gluconeogenic capacity, reduced hepatic lipolysis, and reduced fatty acid oxidation [4, 42]. In a rat model, intravenous glucose tolerance test (i.v.GTT) indicates that IUGR will impairs its capacity of glucose tolerance [43]. Previous studies have shown that IUGR impairs hepatic mitochondrial oxidative phosphorylation and influences hepatic energy homeostasis [44]. Therefore, combining previous studies and our findings in this work, we confirm that pig can also act as a typical IUGR modle.
In this study, a total of 3368 lncRNAs were dientified, which have differential features with mRNAs, such as transcript length, exon number, and length of ORF. A previous study reported that the average lncRNAs expression level was lower than the mRNAs average expression level in pig adipose tissue [45]. Additionally, several studies have shown that lncRNAs were shorter in length and had fewer exons than protein coding transcripts [46, 47]. In a goat model , the exon size in the lncRNA genes was greater than the exon size in the protein-coding genes [48]. Interestingly, we found that both the lengths and exon numbers of lncRNAs of pigs were different from other mammals, including humans [49], mice [50], and goats [48]. These results suggest that these lncRNAs identified in the present study was consistent with previous reports on the typical characterization of lncRNAs [51].
In the cis prediction, we searched for coding genes 10-kb upstream and downstream of all the identified lncRNAs. GO and KEGG analyses of the neighboring protein-coding genes revealed that major enriched pathways were associated with endopeptidase, amino acid metabolism, oxidoreductase activity, insulin signaling pathway, PPAR signaling pathway, and pyruvate metabolism. These results indicate the possible role of lncRNAs in transcriptional regulation of gene expression. Data from several previous studies suggest that IUGR increases fetal hepatic gluconeogenic capacity [4, 52]. One of the major sources of carbon for hepatic gluconeogenesis is from amino acids, including alanine, serine, and threonine [53, 54]. For example, PCK1 and G6PC genes were annotated with the hepatic gluconeogenic related GO term: d-amino-acid oxidase activity (GO: 0003884). Previous study have found PCK1 and G6PC were the predominant genes that encode enzymes involved in the regulation of gluconeogenesis [55]. These results suggest that hepatic gluconeogenesis may be regulated by the influence of lncRNAs on these neighboring protein-coding genes. Recently, hepatic mitochondrial dysfunction was linked to the mechanism of IUGR development. Zhang et al. (2017) reported that IUGR impairs hepatic mitochondrial biogenesis and the efficiency of oxidative phosphorylation (OXPHOS) in suckling piglets [44]. This results in the excess production of mitochondrial superoxide radicals and generates poor hepatic antioxidant defense systems [32]. Insulin signaling is the major pathway that controls glucose homeostasis by coordinating important metabolic processes, including glucose uptake, glycolysis, glucose oxidation, and glycogen synthesis in nearly all body tissues [56]. Previous studies have shown IUGR increases the mRNA levels of hepatic gluconeogenic genes and reduces the mRNA levels of insulin signaling genes in the liver [42]. In the current study, we found genes, including PCK1, GYS2, and GCK, which were indicated in the insulin signaling pathway. All of these genes have been shown to regulate the aberrant glucose metabolism in IUGR [57-59]. The PPAR signaling pathway is involved in the regulation of the AMPK-SIRT1-PGC1α-PPAR signaling cascade, which is responsible for mitochondrial biogenesis and oxidative phosphorylation [60]. Additionally, pyruvate serves as the substrate for gluconeogenesis, which agrees with the high hepatic gluconeogenesis capacity observed in IUGR pigs. Therefore, these findings indicate that lncRNAs play an important role in the cis-regulation of metabolic diseases arising from IUGR. However, even the bioinformatics analysis of lncRNAs in IUGR was highly consistent with the phenotypes in this study, genetic mechanism experiments are needed in the future to verify the regulatory networks.
In the circRNAs prediction, a total of 403 circRNAs were identified, including 44 differentially expressed circRNAs between IUGR and normal pigs. Based on ceRNA function of circRNAs, we found two networks of ATF4-miR214-circRNA7964 and TCF7-miR22-3P-circRNA16347 maybe involve in the biological process of hepatic gluconeogenesis. Previous studies have shown that miR-214 suppresses hepatic gluconeogenesis by targeting Activating Transcriptional Factor 4 (ATF4) [39]. Additionally, Kaur et al. (2015) support the theory that elevated hepatic miR-22-3p expression decreases the expression of enzymes in the gluconeogenic pathway and impairs the gluconeogenesis capacity by silencing the Wnt-responsive transcription factor TCF7 [40]. Interestingly, the expression pattern of the two networks in IUGR and normal pigs was consistent with the observed phenotype of increased hepatic gluconeogenesis in IUGR pigs. Overall, these results indicate that circRNAs may be taken as a key factor in regulating hepatic gluconeogenesis in IUGR pigs by acting as competing endogenous RNAs. Although the two circRNAs’ sponge functions were based on the bioinformatics analysis and agreed with their expression patterns in this study, the definite regulatory networks should be proved by genetic mechanism experiments in the furture.
Conclusion
In summary, we first conducted a genome-wide diversity analysis of lncRNAs and circRNAs expression between adult normal and IUGR pigs. We identified various differentially expressed lncRNAs and circRNAs that are potentially associated with amino acid metabolism, oxidoreductase activity, PPAR signaling pathway, and insulin signaling pathway. Interestingly, these enriched pathways were consistent with the observed phenotypes of metabolic syndrome in IUGR pigs, including increased gluconeogenic activities, reduced mitochondrial biogenesis, and OXPHOX. We believe this study will contribute to further development of novel diagnostic and therapeutic strategies for metabolic syndrome.
Acknowledgements
The study was supported by the Sichuan Science & Tech Support Program (Pig Genetic Resources Exploitation and Utilization for the 13th Five-Year-Project, No. 16ZC2838) the Chinese National Science & Tech Support Program (No. 2015BAD03B01), the earmarked fund for China Agriculture Research System (No. CARS-36-05B). This work was also partially supported by the China Scholarship Council (CSC) under the PhD visiting student exchange program.
Disclosure Statement
The authors declare that they do not have anything to disclose regarding funding or conflicts of interests with respect to this manuscript.