Changes of liver transcriptome profiles following oxidative stress in streptozotocin-induced diabetes in mice

Animals

Specific pathogen-free male C57BL/6 mice weighing 20–22 g were purchased from the Organ transplantation center, Tongji hospital affiliated with Tongji Medical College, Huazhong University of Science and Technology. The Zhengzhou University Animal Care and Use Committee approved all animal experiments (the approval number 2017051805), which were performed in accordance with ‘Animal Research: Reporting In Vivo Experiments’ (ARRIVE) guidelines. After 12 h of fasting, mice received one intra-peritoneal injection of 130mg/kg streptozotocin (STZ, Sigma, St. Louis, MO, USA) solution in 0.05 M citrate buffer (pH 4.5) to induce diabetes (DM, n = 20) (Deeds et al., 2011; Huang Xin & Cao, 2007; Mallek et al., 2018). Normal Control group were injected with citrate buffer (NC, n = 10). Blood glucose (BG) was measured to confirm diabetes, which was defined as glycemia higher than 16.7 mmol/L.

After injection, the mice continued to receive a high-fat diet for another 2 weeks. During the 2 weeks, three animals in the DM group died. Blood glucose level of mice was tested from the tip of the tail. Four weeks after the injection, the mice were euthanized by intraperitoneal injection of 250 mg/kg body weight pentobarbital (Sigma, P3761, under sterile conditions) to harvest liver samples, and blood was collected from the orbital vein to measure serum biochemical markers.

Biochemical marker measurement

Then fast plasma glucose and insulin level of diabetic mice (n = 17) and normal control mice (n = 10) were tested and compared to each other. Serum glucose, insulin levels and total SOD activity were tested by the Automatic biochemical analyzer (cobas 8000 series) using Roche regents according to the manufacture instruction. MDA formed from the breakdown of polyunsaturated fatty acids serves as a convenient index for the determination of the extent of peroxidation reaction. MDA, a product of lipid peroxidation, reacts with thiobarbituric acid to give a pink-colored product, having a maximum absorption at 535 nm (Nair & Nair, 2017). 8-iso-PGF2α was determined with a competitive enzyme-linked immunosorbent assay (ELISA) (Stressgen Biotechnologies Inc., San Diego, CA, USA). HbA1c was detected by High Performance Liquid Chromatography with Borate Affinity Chromatography.

Total RNA extraction and purification

Total RNA from liver tissues of normal and diabetic mice was isolated using the NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ispawich, USA) and quantified using Agilent 2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA). 3 µg of total RNA was used for sequencing preparation using NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ispawich, USA) kit along with Ribo-Zero Gold rRNA (Illumina Inc., CA, USA) to remove rRNA according to the previous study (Zhang et al., 2014). The resulting libraries were sequenced on a HiSeq 2000 (Illumina Inc., CA, USA) instrument that generated paired-end reads of 100 nucleotides.

Illumina HiSeq2000 analysis

RNA extracted from the liver tissues of three control mice were pooled together for sequencing. The sequencing reads were obtained from control pools and STZ-induced diabetic mice (n = 3). Raw sequencing reads were further processed with Perl scripts to exclude the adaptor-polluted reads, low-quality reads and reads with the number of N bases accounting for more than 5%. Q30 statistics was performed to test the data quantity and quality. And Clean Data were mapped to the reference genome (http://www.ensembl.org/index.html) using HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) (Kim, Langmead & Salzberg, 2015). The liver transcriptome was reconstructed from all of the RNA-seq datasets using StringTie 1.3.2.d (http://ccb.jhu.edu/software/stringtie/). DESeq (http://www.bioconductor.org/packages/release/bioc/html/DESeq. html) was used for differential expression analysis between diabetic and normal mice liver transcriptomes. Differentially expressed genes were identified based on threshold changes of ≥2-fold or ≤-2-fold and q values ≤ 0.05. The data were normalized and hierarchically clustered with R software 3.1.1. The potential function of the differentially expressed genes was analyzed by Gene ontology and Pathway analysis. The enriched genes in Kyoto Encyclopedia of Genes and Genomes (KEGG) were calculated by hypergeometric distribution.

Quantitative real-time PCR (qRT-PCR) and statistical analysis

Total RNA was extracted using the RNeasy kit (Qiagen, Inc., Valencia, CA, USA) according to the manufacturer’s instructions, and 2 µg purified RNA was reverse transcribed into cDNA (37 °C for 15 min, followed by 85 °C for 5 s using RT kit; Fermentas; Thermo Fisher Scientific, Inc.). Primers for qRT-PCR were designed based on the sequences from ensembl (http://asia.ensembl.org/index.html). qRT-PCR was performed using the ABI 7500 Real-Time PCR using a QuantiTect SYBR Green PCR kit (Qiagen, Inc.). The qRT-PCR cycle was pre-denaturation at 95 °C for 3 min, followed by 35 cycles of denaturation at 95 °C for 5 s and annealing at 60 °C for 30 s, and a final analysis from 60−95 °C. qRT-PCR results were quantified using the 2^−ΔΔct method. β-actin was chosen as a reference gene. All the gene expression levels were normalized to β-actin measured in parallel.

qRT-PCR assays were performed in triplicate and the data represented the means of three experiments. All data were represented as mean  ± standard deviation. Comparison between groups was performed using the independent sample Student t test with P < 0.05 as the criterion for statistical significance. All analyses were done using SPSS statistics (version 17.0) and GraphPad Prism (version 7).

Biochemical parameters of two group mice after 4 weeks of injection

Four weeks after the injection, control mice were weighing 24.02  ± 4.35 g with blood glucose levels 5.4  ± 0.82 mmol/L, and diabetic mice were weighing 18.02  ± 3.86 g with blood glucose levels 17.4  ± 3.28 mmol/L (Fig. 1). Four weeks after STZ injection, the distribution of biochemical parameters in serum of control and diabetic groups was shown in Table 1. The levels of serum Fasting insulin, HbA1c, MDA and 8-iso-PGF2α were significant higher in DM groups than NC group, while SOD activity decreased significantly in DM groups. Fast insulin level and Body weight were significant lower in DM groups compared to NC group.

RNA sequence data generation and quality control

We observed that administration of streptozotocin caused a significant increase in plasmatic glucose and a decrease in insulin levels. Whole transcriptome RNA sequencing of liver tissue was performed to identify differentially expressed RNAs related to the OS. RNA extracted from the liver tissues of three control mice were pooled together as control group, while the RNA samples from liver tissues of three diabetic mice were sequenced separately as diabetic group. We obtained a total of 6.0746 × 10⁸ raw reads (Table S1, Fig. S1A). The raw data is available at NCBI (accession number PRJNA562053). The Q30 Bases Rate was more than 95% by Q30 statistics (Fig. S1B). We found that a large fraction (median percentage, 68.86%) of the sequence was overlapped by exon regions and that only a small fraction (median percentage, 3.305%) was mapped to the intergenic region (Fig. S1C). Interestingly, long intergenic non-coding RNAs (lincRNA) always located in these areas.

Differential expression analysis of liver transcriptomes

The expression of genes was quantified as Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values (Figs. S1D and S1E). The distribution of the gene expression pattern was similar between the diabetic and normal control mice, only a small fraction of genes were differentially expressed. We identified a total of 2376 novel lncRNAs (Figs. 2A, 2B) and 1326 differentially expressed genes (Table S2). Of which there were 287 up regulated mRNAs and 623 down regulated mRNAs, 161 up regulated lncRNAs and 255 down regulated lncRNAs in the STZ-induced diabetic mice compared to the normal control mice (Table S2). The average percentages of SNP variations in the control and DM group were 93% and 83% (Fig. 2C). The alternative splice statistics showed that the splice occurred mainly in Transcription Start Site (TSS) and Transcription Terminal Site (TTS) (Fig. 2D).

To understand the biological pathways and functions altered in STZ-induced diabetic mouse liver, gene ontology and pathway enrichment analysis were utilized to annotate the target genes. GO analyses found that the dysregulated lncRNAs associated with diabetes mellitus were associated with regulation of apoptotic signaling pathway, negative regulation of transcription from RNA polymerase II promoter, fatty acid catabolic and oxidation process, protein modification and localization (Fig. 3).

KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.kegg.jp/) pathway analysis showed the dysregulated lncRNAs mainly involved in four categories (Fig. 4, Table 2). The first category is inflammation, such as hepatitis B, epstein-Barr virus infection, protein processing in endoplasmic reticulum, lysosome and Toll-like receptor signaling pathway. The other one is peroxisome, which is closely related to OS. Another category is cell cycle and cell apoptosis, including eurotrophin signaling pathway, TNF signaling pathway and ubiquitin mediated proteolysis. The fourth category is insulin signaling pathway.

Heat shock proteins (Hsp) were increased under various environment stimulus. In this study, we found the expressions of Hsp1, Hspb2 Hsp25-ps1 and Hsp1-like mRNA were significantly upregulated in the diabetic group as compared to the control group. Meanwhile, other inflammation related genes including Interferon regulatory factor 7 (Irf7), cyclin-dependent kinase 1 (Cdk1), cyclin D1. proto-oncogene protein c-fos (FOS), interleukin 1 beta (IL1B), interleukin 1 receptor, type I (IL1r1), monocyte differentiation antigen (CD14), cathepsin K (CTSK), phosphoinositide-3-kinase regulatory subunit 5 (Pik3r5), histocompatibility 2, Q region locus 1 (H2-Q1) were down-regulated. Inflammation is closely associated cell apoptosis, some of these genes were also enriched in cell cycle and apoptosis pathway, such as IL1B.

Peroxisomes are essential organelles exerting key functions in fatty acid metabolism such as the degradation of very long-chain fatty acids. Our results showed that two genes involved in lipid metabolism were significantly up-regulated. Acyl-CoA synthetase long-chain family member 4 (Acsl4) and acyl-CoA synthetase long-chain family member 6 (Acsl6) were isozyme of the long-chain fatty-acid-coenzyme A ligase family. Although differing in substrate specificity, subcellular localization, and tissue distribution, all isozymes of this family converted free long-chain fatty acids into fatty acyl-CoA esters, and thereby played a key role in lipid biosynthesis and fatty acid degradation. ACSL4 has a unique substrate specificity for arachidonic acid and modified membrane lipid composition in a manner favourable to lipid peroxidation. Hepatic ACSL4 is coregulated with the phospholipid (PL)-remodeling enzyme lysophosphatidylcholine (LPC) acyltransferase 3 to modulate the plasma triglyceride (TG) metabolism. Liver-specific knockdown of ACSL4 revealed a substantial decrease in circulating VLDL-TG levels and lipid peroxidation in mice fed a high-fat diet (Singh et al., 2019).

We also found one lipid metabolism related gene was significantly down-regulated. ATP-binding cassette, sub-family D (ALD), member 2 (ABCD2) is a member of the ALD subfamily, which is involved in peroxisomal import of fatty acids and/or fatty acyl-CoAs in the organelle. ABCD2 plays a role in the degradation of long-chain saturated and omega 9-monounsaturated fatty acids and in the synthesis of docosahexanoic acid (DHA) (Fourcade et al., 2009). The absence of ABCD2 altered expression of gene clusters associated with lipid metabolism, including PPARα signaling (Liu et al., 2014). Overexpression of ABCD2 alone prevented oxidative lesions to proteins in a mouse X-linked Adrenoleukodystrophy model (Fourcade et al., 2010) .

Our results showed that the expression of several genes involved in the glucose metabolism were significantly changed (Table S3). Insulin interacts with the insulin receptor, and the activated receptor promotes activity of the phosphoinositide-3 kinase (PI3K) enzyme. The function of differentially expressed lncRNAs were not fully understood, but their predicted target genes such as acetyl-CoA carboxylase beta (Acacb) and fructose bisphosphatase 2 (Fbp2) were involved in the glucose metabolisms.

qRT-PCR validation of the differentially expressed genes

Since OS is of significance in hepatic metabolism, a detailed inspection of genes involved in peroxosome pathway was chosen for qRT-PCR analysis to validation the RNA-seq data. qRT-PCR primers were designed based on the lncRNA sequences from mapview (https://www.ncbi.nlm.nih.gov/mapview/) (Table S4). The represented DE genes included three mRNAs, D-aspartate oxidase (Ddo), Alanine-glyoxylate aminotransferase (Agxt) and Hydroxyacid oxidase 1(Hao1), and one lncRNA Hnf4aos. qRT-PCR results were shown in Fig. 5. We found the qRT-PCR results were nearly perfect concordance with the RNA-seq results. These findings confirmed the accuracy of microarray data obtained from RNA-seq results.