Comprehensive transcriptional profiling of aging porcine liver

Animals

In this experiment, four Yanan sows were used to investigate the porcine liver aging process. For this study, two young sows (180-days-old) and two old sows (8-years-old) were selected. Moreover, no direct and collateral blood relationship existed among the four pigs within the last three generations of these experimental animals. All experiments involving animals were conducted according to the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004), and approved by the Institutional Animal Care and Use Committee in College of Animal Science and Technology, Sichuan Agricultural University (Sichuan, China) under permit No. SKY-B20160201.

Feeding and slaughtering standards of the pigs were performed as described in previous literature (Chen et al., 2018). As in previous research, piglet weaning was executed at the age of 28 ± 1 days. After weaning, an initial diet consisting of 3.40 Mcal kg⁻¹ metabolizable energy having 20.0% crude protein (11.5 g kg⁻¹ lysine) was provided to piglets, which continued from day 30 to 60. The pigs were shifted to a different diet, composed of 14.0 MJ kg⁻¹ of metabolizable energy and 18.0% crude protein (9.0 g kg⁻¹ lysine) from the 61st to the 120th day. On the 121st day, pigs were subjected to a feed constituting 16.0% crude protein (8.0 g kg⁻¹ lysine) and 13.5 MJ kg⁻¹ of metabolizable energy. Similar conditions were provided to pigs with the permission to freely access food and water. Pigs were not allowed to access food one night before slaughtering and had a resting period of 2 h after transportation. To exsanguinate and relieve the pain, a sudden electric shock of 90 V and 50 Hz for 10 s was administered (Chen et al., 2018).

Sample preparation

As previous study described, all samples were investigated following the rules and regulations for the handling and establishment of research animals by the Ministry of Science and Technology of China (Chen et al., 2018). Porcine liver tissue was harvested from each pig and frozen in liquid nitrogen. After that, samples were kept in a freezer at −80 °C until RNA extraction. The TRIzol Reagent (Invitrogen, CA, USA) was used to extract the total RNA, which was subsequently treated with DNase and purified by using the RNeasy Mini Kit (Qiagen, CA, USA). The quality and concentration of RNA were determined by an Agilent Bioanalyzer 2100 system.

RNA sequencing

RNA samples were prepared by using about 5 µg RNA per sample. According to the manufacturer’s information, the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs Ipswich, MA, USA) was used to obtain sequencing libraries. The Agilent Bioanalyzer 2100 system was utilized to analyze the library quality. After clusters had been generated, RNA sequencing was performed using Illumina HiSeq 4000 and 150 bp paired-end reads were obtained.

Small RNA sequencing

A small RNA library was prepared by taking a total of 5 µg of RNA per sample as input material. To generate sequencing libraries, NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (New England Biolabs Ipswich, MA, USA) was used, according to the manufacturer’s information. To attribute sequences to each sample, index codes were supplemented. Here, the Agilent Bioanalyzer 2100 system was also utilized to test the quality of the library. According to the information provided by the manufacturer, a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina, Inc, CA, USA) was used to cluster the index-coded samples. After the clusters had been generated, the Illumina MiSeq platform was used to sequence the obtained library preparations, and 50 bp single-end reads were obtained.

Statistical analysis

Clean reads of RNA sequencing were obtained from raw data after eliminating the adapter, poly-N, and low-quality reads. The clean reads were aligned to the reference genome database (Ensemble Susscrofa 10.2) using TopHat2 (v2.0.14) and default parameters were used (Kim et al., 2013). Mapped reads per sample were assembled by using the StringTie software (Pertea et al., 2015), which at least presented in one of the two replicates. Mapped transcripts blasted (e-value = 1e−10) to Ensemble were directly described as known mRNA or lncRNA. To estimate the transcripts per million (TPMs) of both lncRNAs and mRNA, Salmon (v0.6.0) was utilized (Patro, Duggal & Kingsford, 2015). To examine the transcript coding potential, Coding Potential Calculator (CPC) (0.9) (Kong et al., 2007) and Pfam Scan (v1.5) (Punta et al., 2012) were utilized. Coding potential transcripts predicted by one or both tools were filtered out, and non-coding potential transcripts were identified as the candidate set of novel lncRNAs. The circRNAs were predicted by CIRCexplorer2 (junction reads ≥ 2) (Zhang et al., 2016).

The miRBase21 was used as reference and the software mirdeep2 (Friedländer et al., 2011) was utilized to identify the known miRNA and to predict novel miRNA. The miRanda (v3.3a) software (Betel et al., 2008) was used to predict the target genes of miRNAs and default parameters and cutoffs (Score S ≥ 140 and Energy E ≤−20.0) were used. Later, the expression levels of miRNAs were also measured by TPMs.

The mRNAs, lncRNAs, miRNAs, and circRNAs were identified to be differentially expressed by using the edgeR package (Robinson, McCarthy & Smyth, 2010) in the R programming environment. Moreover, in the “classic analysis” section, the edgeR user guidelines were adopted in detail. A q-value <0.05 and fold change >2 were used to determine differentially expressed mRNAs, lncRNAs, and miRNAs in the old group vs. young group. For differentially expressed circRNAs, the standard p-value <0.05 and a fold change >2 was used. Here, the fold change of mRNAs, lncRNAs, miRNAs, and circRNAs was the log2 transformed ratio of the average expression between two groups, and was estimated by edgeR.

The GOseq R package was used to perform the Gene Ontology (GO) enrichment analysis, and GO terms with a corrected p-value <0.05 were considered to be significantly enriched by the genes. Principal Component Analysis (PCA) was performed by the prcomp function in R. All the Pearson’s correlation coefficient analysis in our research used the cor function in R programming.

CircRNA-miRNA co-expression network

The correlation analysis was the base for the construction of a circRNA-miRNA co-expression network between the differentially expressed circRNA and miRNAs. The Pearson’s correlation coefficient (r) was used for the expression analysis of differentially expressed circRNAs and miRNAs. The r >0.8 and a p-value <0.05 were considered relevant for network construction between a circRNA and a miRNA.

Quantitative PCR validation

The oligo (dT) and random 6-mer primers given in the PrimeScript RT Master Mix kit (TaKaRa, Shiga, Japan) were utilized for cDNA synthesis. The CFX96 Real-Time PCR detection system (Bio-Rad Laboratories, Inc., Hercules, CA, USA) was used to conduct q-PCR using the SYBR Premix Ex Taq kit (Takara, Shiga, Japan). Three biological replicates were taken to perform the q-PCR analysis. In this assay, three endogenous control genes (porcine GAPDH, ACTB, and U6 snRNA) were used. The expression levels of objective mRNAs, miRNAs, lncRNAs, and circRNAs were determined by the 2^−ΔΔCt method. The Pearson’s correlation coefficient of the log2 fold change values between RNA-seq and q-PCR was used to validate the reliability of the obtained RNA-seq data.

Transcriptome profile of the aging liver

To determine the changes in transcriptome level due to aging in the liver, four porcine liver tissues were obtained from two age groups. Here, RNA sequencing (RNA-seq) was used to explore transcriptome-wide profiling including mRNA, miRNA, lncRNA, and circRNA. An estimate ∼71.66 million clean reads were achieved per RNA-seq library. An alignment of ∼73.40–79.40% of clean reads were found with the porcine reference genome (Ensemble Susscrofa 10.2) (Table S1). Moreover, 8.00–9.43 million clean reads were gained per small RNA-seq library. 74.39–80.24% alignment of clean reads were mapped reads (Table S2). Among these four samples (including two young livers (YL) and two old livers (OL)), a total of 21,607 mRNAs, 4028 lncRNAs, 888 miRNAs, and 6366 circRNAs were detected (Fig. 1, Data S1). Out of 888 miRNAs, 336 miRNAs were identified as novel miRNAs (Fig. 1E). A total of 84.5%, 64.4%, and 95.2% of these identified mRNAs, lncRNAs, and miRNAs, respectively, were expressed in both young and old groups; however, only circRNAs had an expression of 16.9%. Moreover, the number of only expressed mRNAs, lncRNAs, miRNAs, and circRNAs in OL exceeded those only expressed in YL (Fig. 1, Data S2).

Unsupervised Euclidean matrix plots for expressed mRNA and lncRNA showed that both the old individuals and younger individuals were grouped separately (Fig. 2A), indicating an age-dependent expression pattern. Moreover, a closer distance was found between old individuals compared to that of younger individuals. PCA was used for the identified mRNA and lncRNA transcripts, and the results showed that the young group had a high degree of variance between each biological replicate, while old individuals clustered (Fig. 2B). Furthermore, four samples were investigated using unsupervised Euclidean matrix plots based on miRNAs (Fig. 2C). The result also indicated age-dependent expression patterns. However, the result of PCA for miRNA has shown that there was a higher degree of variance between each biological replicate in the old group than in the young group (Fig. 2D), which contrasted with the PCA result based on mRNA and lncRNA.

Differentially expressed mRNAs during liver aging

Form our results, though most genes had been shared for both ages during liver development, we also detected a proportion of age-dependent differentially expressed genes (DEGs). First, DEGs were screened in OL vs. YL (OYL), which identified 273 DEGs, of which 238 were up-regulated and 35 were down-regulated (Fig. 3A, Data S3). To explore the biological functions of these age-related genes, GO analysis was performed. Compared to the young group, the results showed that the up-regulated genes were significantly enriched for immune response, the regulation of cell differentiation, stress response, both transport and storage of lipids, such as defense response to virus (GO:0051607), and ISG15-protein conjugation (GO:0032020); positive regulation was found for cell differentiation (GO:0045597), response to biotic stimulus (GO:0009607), positive regulation of lipid storage (GO:0010884), and chylomicron remnant clearance (GO:0034382) (Fig. 3C, Data S4). For example, the gene nuclear receptor subfamily 1 group H member 2 (NR1H2, ENSSSCG00000003211) was identified to be up-regulated (log FC = 2.026481405) (Data S3) during liver aging, and it was GO enriched for the positive regulation of lipid storage (Data S4). The NR1H2 encoding liver X receptor (LXR)-β was found to be universally expressed (Kalaany & Mangelsdorf, 2006), and the variation within the NR1H2 gene may facilitate the development of type 2 diabetes (Ketterer et al., 2011). Moreover, the genetic variability at NR1H2 may contribute to an increased risk of Alzheimer’s disease (Adighibe et al., 2006), and the preeclampsia was associated with a polymorphism in NR1H2 (Mouzat et al., 2011).

The down-regulated genes were mainly related to diphosphatase activity, signaling pathway, and metabolic processes, such as phosphodiesterase I activity (GO:0004528), nucleotide diphosphatase activity (GO:0004551), regulation of nodal signaling pathway (GO:1900107), ethanol metabolic process (GO:0006067), and lipid metabolic processes (GO:00066295) (Fig. 3C, Data S4). The ELOVL fatty acid elongase 2 gene (ELOVL2, ENSSSCG00000001045) was found to be closely related to lipid metabolic processes (Data S4), and was found to be significantly down-regulated (log FC = −3.306051731) (Data S3). A previous study has reported ELOVL2 as a highly expressed gene in the liver that was closely related to the elongation process of PUFAs with 22 carbons to produce 24-carbon precursors for DHA and DPAn-6 formation, indicating that it was essential for lipid homeostasis (Pauter et al., 2014). In addition, in testis, the ELOVL2 can also control the level of n-6 28:5 and 30:5 fatty acids, which had been identified as indispensable for normal sperm formation and fertility in male mice (Zadravec et al., 2011).

Differentially expressed noncoding RNA during liver aging

To detect the expression profiling of age-related lncRNA, the identified lncRNAs were screened out during porcine liver aging. Here, 19 lncRNAs were identified as differentially expressed lncRNAs (DELs) in the aging of porcine liver. Seventeen out of these 19 lncRNAs presented up-regulation, and two exhibited down-regulation (Fig. 3B, Data S3). Additionally, GO analysis was performed to investigate the target genes of DELs bearing biological functions. The results showed that the target genes of up-regulated lncRNAs were significantly enriched for the immune response, signaling pathway, and protein transport, such as positive regulation of type I interferon production (GO:0032481), defense response to virus (GO:0051607), positive regulation of NF-kappaB transcription factor activity (GO:0051092), response to virus (GO:0009615), positive regulation of RIG-I signaling pathway (GO:1900246), epoxygenase P450 pathway (GO:0019373), and SMAD protein import into nuclei (GO:0007184) (Fig. 3D, Data S5). Moreover, the target genes of down-regulated lncRNAs were significantly related to oxidation–reduction processes, metabolic processes, diphosphatase activity, iron ion and heme binding, and transport processes, such as oxidation–reduction process (GO:0055114), oxidoreductase activity (GO:0016491), metabolic process (GO:0008152), phosphodiesterase I activity (GO:0004528), nucleotide diphosphatase activity (GO:0004551), iron ion binding (GO:0005506), heme binding (GO:0020037), lipid transport (GO:0006869), and both bile acid and bile salt transport (GO:0015721) (Fig. 3D, Data S5).

To investigate the co-regulated relationships of lncRNA-mRNA during liver aging, the DELs and their predicted target genes were screened, and three significantly up-regulated lncRNAs (MSTRG.58269, MSTRG.85782, and MSTRG.177174; logFC >9.0) as well as one down-regulated lncRNA (MSTRG.205482; logFC <-2.0) were found. Several of their target genes were differentially expressed during liver aging (Table 1). Here, MSTRG.58269 was found to have biggest change in expression (logFC = 9.880537444), and its target genes DMBX1 (enriched for DNA binding) and GUCA1A (enriched for Calcium sensitive guanylate cyclase activator activity) were both up-regulated in response to liver aging. Furthermore, most of the target genes of MSTRG.85782 were up-regulated during liver aging, and these genes were significantly enriched for immune response, such as IRG6 and ISG15 (enriched for the defense response to virus), S1PR3 (regulation of interleukin-1 beta production), and IRG6 (negative regulation of viral genome replication). In addition, the target genes of MSTRG.205482 were down-regulated in response to liver aging, and were significantly enriched for metabolic process, such as LIPK (related to lipid metabolic process), SULT1B1 (related to flavonoid metabolic process), and TLR4 (related to lipopolysaccharide receptor activity) (Table 1).

A total of 26 differentially expressed miRNAs (DEMs) were identified in OYL, of which 14 were up-regulated and 12 were down-regulated (Fig. S1A, Data S3). Moreover, 103 differentially expressed circRNAs (DECs)were identified in OYL, of which 81 were up-regulated and 22 were down-regulated (Fig. S1B, Data S3). Here, the miR-206 (logFC = 5.300862103) was identified to be significantly up-regulated in response to liver aging (Data S3). A previous study has shown that miR-206 can regulate the glucokinase activity that knocked out the miR-206 and produced the elevated glycogen content in the mice liver (Manjula et al., 2016). Moreover, it was found that miR-206 represses the LXR-α activity and expression in hepatic cells (Vinod et al., 2014).

The circRNA as competitive endogenous RNA (ceRNA) regulate the function of miRNA acting as a ‘sponge’, thus it could regulate the expression of target genes. This indicates that circRNAs and their target miRNAs may exert a co-expression during liver aging (Li et al., 2015; Sanger et al., 1976). Here, the Pearson correlation coefficien was analyzed between circRNA and miRNA based on the expression levels of DECs and DEMs, and a co-expression network map was constructed (Dou et al., 2016). The network map included 5 miRNAs, 13 circRNAs, and 13 relationships (Fig. 4, Data S6). In this network, 11 circRNAs associated with three miRNAs were found to be synchronously up-regulated. However, two circRNAs associated with two miRNAs were found to be synchronously down-regulated. In addition, seven circRNAs with GL893173.1_33471 (novel miRNA) presented most relationships, while circRNA003190 and circRNA017214 showed single relationship with miR-299 and miR-210, respectively (Fig. 4, Data S6).

The expression of genes validated by q-PCR

q-PCR was performed to verify the RNA-seq results. A certain number of genes were chosen that were found to play a significant role in specific diseases or in vital biology functions. Moreover, several DEGs were randomly selected to conduct q-PCR to determine the expression levels of these gene in the four samples (two 180-days pigs and two 8-years pigs). Selected genes were composed of two mRNAs (ELOVL2 and NR1H2), two lncRNAs (MSTRG.6003 and MSTRG.205482), two miRNAs (miR-210 and miR-299), and two circRNAs (circRNA005502 and circRNA015222) (Fig. S2). The results showed that the expression patterns were highly consistent between the two methods (r = 0.887581, p = 0.003259) (Fig. S2).In this study, the q-PCR results affirmed the high reproducibility and reliability of the gene expression profiling from RNA-seq results. The primer pairs used here are listed in Table S3.