Experimental male SAMP8 (n = 25 animals, pathogen- and virus-free) and SAMR1 mice (n = 25 animals, pathogen- and virus-free) at the age of 3 months were purchased from Beijing WTLH Biotechnology, Co., Ltd., Beijing, China. The animals were kept in separate cages with standard environment (23°C, 50–60% humidity, and 12 h light/dark cycle) in the laboratory animal center and fed ad libitum. The feed weight of the mice was recorded every day, and their body weight was recorded weekly until they reach 8 months. The mice were anesthetized with sodium pentobarbital (50 mg/kg) via intraperitoneal injection (i.p.). Blood was collected through eyeball removal, kept at 4°C temperature for 4 h, and centrifuged at 4000 rpm for 10 min at 4°C for serum collection. The mice were then euthanized through cervical dislocation and dissected to extract their livers. The serum samples were immediately used for blood chemistry test, and the liver tissues were stored in a liquid nitrogen jar at −196°C for RNA sequencing and other experiments.

All animal experiments were conducted in accordance with the “Guide for the care and use of laboratory animals” (National Research Council, 2011) and approved by the Institutional Animal Care and Use Committee of Beijing Normal University (BNU NO. 2020).

Malondialdehyde (MDA) content was determined following the manufacturer’s instructions (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). The liver samples were weighed and homogenized in phosphate buffered saline, and the supernatant was collected for further analysis. The MDA content was determined via thiobarbituric acid (TBA) reaction following Yun’s method (Yun et al., 2020). The results were expressed in nmol per mg of hepatic tissue protein.

The mice were deprived of food and water overnight prior to the test. Serum concentrations of total triglycerides (TG), total cholesterol (TC), high-density lipoprotein cholesterol (HDL-C), and low-density lipoprotein cholesterol (LDL-C) were measured using commercially available assay kits (Nanjing Jiancheng Bioengineering Institute, Nanjing, China) in accordance with the manufacturer’s instructions.

Total RNA was extracted from the liver samples by using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the standard protocol proposed by Rio et al. (2010). The extracted RNA was then subjected to a strict quality control. RNA integrity was assessed by 1% agarose gel electrophoresis and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, United States), RNA purity was determined using a NanoPhotometer R spectrophotometer (IMPLEN, Westlake Village, CA, United States), and RNA concentration was measured using the Qubit RNA Assay Kit in the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States).

RNA sequencing was conducted as previously described (Zhang et al., 2016). Six high-quality cDNA libraries were constructed using the total RNAs isolated from liver samples: three for the SAMP8 mice and three for the SAMR1 mice. In brief, 3 μg of high-quality RNA per sample was used as input material. Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, United States) was used to remove the ribosomal RNA (rRNA), and ethanol precipitation was conducted to clean up and obtain rRNA-free residues. NEBNext^® Ultra^TM II Directional RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, United States) was applied to the rRNA-depleted RNA to generate sequencing libraries. Fragmentation was conducted at 94°C for 15 min in NEBNext First Strand Synthesis Reaction Buffer (5×). Random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH−) were utilized to synthesize the first-strand cDNA. DNA Polymerase I and RNase H were employed to synthesize the second-strand cDNA, which was then purified using 1.8 × Agencourt^® AMPure^® XP Beads (Beckman Coulter, Brea, CA, United States). NEBNext End Prep Enzyme Mix was used to perform end repair for cDNA library at 20°C for 30 min, followed by 65°C for 30 min. NEBNext Adaptor with hairpin loop structure was ligated to the cDNA using Blunt/TA Ligase Master Mix at 20°C for 20 min. The library fragments measuring 150–200 bp were purified with AMPure XP System. The size-selected and adaptor-ligated cDNA was treated with 3 μl of NEBNext USER^® Enzyme at 37°C for 15 min, followed by 5 min at 95°C prior to PCR amplification. PCR was performed with Universal PCR Primer, Index (X) Primer, and Phusion^® High-Fidelity DNA Polymerase. The PCR products were purified using the AMPure XP system, and library quality was assessed by the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, United States) was employed to cluster the index-coded samples. The libraries were forwarded to a sequencing run on the Illumina HiSeq 4000 system at the Novogene Bioinformatics Institute (Beijing, China) to generate 150 bp paired-end reads.

High-quality clean reads were obtained using in-house perl scripts by removing the reads with adaptors, low-quality reads, and reads poly-N >10% from raw data. The GC%, Q30, and Q20 of the clean data were calculated. GRCm38 was used as a mouse reference genome, and the paired-end clean reads were mapped to this reference sequence¹ using TopHat v.2.0.9 (Kim et al., 2013). The mapped reads of each sample were assembled to transcripts by running Cufflinks v2.1.1 software (Trapnell et al., 2010).

The pipeline for lncRNA discovery and identification was prepared as previously described (Zhang et al., 2016). Transcripts with <200 bp in length, less than two exons, and less than three reads coverage were first removed. The remaining transcripts were then compared with known mouse lncRNAs, mRNAs, snoRNA, snRNA, pre-miRNA, tRNA, rRNA, and pseudogenes and then assessed using CNCI, CPC, PFAM, and phyloCSF software (Kong et al., 2007; Mistry et al., 2007; Lin et al., 2011; Sun et al., 2013). Finally, the qualifying lncRNAs were selected and designated as novel or known lncRNAs.

PhyloFit and phastCons from Phast package were used to compute a set of conservation scores for lncRNAs and protein-coding transcripts (Siepel et al., 2005).

Long non-coding RNA and mRNA expression levels in each sample were measured as fragments per kilobase of exon model per million mapped fragments (FPKM) by using the Cuffdiff. Transcripts with p adjust value < 0.05 (q value < 0.05) were considered as statistically significant between SAMP8 and SAMR1.

Long non-coding RNAs may affect the gene expression by playing cis and trans regulating roles (Guil and Esteller, 2012; Yan et al., 2017). Here, only the differentially expressed lncRNAs and mRNAs were used in the prediction to explore the potential function of lncRNAs. The mRNAs within 100 kb upstream and downstream of lncRNAs were selected as cis results. For the trans role, the lncRNAs identify each other based on expression levels, and the absolute values of Pearson’s correlation coefficients >0.95 (|r| > 0.95) were selected.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were applied to the differentially expressed protein-coding genes and target genes of lncRNAs. GOseq R package was used to perform GO analysis. GO terms with p value < 0.05 were recognized as significantly enriched. KOBAS software was used to detect the enriched KEGG pathways. Hypergeometric p value < 0.05 was considered significant.

qPCR was performed on a 7500 Real-time PCR machine (Applied Biosystems, Foster City, CA, United States) by using the SYBR Green qPCR Kit (GenePharma, Shanghai, China) following the manufacturer’s protocol. The 20 μl reaction volume contained 8.2 μl of H₂O, 0.4 μl of each primer, 1 μl of cDNA, and 10 μl of 2 × SYBR Green qPCR Mix. The conditions were as follows: 95°C for 5 min, followed by 40 cycles (95°C for 15 s and 60°C for 30 s). Specific quantitative primers were designed using the primer-BLAST tool on the NCBI website as listed in Table 1. Mouse housekeeping gene β-actin was selected as an internal control. Each experiment was performed in triplicate.

Data were analyzed using Graph pad prism 8 software and represented as mean ± standard error of the mean (SEM). Body weight data were compared by applying two-way repeated measures ANOVA. Two-tailed unpaired Student’s t-test was employed to compare other types of data. Statistical significance was accepted at p value < 0.05.

Figures 1A–C show no remarkable difference in the body weight, liver weight, and food intake between SAMP8 and SAMR1 mice at 8 months of age (p value > 0.05). This finding implied that aging does not alter the body weight, liver weight, and food intake of these mice. Oxidative stress marker MDA is one of the major end products of lipid peroxidation and is traditionally used as an indicator of liver tissue damage (Wang et al., 2019). Figure 1D shows that the levels of liver MDA in 8-month-old SAMP8 mice were significantly increased compared with those in SAMR1 mice (p value < 0.05). In addition, the serum LDL-C levels of SAMP8 mice were dramatically higher than those of SAMR1 mice (p value < 0.05, Figure 1E). No differences in TG, TC, and HDL-C were found between the two groups (p value > 0.05, Figure 1E). These findings revealed that aging is a potential cause of impaired hepatic cholesterol synthesis and metabolism, abnormal function of low-density lipoprotein receptor, and even fatty liver diseases. The livers of the 8-month-old SAMP8 mice initially exhibited impairment on oxidative stress and lipid metabolism, the common clinical characteristics of aging human liver.

A total of 628,530,304 raw reads were generated using the Illumina HiSeq 4000 system. After the low-quality reads were removed, 605,801,688 clean reads were obtained. The high-quality clean reads were mapped to the mouse reference genome sequence by TopHat v.2.0.9. The mapping rates were approximately 87.11 and 84.01% in SAMP8 and SAMR1 mice, respectively. In accordance with the Cufflinks v2.1.1 results, 137,750 transcripts were assembled and used for the following analysis.

After the initial screening, 1,736 known mouse lncRNAs corresponding to 1,399 lncRNA genes and 797 presumed lncRNAs were detected. CNCI, CPC, PFAM, and phyloCSF tools were used to further confirm these 797 presumed lncRNAs. A total of 446 novel lncRNAs (defined as LNC_000XXX) corresponding to 317 lncRNA genes were identified (Figure 2). In addition, 54,275 mRNA transcripts corresponding to 21,959 protein-coding genes were also detected. In summary, 2,182 lncRNAs and 54,275 mRNAs were detected in SAMP8 and SAMR1 mouse livers.

The characteristics of the obtained 2,182 lncRNAs and 54,275 mRNAs were analyzed. Most of the known and novel lncRNAs contained fewer exons than mRNAs (Figure 3A). The average full length and open reading frame (ORF) length of the known and novel lncRNAs were shorter than those of mRNAs (Figures 3B,C). The conservatism of the known and novel lncRNAs was significantly lower than that of mRNAs (Figure 3D). In addition, the known and novel lncRNAs had significantly lower expression than mRNAs (Figure 3E). These results are consistent with previous studies (Singh et al., 2017; Wu et al., 2018).

FPKM was used to estimate the expression of lncRNA and mRNA transcripts. Among the 28 significantly dysregulated lncRNA transcripts (corresponded to 28 lncRNA genes) identified between the two groups (Table 2, q value < 0.05), 12 were upregulated and 16 were downregulated in SAMP8 mice. Among the 210 remarkably dysregulated mRNA transcripts (corresponding to 205 protein-coding genes) identified, 107 were upregulated and 103 were downregulated in SAMP8 mice (Supplementary Table 1, q value < 0.05). Cluster analysis and principal component analysis (PCA) were performed for the 28 differentially expressed lncRNAs (Figures 4A,B) and 210 significantly dysregulated mRNAs (Figures 4C,D). The three replicates of the SAMP8 group were clustered together, and the same was performed on the SAMR1 group.

Cis and trans analyses were performed to predict the underlying relationships between altered lncRNAs and mRNAs. With 100 kb as the cutoff in cis, 2 out of 28 differentially expressed lncRNA genes corresponded to two protein-coding genes (Supplementary Table 2). According to Pearson’s correlation coefficients (|r| > 0.95) in trans, 12 out of 28 differentially expressed lncRNA genes corresponded to 41 protein-coding genes (Supplementary Table 2). The abovementioned target genes were subjected to GO and KEGG analyses. A total of 147 GO terms (Supplementary Table 3, p value < 0.05) and 26 KEGG pathways (Supplementary Table 4, p value < 0.05) were significantly enriched. The top five GO terms were arachidonic acid epoxygenase activity (GO: 0008392), oxidoreductase activity (GO: 0016491), chitin catabolism (GO: 0006032), intracellular membrane-bounded organelle (GO: 0043231), and icosanoid biosynthesis (GO: 0046456). The top five KEGG pathways were retinol metabolism, inflammatory mediator regulation of TRP channels, metabolic pathways, phagosome, and allograft rejection. Several liver aging-associated terms and pathways were also observed, such as oxidoreductase activity (GO: 0016491), lipid catabolism (GO: 0016042), steroid hydroxylase activity (GO: 0008395), retinol metabolism, and metabolic pathways. Therefore, lncRNAs could regulate protein-coding genes associated with liver aging through their cis and trans roles.

Gene Ontology and KEGG surveys were also performed on 205 significantly dysregulated protein-coding genes. A total of 237 GO terms (Supplementary Table 5, p value < 0.05) and 33 KEGG pathways (Supplementary Table 6, p value < 0.05) were significantly enriched. The top five GO terms were oxidoreductase activity (GO: 0016491), arachidonic acid epoxygenase activity (GO: 0008392), heme binding (GO: 0020037), oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen (GO: 0016705), and oxidation–reduction (GO: 0055114). The top five KEGG pathways were retinol metabolism, metabolic pathways, steroid hormone biosynthesis, arachidonic acid metabolism, and chemical carcinogenesis. These results were reflected in liver aging, such as several GO terms (e.g., GO: 0016491, GO: 0055114, and GO: 0004806) and KEGG pathways (e.g., retinol metabolism, metabolic pathways, and fatty acid degradation).

Two limiting conditions were imposed to deeply understand the relationship of lncRNAs with liver aging and related diseases. The first was that the lncRNAs and their target genes must be differentially expressed in the liver tissues of SAMP8 and SAMR1 mice, and the other was the selected pairs (lncRNA–mRNA) must be associated with liver aging and related diseases. A pair matching the above two conditions would be selected. Ces1g reduces hepatic steatosis and counteracts dyslipidemia (Bahitham et al., 2016) and is regulated by LNC_000027, ENSMUST00000144661.2, ENSMUST00000181906.1, and LNC_000204. Bco2 maintains normal hepatic lipid and cholesterol homeostasis (Lim et al., 2018) and is acted upon by LNC_000027 and LNC_000204. Neu1, a gene involved in hepatic glucose and lipid metabolism (Hu Y. et al., 2019), is targeted by ENSMUST00000144661.2, LNC_000366, and LNC_000150. Dhtkd1 possibly modulates mitochondrial biogenesis, enhances energy expenditure to control liver steatosis (Lim et al., 2014), and is acted upon by LNC_000204, LNC_000150, ENSMUST00000144661.2, and ENSMUST00000181906.1. Card6 protects against hepatic steatosis by suppressing Ask1 (Sun et al., 2018) and is regulated by ENSMUST00000180730.3. Gstm3 is an anti-oxidative gene that acts against liver injury (Sun et al., 2017; Zhao et al., 2017) and is controlled by ENSMUST00000181906.1, LNC_000204, and LNC_000150. VLDLR plays an important role in hepatic steatosis (Zarei et al., 2018) and is governed by LNC_000027 and LNC_000426. The proliferation of hepatocellular carcinoma cells is inhibited by the knockdown of BLVRB (Huan et al., 2016), which is regulated by LNC_000204. Sort1 overexpression in mouse liver reduces plasma triglycerides and cholesterol (Strong et al., 2012; Conlon, 2019). Cyp2c44, a member of cytochrome P450 (P450), mediates arachidonic acid metabolism and further regulates inflammation in hepatic tissues (Theken et al., 2011). Sort1 and Cyp2c44 are targeted by LNC_000150. The detailed results are listed in Table 3. The expression of these selected pairs was examined by qPCR to further prove their reliability. As shown in Figure 5, all the mRNA and lncRNA transcripts were detected and exhibited significantly different expression between the SAMP8 and SAMR1 mice. The qPCR results were consistent with the RNA sequencing data. Therefore, these lncRNAs are most likely to participate in liver aging and related diseases.