Comprehensive analysis of long non-coding RNAs and mRNAs in skeletal muscle of diabetic Goto-Kakizaki rats during the early stage of type 2 diabetes

Ethical approval

The study was approved by the institutional review board of the Guangdong Key Laboratory of Laboratory Animals. All protocols were carried out in accordance with the guidelines of the Institutional Animal Care and Use Committee (IACUC) (Ethics certificate No.: IACUC2014029).

Animal breeding and tissues samples collection

Four groups of rats (diabetic male GK rats and control male Wistar rats at 3 weeks of age, diabetic male GK rats and diabetic male GK rats at 4 weeks of age, n = 10 each group), totally 40 subjects were used in this study. Rats were raised in a room with 12 h dark: 12 h light cycle, 20 to 25 °C temperature and 60 ± 5% humidity, at the SLAC Laboratory Animal Co., Ltd. (Shanghai, China) (Almon et al., 2012; Nie et al., 2017; Nie et al., 2011; Xue et al., 2011). All animals were free access to food and water. Body weight of each rat was measured weekly by weighing. Food disappearance was measured by weighing the difference in the weight of feed added and the feed remaining. The behavior of rats including feeding, drinking, sleeping and digging were observed. Blood samples were collected from the orbital plexus veins behind the eyeball using EDTA (4 mM final concentration) as an anticoagulant. Plasma was obtained from blood after centrifugation (2000 ×g, 4 °C, 15 min), divided into aliquots, and then stored at −80 °C. All rats were administered anesthesia with pentobarbital sodium (intraperitoneal, 50 mg/kg body weight), then were killed by cervical dislocation. Samples of gastrocnemius muscle of each rat were harvested, followed by rapidly frozen in liquid nitrogen, and stored at −80 °C for future studies (Nie et al., 2017; Nie et al., 2011). Six gastrocnemius muscle samples from six rats each group were selected randomly for RNA-sequencing in the present study.

Measurement of plasma glucose and insulin concentration

The automatic Dry Biochemical Analyzer FUJIFILM DRI-CHEM 7000i with GLU-PIII slides (Fujifilm, Saitama, Japan) was used to measure random plasma glucose concentration. And Thermo scientific Rat Insulin ELISA Kit (Cat#ERINS, Invitrogen, Waltham, MA, USA) was used to measured plasma insulin concentration. Assays were conducted according to the manufacturer’s instructions.

RNA extraction and sequencing

Total RNA for RNA-sequencing was extracted from red part of each gastrocnemius muscle using TRIzol Reagent (Cat#15596-018, Life Technologies, Carlsbad, CA, USA) following the manufacturer’s instructions. RNA integrity and concentration were measured by the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Ribosomal RNA was removed using Epicentre Ribo-Zero™ Gold Kits (Epicentre, Madison, WI, USA) according to the manufacturer’s instructions. RNA-sequencing was performed on Illumina HiSeq X Ten system (Illumina) following the HiSeq X Ten User Guide to generate 150 bp paired-end reads.

Analysis of differentially expressed mRNAs and lncRNAs

After quality control and filtering of low quality reads, we used STAR (Dobin & Gingeras, 2015) version 020201 to align the cleaned reads of each sample to the Rattus norvegicus reference genome (Ensembl Rnor_6.0 version 92) with the parameters of –quantMode GeneCounts –outSAMstrandField intronMotif –outSAMtype BAM SortedByCoordinate –outSAMtype BAM SortedByCoordinate –twopassMode Basic. All the corresponding annotation files of Rattus_norvegicus.Rnor_6.0.92.gtf (https://doi.org/10.6084/m9.figshare.11786277) and NONCODEv5_rat_rn6_lncRNA.gtf were downloaded from the Ensembl database54 and NONCODE version v5.0 (http://www.noncode.org/datadownload/NONCODEv5_rat_rn6_lncRNA.gtf.gz), respectively. Cufflinks were used for alignment of novel transcripts. Then the coding-probability of novel transcripts were identified by CPC2 (Kang et al., 2017), CPAT (Wang et al., 2013b) and CNCI (Sun et al., 2013). The novel transcripts with low coding-probability, or without coding-probability should meet the criteria: coding_probability score less than 0.5 in CPC2 and CPAT, and identified as noncoding by CNCI. Those novel transcripts meet criteria above with ≥ 200 bp in length and at least two exons, were defined as novel lncRNAs. Stringtie (Pertea et al., 2015) version 1.3.0 was used to assemble novel lncRNAs, annotated lncRNAs and annotated mRNAs transcripts. The novel lncRNAs were shown in Table S1. Ballgown R package (Frazee et al., 2015) version 2.10.0 was used to estimate the fragments per kilobase of exon per million fragments mapped (FPKM) of lncRNAs and mRNAs. The lncRNAs and mRNAs were filtered with FPKM < 0.5 (Moran et al., 2012). The FPKM from four groups of rats correspond to normal distribution based on the shapiro.test of Shapiro–Wilk test. The normal distribution of FPKM of four groups of rats were shown in Table S2. Next, the FPKM of GK rats at the age of 3 weeks were compared to Wistar rats at the age of 3 weeks, and GK rats at the age of 4 weeks were compared to Wistar rats at the age of 4 weeks. Thus, the differentially expressed lncRNAs (DELs) and differentially expressed mRNAs (DEGs) were obtained by Bayes-regularized t-test with an false discovery rate (FDR) correction using Cyber-T bayesreg (Kayala & Baldi, 2012). FDR < 0.05 was regarded as statistically significant. The power of test was calculated by pwr.t.test in R package pwr. The flowchart of data analysis was shown in Fig. 1.

Analysis of KEGG pathways and GO

A Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 was used to obtain the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and biological process in Gene Ontology (GO). The statistical significance threshold was P < 0.05.

Protein–protein interaction

An online database resource Search Tool for the Retrieval of Interacting Genes (STRING) (Szklarczyk et al., 2011) version 11.0 was performed to analyze protein–protein interaction of overlapping upregulated mRNAs and downregulated mRNAs between 3 and 4 weeks, respectively. After filtering disconnected nodes, we selected the minimum confidence score above 0.4 of the interaction. The confidence score was a combined score of neighborhood on chromosome, gene fusion, phylogenetic cooccurrence, homology, co-expression, experimentally determined interaction, database annotated, and automated text-mining. Those connected nodes with confidence score were downloaded for constructing the networks of protein–protein interaction.

The correlation analysis of lncRNAs and mRNAs

Python version 3.6.4 was conducted to calculate the Pearson correlation coefficients of lncRNAs and mRNAs. Those selected co-expressed lncRNA-mRNA pairs met the following criteria: correlation coefficient value > 0.9, and the absolute fold change of these DEGs and DELs ≥ 1.5. Then, pwr.r.test in R package pwr was carried out to calculate the power of the correlation.

Prediction of target genes of DELs

Since over 65% of lncRNAs were located within 10 kb of genes (Jia et al., 2010), we utilized the University of California Santa Cruz (UCSC) genome browser to identify potential cis-target genes located within 10 kb of lncRNAs (Liang et al., 2017). Then, the Basic Local Alignment Search Tool (BLAST) was applied to screen mRNAs that have complementary sequences to lncRNAs, followed by RNAplex (Liang et al., 2017) to identify trans-regulated target genes of lncRNAs. Subsequently, DELs and their corresponding target DEGs were obtained.

Prediction of microRNAs both targeted to key DELs and DEGs

We predicted the target microRNAs (miRNAs) of DELs and DEGs in key lncRNA-mRNA pairs. TargetScan (Agarwal et al., 2015) was carried out to predict the target miRNAs of key DEGs, and the RNAhybrid (Rehmsmeier et al., 2004) was applied to predict the target miRNAs of key DELs.

Construction of interaction network

Cytoscape version 3.6.1 was exerted to construct the networks of lncRNA-mRNA pairs and protein–protein interaction networks with those downloaded files including connected nodes with confidence score.

Statistical analysis

All data were expressed as mean ± standard deviation (SD) unless otherwise noted. The significant difference was measured using a two-tailed student t-test. P < 0.05 was considered statistically significant.

The characteristics of rats

The plasma glucose concentration of GK rats was significantly higher than that of control Wistar rats at 3 and 4 weeks of age (P < 0.001, Table 1), which was in accordance with previous research (Ando et al., 2018). Besides, the plasma glucose concentration of 3-week-old GK rats was significantly higher than that of 4-week-old GK rats (P < 0.001, Table 1).

Differentially expressed lncRNAs and mRNAs

In total, we got 438 and 1,000 differentially expressed mRNAs (DEGs) between GK and Wistar rats at 3 and 4 weeks of age, respectively (false discovery rate, FDR < 0.05) (Fig. 2A, Table S3). There were 401 and 746 differentially expressed lncRNAs (DELs) in GK rats compared with Wistar rats at 3 and 4 weeks of age, respectively (FDR < 0.05) (Fig. 2B, Table S4). Among the DEGs, 141 overlapping upregulated mRNAs and 103 overlapping downregulated mRNAs were detected between 3 and 4 weeks (Fig. 2A). A total of 91 overlapping upregulated lncRNAs and 124 overlapping downregulated lncRNAs were found between 3 and 4 weeks (Fig. 2B). From the results of enrichment pathway analysis of DEGs at 3 and 4 weeks, we found the insulin resistance pathway was the only one overlapping pathway between 3 and 4 weeks among top 10 KEGG pathways (Figs. 2C and 2D). The DEGs enriched in insulin resistance pathway were sterol regulatory element binding transcription factor 1 (Srebf1, also known as Srebp1c), solute carrier family 27 member 1 (Slc27a1), protein kinase C, theta (Prkcq), cAMP responsive element binding protein 3-like 1 (Creb3l1), forkhead box O1 (Foxo1), TBC1 domain family, member 4 (Tbc1d4, also termed as AS160), and carnitine palmitoyltransferase 1A (Cpt1a).

To obtain the interaction of the proteins encoded by 141 overlapping upregulated genes and 103 overlapping downregulated genes between 3 and 4 weeks, we analyzed protein–protein interaction of these proteins using STRING. Next, we constructed the protein–protein interaction network (Figs. 3A and 3B). Then, the top 10 mRNAs according to the node degree among network and their corresponding node degrees were listed in Table S5 . The network among top 10 upregulated and downregulated node mRNAs were shown in Figs. 3C and 3D. The dysregulated genes Srebf1, Slc27a1, Foxo1 and Cpt1a that enriched in insulin resistance pathway also existed in the network of top 10 upregulated and downregulated node mRNAs (Figs. 3C and 3D), which indicating that these four genes might be important for the development of hyperglycemia and T2D in GK rats at the age of 3 and 4 weeks.

The co-expressed lncRNA-mRNA pairs with high correlation coefficients

To investigate the potential function of these DELs, we performed lncRNA-mRNA co-expression network analysis. After filtering, a total of 901 co-expressed lncRNA-mRNA pairs with high correlation coefficients were selected, including 136 DEGs and 120 DELs (Table S6). 136 DEGs were enriched in two KEGG pathways (P < 0.05), including transcriptional misregulation in cancer and pathways in cancer. But both pathways were not related to T2D. 136 DEGs were enriched in biological processes (Fig. 4A). Among these DEGs, 2 DEGs (pyruvate dehydrogenase kinase 4, Pdk4 and Cpt1a) were enriched in “regulation of fatty acid oxidation”, and 3 DEGs (Foxo1, Pdk4 and SH2B adaptor protein 2, Sh2b2) were enriched in “insulin receptor signaling pathway”. The top 10 nodes ranked by degrees in co-expressed lncRNA-mRNA network were consist of 7 mRNAs and 3 lncRNAs (Table 2). The network of these dysregulated mRNAs (Cep19, Cpt1a, Ephx2, Foxo1, Pdk4, Sh2b2 and Stc2) and their co-expressed lncRNAs were shown in Fig. 4B. And the expression of key co-expressed lncRNA-mRNA pairs were shown in Figs. 4C and 4D. Notably, the dysregulated genes Cpt1a, Foxo1 and Pdk4 also appeared in the network of top 10 upregulated and downregulated mRNAs (Figs. 4C and 4D), indicating that these mRNAs might associate with hyperglycemia and T2D in GK rats at the age of 3 and 4 weeks.

However, other genes involved in fatty acid transport and β-oxidation were not significantly dysregulated in GK rats at the age of 3 and 4 weeks (Table S7). In our study, genes related to glycolysis and glycogen synthesis were not dysregulated in GK rats at 3 and 4 weeks of age (Table S7).

The predicated target mRNAs of differentially expressed lncRNAs

To identify the potential role of dysregulated lncRNAs in the development of hyperglycemia and T2D in GK rats at the age of 3 and 4 weeks, we predicted their cis- and trans-target mRNAs. A total of 15 predicted cis-target DEGs and 88 predicted trans-target DEGs were obtained in DELs at 3 weeks (Table S8). There were 31 predicted cis-target DEGs and 382 predicted trans-target DEGs in DELs at 4 weeks (Table S9). There were 32 overlapping DEL-target DEGs between 3 and 4 weeks (Fig. 5A, Table S10). Network analysis for these overlapping DEGs found out 32 lncRNA-mRNA pairs, including 18 DEGs and 19 DELs (Fig. 5B). Among the 18 DEGs, interleukin 15 (Il15), F-box and WD repeat domain containing 7 (Fbxw7) and uncoupling protein 3 (Ucp3) were related to increased glycaemia (Gray & Kamolrat, 2011; Zhao et al., 2018), glucose intolerance (Fujimoto et al., 2019), and increased fatty acid oxidation (Bezaire et al., 2005). Among those 32 lncRNA-target mRNA pairs, 5 lncRNA-target mRNA pairs (NONRATG014028.2-Pim1, NONRATG011882.2-Il15, NONRATG013497.2-Fbxw7, NONRATG011747.2-Mrps35, and MSTRG.1662-Ucp3) appeared in 901 co-expressed lncRNA-mRNA pairs with high correlation coefficients, suggesting these dysregulated lncRNA-mRNA pairs might involve in the hyperglycemia and T2D of GK rats at the age of 3 and 4 weeks. And the relative expression of NONRATG011882.2-Il15, NONRATG013497.2-Fbxw7, and MSTRG.1662-Ucp3 pairs were shown in Fig. 5C.

The target microRNAs (miRNAs) of key DEGs and DELs

To explore the role of lncRNAs in the expression of mRNAs, we predicted the target miRNA of key DELs and DEGs. Then we obtained the overlapping miRNAs targeted both DEGs and DELs in key lncRNA-mRNA pairs (Fig. 6), which provided miRNAs linkers between DELs and DEGs. We found that rno-miR-139-5p, rno-miR-486 and rno-miR-93-5p target both MSTRG.14356 and Foxo1 (Fig. 6). We got three target miRNAs rno-miR-20b-5p, rno-miR-27a-3p and rno-miR-17-5p of MSTRG.2584 and Foxo1 (Fig. 6). The target miRNAs of Stc2 and MSTRG.2584 were rno-miR-24-3p, rno-miR-532-5p, rno-miR-181a-5p and rno-miR-181b-5p (Fig. 6). MiRNAs rno-miR-195-5p, rno-miR-181b-5p, rno-miR-23b-3p, rno-miR-139-5p and rno-miR-23a-3p were the target miRNAs of MSTRG (Fig. 6).12678 and Pdk4. Three miRNAs including rno-miR-34a-5p, rno-miR-125a-5p and rno-miR-125b-5p targeted MSTRG.1662 and Ucp3 (Fig. 6). Additionally, the target miRNA of NONRATG013497.2 and Fbxw7 was rno-miR-24-3p (Fig. 6). And rno-miR-326-3p was the target of NONRATG011882.2 and Il15 (Fig. 6).