Identification of Potential Therapeutic Targets and Pathways of Liraglutide Against Type 2 Diabetes Mellitus (T2DM) Based on Long Non-Coding RNA (lncRNA) Sequencing.

BACKGROUND The aim of this study was to explore the potential therapeutic targets and pathways of liraglutide against type 2 diabetes mellitus (T2DM) in streptozotocin-induced diabetic rats based on lncRNA sequencing. MATERIAL AND METHODS Male Wistar rats were randomly divided into 3 groups: the control group (n=10), the T2DM model group (high-sugar and high-fat diet, and streptozotocin-induced, n=11), and the liraglutide group (model plus liraglutide, n=10). After 8 weeks of drug treatment, lncRNA sequencing was used to identify the lncRNA therapeutic targets and their related protein-coding genes of liraglutide against T2DM, which were further studied by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis to determine the major biological processes and pathways involved in the action of liraglutide treatment. Lastly, several lncRNA targets were randomly detected based on quantitative real-time polymerase chain reaction (QRT-PCR) to verify the accuracy of sequencing results. RESULTS A total of 104 lncRNA targets of liraglutide against T2DM were screened, with 27 upregulated and 77 downregulated, including NONRATT030354.2, MSTRG.1456.6, and NONRATT011758.2. The major biological processes involved were glucose and lipid metabolism and amino acid metabolism. Liraglutide had a therapeutic effect in T2DM, mainly through the Wnt, PPAR, amino acid metabolism signaling, mTOR, and lipid metabolism-related pathways. CONCLUSIONS In this study, we screened 104 lncRNA therapeutic targets and several signaling pathways (Wnt, PPAR, amino acid metabolism signaling pathway, mTOR, and lipid metabolism-related pathways) of liraglutide against T2DM based on lncRNA sequencing.

Background

Type 2 diabetes mellitus (T2DM) is been a chronic metabolic disease with high incidence and high disability rates and involves a series of complications. T2DM has become a severe public health issue throughout the world, especially in developing countries [1]. Nevertheless, the pathogenesis of T2DM is extremely complicated and unclear, and there are few currently available clinical therapeutic drugs. There have been various novel drugs for the treatment of diabetes, including GLP-1 analogues such as liraglutide and exenatide, the ultra-long-acting basic insulin analogues such as Degu insulin, and the DPP-4 inhibitors such as Zafatek (trelagliptin succinate) [2-4]. As the first human glucagon-like peptide-1 (GLP-1) analogue worldwide, liraglutide (molecular formula: C172H265N43O51) shares 97% homology with the natural human GLP-1 [5], which is an endogenous incretin hormone that promotes the glucose concentration-dependent secretion of insulin in pancreatic β cells. Liraglutide replaces the 34th lysine of native GLP-1 with chlorinated acid and inserts a glutamate-mediated 16-carbon palmitoyl fatty acid side chain to lysine at position 26 [5]. Consequently, liraglutide not only maintains the efficacy of natural GLP-1, but also possesses strong chemical stability due to the presence of fatty acid side chains that resist degradation by DPP-IV and with the half-life of 12–14 h [6]. Owing to its unique chemical structure, liraglutide has an excellent hypoglycemic effect on T2DM via once-daily injection, which also reduces the endothelial endoplasmic reticulum stress and insulin resistance and helps with weight loss and cardiovascular protection [5-7]. To assess the multifaceted regulatory mechanisms underlying the therapeutic effect on T2DM, we focussed on long noncoding RNA (lncRNA) sequencing in this study. lncRNAs, which are noncoding RNAs more than 200 nucleotides in length, have become an important topic in genetics research [8]. Many studies have revealed that lncRNAs are closely involved in T2DM, endothelial endoplasmic reticulum stress, insulin secretion, and islet cells [9,10]. Exploring the lncRNA therapeutic targets, biological processes, and pathways of liraglutide can widen the horizons for the diagnosis and treatment of T2DM, which has remarkable clinical significance.

We first investigated various lncRNA therapeutic targets of liraglutide based on lncRNA transcriptomics in streptozotocin-induced T2DM model rats. Rats were fed a high-fat and high-sugar diet for 8 weeks to induce insulin resistance, and then injected with STZ to destroy the function of pancreatic cells. The T2DM rat model produced using this method is very similar to the human disease model. We assessed differential expression (DE) of lncRNAs in 3 groups for further research. Additionally, GO Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to identify the major pathways of biological processes and functions involved in treatment with liraglutide. We used QRT-PCR to evaluate the expression of lncRNA therapeutic targets to verify the accuracy of sequencing. The study workflow is shown in Figure 1.

Figure 1

Study workflow.

Material and Methods

The ethics of animal experiments and model establishment

All experiments in the present study were approved by the Animal Ethics Committee of the Affiliated Hospital of Shandong University of Traditional Chinese Medicine (approval no. AWE-2-19-001). Thirty-four male Wistar rats, 4 weeks old, weighing 130±10 g, were purchased from Beijing Vital River Company (SCXK 2016-0006, no. 11400700377281). The ambient temperature was maintained at 18–22°C and humidity was 40–60%, with a 12-h light/dark cycle. Ten rats were randomly selected as the control group (C group). The remaining rats were fed an adaptive diet for 1 week, followed by 8 weeks of high-sugar and high-fat diet. Then, these rats were intraperitoneally injected with 35 mg/kg STZ (sigma, # 18883-66-4, USA). At 72 h after the intervention, when the concentration of fasting blood glucose was more than 16.7 mmol/L as measured by a Roche glucometer, it was considered that the model was successfully established. Excluding rats that were not successfully modeled, the remaining model rats were further randomly separated into 2 groups on the same day: the model group (M group, n=11) and the liraglutide group (L group, n=10). On the next day, we began treatment. Liraglutide (0.11 mg/kg·d, Victoza, Denmark) was subcutaneously administered to rats in the liraglutide group, and the rats of the control group and the model group were treated with the same amount of physiological saline. After 8 weeks of treatment, rats were dissected and the pancreases were obtained for further analysis.

HE staining and immunofluorescence

The pancreatic tissue of each group was fixed with 10% poly-methyl, routinely embedded in paraffin, sliced 4-μm thick, dewaxed, stained, sealed with neutral gum, observed under light microscope, and photographed.

The paraffin sections were deparaffinized, antigen-repaired, and rinsed, followed by treating with an auto-fluorescence quencher. After washing, BSA was incubated with the sections for 30 min and the blocking solution was discarded. Next, the insulin and glucagon primary antibody (Servicebio, #GB13121; #GB13097, China) were added dropwise to the system, followed by overnight incubation at 4°C. Afterwards, the secondary antibody (Servicebio, #GB21301; #GB25303, China) was added to the system, cultured at room temperature for 50 min, then treated with the DAPI dye solution and incubated at room temperature for 10 min. The anti-fluorescence quenched the sealer, and the images were visualized under a fluorescence microscope.

lncRNA sequencing

RNA extraction and quality control

The samples of rats were randomly collected from each group in triplicate. RNA was extracted using the miRNeasy Mini Kit (Qiagen, #74106, Germany) in strict accordance with the manufacturer’s instructions. The obtained total RNA was quality-checked with an Agilent Bioanalyzer 2100 (Agilent technologies, USA) and quantified with a Qubit®3.0 Fluorometer and NanoDrop One spectrophotometer.

Library construction and transcriptome sequencing

After purification, divalent cations were used for fragmentation of mRNA at 94°C for 8 min. Reverse transcriptase and random primers were used for the replication of the obtained RNA fragments to first-strand cDNA, followed by synthesis of second-strand cDNA using DNA Polymerase I and RNase H. Subsequently, the cDNA fragments underwent end repair, in which a single ‘A’ base was added, and then ligated to the adapters. The products were purified and enriched bases on PCR for the establishment of the cDNA library. The cDNA in the developed libraries were further quantified based on Qubit® 2.0 Fluorometer (Life Technologies, USA) and verified with an Agilent 2100 bioanalyzer (Agilent Technologies, USA) for the confirmation of insert size and evaluation of the molar concentration. cBot was used for the separation of clusters when the library was diluted to 10 pM, and were subsequently sequenced with an Illumina NovaSeq 6000 (Illumina, USA). The library development and sequencing were conducted at Shanghai Sinomics Corporation.

Data acquisition, lncRNA identification, and expression analysis

For data analysis, the raw data were first filtered, removing the low-quality reads. To acquire the novel transcripts, all the assembled transcript isoforms were compared with the known protein-encoding transcripts in rats by cuffcompare. Putative lncRNAs were defined as the novel transcripts satisfying the following criteria: transcript length ≥200 nt; the ability of coding potential and encoding proteins; ORF <300 bp; no records in the Pfam database; combining the results of CNCI, Pfam, and CPC, the score of CPC and CNCI <0. Statistically significant DE lncRNAs were screened based on p value and fold change (p value <0.05, FC >2 or FC <0.5). Afterwards, the hierarchical clustering and correlation analysis were performed with scripting. Volcano plots and Venn diagrams were constructed to determine the differentially expressed genes, and the chromosomal localization of these lncRNAs was observed.

Exploration of lncRNA therapeutic targets and functional analysis

Based on the P value and FC, we found differentially expressed (DE) lncRNAs between the 3 sets of samples, and the results were represented by Venn diagrams and volcano plots, and chromosome mapping analysis was performed on these DE lncRNAs. The intersection between the upregulated transcripts in the T2DM model vs. control group and the downregulated transcripts in the liraglutide vs. model groups was investigated. Additionally, the downregulated lncRNAs in T2DM and the upregulated ones after liraglutide treatment were used to generate another set of intersections. Based on the intersections, the mechanisms by which liraglutide reversed the pathophysiological changes in T2DM were identified.

Since the lncRNAs were mainly functionalized by encoding genes with proteins [11], the known related protein-coding genes needed to be explored. A gene transcribed within a 10-kbp window upstream or downstream of the lncRNAs was considered to be a cis-acting target gene [12]. RNAplex software was used for the prediction of trans-acting target genes.

GO and KEGG enrichment analyses were performed to identify the biological process and signaling pathways of liraglutide in the treatment of T2DM. The lncRNAs at the intersections and the relevant biological functions and pathways help in interpreting the potential mechanisms of liraglutide in treating T2DM.

QRT-PCR validation

To validate the sequence data, 6 DE LncRNAs were chosen for the QRT-PCR analysis, including:

NONRATT015614.2

(forward primer: 5GGACCCTGGCCTTCCTCTA3′; reverse primer: 5′GTGGCTGAACTTTGATTTCGTAT3′);

NONRATT004911.2

(forward primer: 5′TGAAGACGCAGAGTAAATCCT3′; reverse primer: 5′TCTACCACTGACCTAAATCCC3′);

NONRATT018630.2

(forward primer: 5′GCTTTCTGGGTATGTCTTCTCC3′; reverse primer: 5′CTGGTCTTCCGTAAGTCTTGTC3′);

NONRATT029906.2

(forward primer: 5′CTGTTGGGACTGTTGGAAA3′; reverse primer: 5′CCCTAAGCGAAATAAAGCA3′);

NONRATT024782.2

(forward primer: 5′ATCTGATGCCCTCTTCTGGTGT3′; reverse primer: 5′ATGTATCCTGAGCTGGCCTTTA3′);

NONRATT026027.2

(GCATCCTACCCACCCTCACT, GCCTCTGATGGCTGGTCTTT).

The primer sequences were designed by Sangon Biotech Co. (Shanghai, China). Ct values were normalized to GAPDH, and ΔΔCt was calculated as (ΔCt sample –ΔCt reference), and the 2−ΔΔCt method was used to show the relative expression.

Statistical analysis

SPSS 18.0 software was used for data analysis, and all results are presented as mean±SD. The two-tailed Student’s t test was used for the data comparison of 2 groups. P<0.05 was considered to be statistically significant.

Results

HE staining and immunofluorescence assay

The results of HE staining in comparison with the control group showed the islet morphology of the pancreatic tissue of the DM group was irregular. Specifically, the islets were obviously atrophied, the contour was less rounded, the number of islet cells were dramatically decreased, and the boundary with the exocrine glands was ambiguous and disordered. Compared to the model group, after treatment with liraglutide, the islet morphology in the pancreas was improved, with a distinct outline, and the morphology was similar to that of normal tissues. In addition, the islets in the island were markedly upregulated, with clearer and more regular boundaries of exocrine glands. The results revealed that the islet morphology was conspicuously improved and the number of islet cells was dramatically increased by the administration of liraglutide (Figure 2A).

Figure 2

(A) HE staining results of 3 groups (magnification ×400). (B) The immunofluorescence analyses of different groups (magnification ×400, the red fluorescence shows insulin and the green show glucagon).

In comparison with the control group, the immunofluorescence results suggested that there were far fewer β cells in the middle of islets of the model group, and the cells were unevenly distributed. There was no significant change in pancreatic α cells. In contrast to the model group, β cells in the pancreas of rats in the liraglutide group were prominently upregulated and gathered in the center of islets, while pancreatic α cells had no obvious changes. Our results show that liraglutide can dramatically increase the number of β cells and improve insulin secretion (Figure 2B).

Identification of DE lncRNAs among 3 groups and chromosomal localization analysis

According to P value and FC, compared to the control group, 104 DE lncRNAs were upregulated and 177 were downregulated in the T2DM model group. According to chromosomal localization, although lncRNAs were abundant and present on every chromosome after T2DM-onset, DE lncRNAs were mainly localized on chromosomes 1, 2, 7, and 12, as revealed in Figure 3. After liraglutide treatment, there were 112 upregulated and 203 downregulated DE lncRNAs. Figure 3 shows that the localizations of DE lncRNAs were generally on chromosomes 1, 2, and 7. In Figure 4A and 4B, the overall distributions of DE lncRNAs among the 3 groups could be clearly observed based on the Venn diagrams and volcano plots.

Figure 3

The chromosomal location of the 3 sets of samples (the abscissa is the position of the chromosome and the vertical coordinate is the number of chromosomes).

Figure 4

(A) Venn diagram of the intersection of the 3 groups. (B) The volcano plots revealed the DE lncRNAs between the model group vs. control group and the liraglutide group vs. model group. The red dots show upregulated lncRNA, blue dots show downregulated lncRNA, and gray shows no significant differences.

The screening of potential lncRNA therapeutic targets

The disease-changing lncRNAs reversed by liraglutide were screened to discover which were the precise therapeutic targets for the treatment of T2DM with liraglutide. Consequently, 104 lncRNAs therapeutic targets were obtained (Table 1). There were 27 upregulations after the administration of liraglutide, including MSTRG.1456.6, MSTRG.18306.1, MSTRG.1530.94, MSTRG.18301.2, MSTRG.11495.5, MSTRG.6009.1, NONRATT030354.2, MSTRG.5180.1, MSTRG.1530.73, NONRATT024782.2, NONRATT029906.2, NONRATT003195.2, NONRATT022556.2, MSTRG.13735.3, NONRATT027557.2, MSTRG.1530.123, NONRATT030577.2, MSTRG.1530.128, NONRATT018153.2, NONRATT025333.2, NONRATT003365.2, ENSRNOT00000085168, NONRATT007884.2, NONRATT003631.2, MSTRG.7668.9, NONRATT016299.2, and NONRATT003632.2. There were 77 downregulations after liraglutide treatment, including MSTRG.5180.22, NONRATT011758.2, NONRATT021220.2, MSTRG.20445.6, NONRATT026723.2, MSTRG.11301.18, MSTRG.9982.1, NONRATT023797.2, MSTRG.11959.1, MSTRG.9983.1, NONRATT015614.2, MSTRG.11301.14, MSTRG.5181.1, NONRATT007668.2, MSTRG.19274.3, MSTRG.19274.4, MSTRG.9984.1, MSTRG.10769.1, MSTRG.5180.59, NONRATT021420.2, MSTRG.8589.1, NONRATT026017.2, MSTRG.1530.117, NONRATT007560.2, NONRATT015190.2, NONRATT030883.2, MSTRG.1456.9, NONRATT023893.2, ENSRNOT00000083718, NONRATT012600.2, MSTRG.3668.1, MSTRG.11246.1, MSTRG.19303.4, NONRATT020829.2, NONRATT018630.2, NONRATT024781.2, MSTRG.19303.2, NONRATT028050.2, NONRATT024530.2, MSTRG.5245.16, MSTRG.13735.4, MSTRG.1351.2, MSTRG.1530.113, NONRATT010566.2, MSTRG.2100.4, NONRATT008075.2, NONRATT013576.2, MSTRG.19303.12, MSTRG.19303.16, NONRATT015129.2, MSTRG.11301.20, MSTRG.1530.9, MSTRG.19303.17, NONRATT004911.2, NONRATT016292.2, MSTRG.5245.14, NONRATT021144.2, MSTRG.6299.1, MSTRG.20673.3, NONRATT008300.2, MSTRG.6007.6, MSTRG.19303.18, ENSRNOT00000087028, NONRATT023036.2, NONRATT008225.2, MSTRG.21671.1, NONRATT024339.2, NONRATT001395.2, MSTRG.2393.1, MSTRG.19188.9, NONRATT020841.2, MSTRG.5180.17, NONRATT002773.2, NONRATT015765.2, NONRATT004048.2, NONRATT014284.2, and NONRATT026027.2. We also further clustered these 104 lncRNAs and produced clustering maps (Figure 5).

Table 1

Specific information on 104 lncRNAs, including FC, P value, regulation, and subtypes, among 3 groups.

lncRNA	log2FC		P value		Regulation		Subtype
	M vs. C	L vs. M	M vs. C	L vs. M	M vs. C	L vs. M
MSTRG.1456.6	−4.743000313	5.513555708	1.48E-11	5.06E-14	Down	Up	Antisense LncRNA
MSTRG.18306.1	−5.813267702	5.569096364	4.87E-11	0.00011411	Down	Up	Antisense LncRNA
MSTRG.1530.94	#NAME?	inf	2.40E-09	5.01E-08	Down	Up	Antisense LncRNA
MSTRG.18301.2	−3.890090199	3.707173112	4.43E-08	0.003776484	Down	Up	Intergenic LncRNA
MSTRG.11495.5	#NAME?	inf	1.25E-07	0.00030484	Down	Up	Intergenic LncRNA
MSTRG.6009.1	#NAME?	inf	8.30E-07	0.000318274	Down	Up	Intergenic LncRNA
NONRATT030354.2	−5.268148594	4.59178329	2.10E-06	0.000425377	Down	Up	Exonic_sense
MSTRG.5180.1	−7.759506305	7.404002656	2.72E-06	7.66E-05	Down	Up	Intergenic LncRNA
MSTRG.1530.73	−5.36827079	5.666529456	3.72E-06	0.000443118	Down	Up	Antisense LncRNA
NONRATT024782.2	−4.627137132	4.702185072	2.90E-05	3.83E-05	Down	Up	Exonic_sense
NONRATT029906.2	−3.423342356	3.664979392	0.000816779	0.001816205	Down	Up	Exonic_sense
NONRATT003195.2	−1.828257785	1.2137236	0.001190049	0.038764384	Down	Up	Exonic_sense
NONRATT022556.2	−2.17175632	1.875589531	0.001977024	0.009853911	Down	Up	Intronic_sense
MSTRG.13735.3	#NAME?	inf	0.003237789	0.00022438	Down	Up	Intergenic LncRNA
NONRATT027557.2	−2.294101628	1.449048513	0.003321136	0.049900577	Down	Up	Intergenic
MSTRG.1530.123	−4.679010037	4.742641573	0.003890266	0.001129404	Down	Up	Antisense LncRNA
NONRATT030577.2	−2.422284192	2.948953489	0.004104782	2.68E-06	Down	Up	Exonic_sense
MSTRG.1530.128	−3.574914702	3.898916074	0.006809672	0.009434631	Down	Up	Antisense LncRNA
NONRATT018153.2	−1.519952453	1.789990773	0.010235007	0.002359126	Down	Up	Intronic_sense
NONRATT025333.2	−2.789473738	3.524107076	0.012374129	0.000351516	Down	Up	Exonic_sense
NONRATT003365.2	−2.171464786	3.172344935	0.012569546	1.02E-05	Down	Up	Exonic_sense
ENSRNOT00000085168	−1.884214285	1.536023487	0.015376251	0.04164753	Down	Up	Bidirectional
NONRATT007884.2	−1.49593349	2.359567666	0.018515518	0.000288795	Down	Up	Exonic_sense
NONRATT003631.2	−1.320304084	2.024315038	0.023267283	0.000607493	Down	Up	Intergenic
MSTRG.7668.9	−1.183499239	1.286391685	0.023999121	0.02178103	Down	Up	Antisense LncRNA
NONRATT016299.2	−1.500122512	1.334528093	0.029757841	0.034135171	Down	Up	Exonic_sense
NONRATT003632.2	−1.105928494	1.325361911	0.030854214	0.012459209	Down	Up	Intergenic
MSTRG.5180.22	7.446809946	−2.382791481	1.91E-17	0.037115032	Up	Down	Intergenic LncRNA
NONRATT011758.2	7.860975308	−2.143691132	1.27E-15	0.007919459	Up	Down	Exonic_sense
NONRATT021220.2	6.43043182	−6.206406186	2.12E-15	3.59E-18	Up	Down	Exonic_sense
MSTRG.20445.6	inf	−3.417817831	4.39E-12	6.67E-05	Up	Down	Intergenic LncRNA
NONRATT026723.2	4.535584232	−4.763048926	7.06E-11	5.62E-11	Up	Down	Exonic_sense
MSTRG.11301.18	9.658751684	#NAME?	8.36E-11	1.57E-11	Up	Down	Intron LncRNA
MSTRG.9982.1	inf	−3.225521198	5.03E-10	0.000756633	Up	Down	Intergenic LncRNA
NONRATT023797.2	4.507660401	−2.549124425	3.92E-09	9.23E-05	Up	Down	Bidirectional
MSTRG.11959.1	4.174470175	−3.778946282	9.00E-09	2.48E-09	Up	Down	Intergenic LncRNA
MSTRG.9983.1	inf	−4.241192717	1.15E-08	0.000275468	Up	Down	Intergenic LncRNA
NONRATT015614.2	4.549936726	−3.350159665	2.65E-08	4.06E-06	Up	Down	Bidirectional
MSTRG.11301.14	inf	#NAME?	2.95E-07	2.97E-07	Up	Down	Intron LncRNA
MSTRG.5181.1	inf	−2.148052314	2.78E-07	0.013877698	Up	Down	Intergenic LncRNA
NONRATT007668.2	4.274962562	−2.90779253	2.88E-07	0.000124786	Up	Down	Intronic_sense
MSTRG.19274.3	inf	−4.377122746	6.74E-07	0.001288275	Up	Down	Intron LncRNA
MSTRG.19274.4	2.713718919	−1.729000444	3.83E-06	0.00372144	Up	Down	Antisense LncRNA
MSTRG.9984.1	7.333129947	−4.913864356	3.95E-06	0.000129424	Up	Down	Intergenic LncRNA
MSTRG.10769.1	3.017684836	−2.136276284	6.11E-06	0.001250674	Up	Down	Intron LncRNA
MSTRG.5180.59	2.691302066	−2.116419082	1.56E-05	0.000891163	Up	Down	Intergenic LncRNA
NONRATT021420.2	inf	−2.482990883	6.74E-07	0.001825953	Up	Down	Intergenic
MSTRG.8589.1	2.428030991	−1.270356705	2.19E-05	0.026715661	Up	Down	Antisense LncRNA
NONRATT026017.2	3.067981167	−1.587203641	2.33E-05	0.024197572	Up	Down	Exonic_sense
MSTRG.1530.117	inf	−4.713575915	3.82E-05	0.000637472	Up	Down	Intron LncRNA
NONRATT007560.2	4.329709516	−5.169573424	4.04E-05	4.91E-06	Up	Down	Exonic_sense
NONRATT015190.2	3.885078028	−3.376452176	0.000116874	0.000788245	Up	Down	Exonic_sense
NONRATT030883.2	2.912494377	−2.974995522	0.000154878	0.00014227	Up	Down	Exonic_sense
MSTRG.1456.9	inf	#NAME?	0.000215436	0.000207305	Up	Down	Antisense LncRNA
NONRATT023893.2	2.140679407	−1.330633497	0.000222945	0.034025848	Up	Down	Exonic_sense
ENSRNOT00000083718	inf	#NAME?	0.00030439	0.000290699	Up	Down	Intergenic
NONRATT012600.2	3.93966141	−3.728798986	0.00034416	0.000440445	Up	Down	Exonic_sense
MSTRG.3668.1	2.603023236	−2.045267464	0.000591145	0.004829552	Up	Down	Antisense LncRNA
MSTRG.11246.1	2.066124224	−3.050106783	0.000605332	0.000225883	Up	Down	Intergenic LncRNA
MSTRG.19303.4	5.60089652	−4.536826565	0.000693249	0.001963455	Up	Down	Intron LncRNA
NONRATT020829.2	4.78318059	−4.614368758	0.000718728	0.000969879	Up	Down	Exonic_sense
NONRATT018630.2	2.22100913	−1.971267322	0.000759833	0.004717386	Up	Down	Intronic_sense
NONRATT024781.2	2.822227607	−1.783103169	0.000773861	0.026243132	Up	Down	Exonic_sense
MSTRG.19303.2	inf	#NAME?	0.000885874	0.000848526	Up	Down	Intron LncRNA
NONRATT028050.2	inf	#NAME?	0.000909817	0.000891584	Up	Down	Intronic_sense
NONRATT024530.2	2.117262597	−2.68507487	0.000950953	0.000178981	Up	Down	Intergenic
MSTRG.5245.16	6.731430255	−6.988093941	0.001157906	0.001004914	Up	Down	Intergenic LncRNA
MSTRG.13735.4	inf	#NAME?	0.001196225	0.0011795	Up	Down	Intergenic LncRNA
MSTRG.1351.2	inf	−6.385972563	0.001235499	0.003415514	Up	Down	Antisense LncRNA
MSTRG.1530.113	4.704157435	−3.313542754	0.001264395	0.005825531	Up	Down	Antisense LncRNA
NONRATT010566.2	3.365138578	−3.490224709	0.001393601	0.001152729	Up	Down	Exonic_sense
MSTRG.2100.4	inf	#NAME?	0.001425832	0.001252196	Up	Down	Antisense LncRNA
NONRATT008075.2	3.032888346	−4.708810304	0.001477835	1.46E-07	Up	Down	Exonic_sense
NONRATT013576.2	4.886580918	−5.075248519	0.001472997	0.001269298	Up	Down	Exonic_sense
MSTRG.19303.12	1.985103384	#NAME?	0.001507159	0.001579596	Up	Down	Intron LncRNA
MSTRG.19303.16	inf	#NAME?	0.001758442	0.001702558	Up	Down	Intron LncRNA
NONRATT015129.2	1.827183471	−2.205818879	0.001834647	0.00061533	Up	Down	Exonic_sense
MSTRG.11301.20	3.525195633	−4.661747539	0.001873192	0.01887575	Up	Down	Intron LncRNA
MSTRG.1530.9	6.696409457	−4.661747539	0.001890656	0.01887575	Up	Down	Antisense LncRNA
MSTRG.19303.17	inf	#NAME?	0.002613953	0.002460638	Up	Down	Intron LncRNA
NONRATT004911.2	1.960857645	−2.298512462	0.002805181	0.000234504	Up	Down	Intronic_sense
NONRATT016292.2	4.135749829	−3.324687904	0.003994089	0.012925679	Up	Down	Exonic_sense
MSTRG.5245.14	inf	#NAME?	0.004746998	0.00476209	Up	Down	Intergenic LncRNA
NONRATT021144.2	3.143234076	−3.290826659	0.006037138	0.005515802	Up	Down	Exonic_sense
MSTRG.6299.1	inf	#NAME?	0.006940124	0.006942518	Up	Down	Intron LncRNA
MSTRG.20673.3	1.54177048	−2.099779452	0.009879226	0.001135172	Up	Down	Intron LncRNA
NONRATT008300.2	1.59131684	−1.423986135	0.011099917	0.010774052	Up	Down	Exonic_sense
MSTRG.6007.6	3.800748566	−3.929021582	0.01191241	0.011230168	Up	Down	Intergenic LncRNA
MSTRG.19303.18	inf	#NAME?	0.012213679	0.012736492	Up	Down	Intron LncRNA
ENSRNOT00000087028	1.387888757	−1.194919927	0.013785164	0.039176636	Up	Down	Intergenic
NONRATT023036.2	1.291422763	−1.436154296	0.018713238	0.030642508	Up	Down	Exonic_sense
NONRATT008225.2	2.277530254	−3.837996506	0.020896786	9.69E-05	Up	Down	Intergenic
MSTRG.21671.1	1.451298167	−3.290826659	0.021739051	0.005515802	Up	Down	Antisense LncRNA
NONRATT024339.2	1.555328009	−2.570512467	0.024104945	0.001393028	Up	Down	Exonic_sense
NONRATT001395.2	1.58713361	−1.347997983	0.026086395	0.046456227	Up	Down	Intronic_sense
MSTRG.2393.1	1.467819139	−1.768429489	0.027831951	0.018507121	Up	Down	Antisense LncRNA
MSTRG.19188.9	2.458563141	−1.984620065	0.028304987	0.036176215	Up	Down	Exonic_sense
NONRATT020841.2	1.633176918	−2.012524762	0.030734688	0.010409861	Up	Down	Exonic_sense
MSTRG.5180.17	3.851982823	−3.230608763	0.031543156	0.038975462	Up	Down	Exonic_sense
NONRATT002773.2	1.404131237	−1.488955374	0.031592231	0.019918394	Up	Down	Intronic_sense
NONRATT015765.2	1.854217567	−2.438880302	0.036617281	0.001020115	Up	Down	Exonic_sense
NONRATT004048.2	1.320307692	−2.02571342	0.040127517	0.002966605	Up	Down	Intronic_sense
NONRATT014284.2	1.536364828	−2.296449909	0.041338987	0.00591903	Up	Down	Antisense LncRNA
NONRATT026027.2	1.739597948	−2.007928677	0.043290989	0.042622989	Up	Down	Exonic_sense

Figure 5

Cluster analysis of 104 lncRNA therapeutic targets. The labels below the hierarchical clustering thermal map indicate sample number and the labels on the right indicate gene number. Each column indicates a sample and each row indicates a gene. Black indicates no change at the gene level, red indicates upregulation, and green indicates downregulation. The brightness of the color indicates an increase or decrease at the gene level. Genes with similar expression are clustered with samples.

Functional analysis

Based on the 104 selected lncRNAs, 623 relevant protein-coding genes were identified by cis and trans-regulation, and these were further investigated in GO and KEGG enrichment analysis, which suggested that the biological processes involved in the treatment of T2DM mainly contained the glycolipid metabolism and amino acid metabolism. These processes could be specialized as the cellular lipid metabolic process, cellular amino acid metabolic process, regulation of glucose metabolic process, fatty acid catabolic process, fat cell differentiation, lipid biosynthetic process, apoptotic mitochondrial changes, cholesterol metabolic process, glycosylation, and DNA repair (Figure 6A). The molecular function (MF) and the cell components (CC) of GO analysis also are shown in Figure 6A. In vivo, the biological functions were performed by the coordination of different genes. Based on the significant enrichment of a pathway, the major biochemical metabolic pathways and signal transduction pathways involved in target genes were determined. The KEGG results (Figure 6B) revealed that the peroxisome proliferator-activated receptor (PPAR) signaling pathways, amino acid metabolic pathways (tyrosine metabolism; glycine, serine, and threonine metabolism; tryptophan metabolism; beta-alanine metabolism; arginine and proline metabolism; valine, leucine, and isoleucine degradation; and lysine degradation), mammalian target of rapamycin (mTOR) signaling pathway, and the lipid related metabolism pathway (Adipocytokine signaling pathway, glycerolipid metabolism, sphingolipid metabolism, Adipocytokine signaling pathway, and fatty acid degradation), were mainly involved in the treatment of liraglutide for T2DM.

Figure 6

(A) GO enrichment analysis results. The ordinate indicates the specific GO entry name. The abscissa indicates the richness factor. The color of dots indicates the significance of GO (q value), and the shape of the colored dots indicates the affiliation of the categories in the GO database. The size indicates the number of genes mapped to this GO entry. (B) KEGG enrichment analysis results. The abscissa indicates the richness factor. The larger the richness factor, the greater the degree of enrichment. The ordinate is the name of the pathway entry. According to the ranking information of richness factor, the 20 most important KEGG pathways were displayed.

The validation of QRT-PCR

The Q-RTPCR results were approximately the same as the sequencing results, and further verified the accuracy of the sequencing results (Figure 7).

Figure 7

QRT-PCR results of 6 DE lncRNAs were consistent with the sequencing results. The abscissa is the name of the lncRNA and the ordinate is relative relationship.

Discussion

As a worldwide epidemic disease, T2DM has become one of the major diseases that threatens human health, with the characteristics of high incidence, high disability, and multiple complications [1]. Liraglutide is a GLP-1 analogue that can reduce endothelial endoplasmic reticulum stress and insulin resistance, reduce blood lipid levels and blood pressure, and protects against cardiovascular diseases [3-5]. In the present study, the lncRNA therapeutic targets of liraglutide in the treatment of diabetes were screened for the first time, and potentially involved biological processes and signaling pathways were explored.

HE staining and immunofluorescence analyses demonstrated that liraglutide can improve the morphology of islet cells and improve the function of β cells in diabetic model rats. We identified 104 lncRNA therapeutic targets of liraglutide against T2DM based on the lncRNA sequencing. Then, we randomly selected several of these 104 lncRNAs for verification based on qRT-PCR, which were consistent with the sequencing outcomes, confirming the accuracy of the sequencing results. The results of the present study lay a strong foundation for further exploration of the mechanism underlying the effect of GLP-1 drugs such as liraglutide in treatment of T2DM, and may improve the diagnosis and treatment of diabetes, which is of crucial clinical significance. Some of the 104 lncRNAs we found have been reported to be highly associated with the onset and treatment of diabetes. For example, Yu [13] indicated that oxidative stress in cardiomyocytes and apoptosis induced by high glucose can be regulated by the lncRNA NONRATT007560.2, suggesting that it could play a vital role in the occurrence and development of diabetic cardiomyopathy. Research on lncRNAs is still at the preliminary stage of exploration and basic research, and the present study may provide information useful for further research [14]. We plan to further explore the specific function of the 104 lncRNAs identified in the present study, as well as the upstream and downstream regulation mechanisms involved.

Our GO enrichment results demonstrated that the major biological processes involved in the effect of liraglutide consisted of glycolipid metabolism and amino acid metabolism. KEGG enrichment analysis also revealed that liraglutide exerts its anti-diabetic effect mainly via the PPAR signaling pathway, amino acid metabolic pathway, mTOR pathway, Wnt pathway, and lipid metabolism-related pathway. Our literature review found that the PPAR signaling pathway, the lipid metabolism-related pathway, the Wnt pathway, and the mTOR pathway are affected by liraglutide [15-17]. However, we found no studies indicating that liraglutide alleviates T2DM through the amino acid metabolic pathway. Hence, we think that our study may provide new research insights into the mechanism underlying the effect of this drug in the treatment of T2DM.

Specifically, as a nuclear hormone receptor activated by fatty acid and its derivatives, PPAR is a ligand-activated receptor in the nuclear hormone receptor family, which can be separated into 3 subtypes: PPARα, PPARβ, and PPARγ [18]. Recently, PPARs were found to be closely correlated with energy (lipid and sugar) metabolism, cell differentiation, proliferation, apoptosis, and inflammatory response, the biological effects of which are achieved through complex signaling pathways [19]. PPAR transcriptional activity is regulated by non-gene crosstalk with phosphatases and kinases, including erk1/2, p38-mapk, PKA, PKC, AMPK, and GSK3 [18]. Zhang [20] revealed that the heart function of diabetic rats could be improved by liraglutide via the PPAR signaling pathway.

mTOR is a serine/threonine protein kinase that is highly conserved in structure and function, and mainly exists in the form of mTORCl and mTORC2 complexes in vivo [21]. As a nutrient sensor, mTORC1 is located at the center of the complex signaling network, and is a regulatory protein for various crucial signaling pathways in cells [22]. The mTORC1 signal was reported to be highly associated with the functions of islet cells, which activated mTORC1 signaling in pancreatic β cells, resulting in the upregulation of insulin levels and glucose, stimulating insulin secretion (GSIS) [23]. Among various mTORC1 regulatory factors, leucine is the most efficient amino acid activating the mTORCI signaling pathway, which stimulates the phosphorylation of p70S6K via the mTOR pathway, thus promoting insulin secretion [20,23]. The KEGG results indicated that liraglutide regulates leucine metabolism. Evidence suggests that liraglutide controls the mTOR pathway by regulating leucine as one of the vital therapeutic mechanisms. Zhang [17] demonstrated that liraglutide can reverse myocardial damage by promoting autophagy via the AMPK-mTOR signaling pathway in the Zucker diabetic fatty rats. Suppression of the Wnt signaling pathway is a risk factor for the development of T2DM [24]. In islet cells, GLP-1 activates cAMP by acting on the GLP-1 receptors on the cell membrane, and cAMP acts as a second messenger, transmitting signals into the cell to activate the Wnt signaling pathway, promoting islet cell regeneration and reducing apoptosis [25].

Our literature review found no studies showing that liraglutide alleviates T2DM through the amino acid metabolic pathway [26]. In the present study, we found that liraglutide can play a role in the treatment of T2DM through the amino acid metabolic pathway. Therefore, we think that our study may provide new research insights into the mechanism underlying the effect of liraglutide against T2DM. This finding may be of great significance and warrants further exploration. The relationship between amino acid metabolism and T2DM has only recently received. There are numerous types of amino acids, and most of them are biogenic sugars and ketogenic amino acids, which can be used as the substrates of the tricarboxylic acid cycle to participate in gluconeogenesis and promote the production of endogenous glucose [27]. In recent years, studies have proved that in insulin-sensitive tissues such as skeletal muscle and liver and fat tissues, several amino acids can activate mTOR-S6K1, phosphorylate irs-1 Ser312 and ser636/639 residues, reduce the activity of insulin-induced downstream factors such as PI3k/Akt, lead to insulin resistance, and inhibit glucose transport [28,29]. It was reported that leucine promotes insulin secretion, inhibits AMPK, and activates the mTOR-S6K1 signaling pathway [30]. Additionally, leucine, phenylalanine, and arginine play roles in promoting insulin secretion [31,32]. Arginine can also suppress the production of vascular endothelial oxygen free radicals and reduce protein kinase activity and plasma triglyceride levels by providing nitric oxide [33]. In addition, the utilization of tissue glucose and the hypoglycemic mechanism is promoted by histidine and glycine, which reduces oxidative stress and chronic inflammation and improve insulin sensitivity [34,35]. In the present study, DE lncRNAs were found in the control group, the T2DM model group, and the liraglutide group in STZ-induced rats. We identified 104 lncRNAs that are therapeutic targets of liraglutide in the treatment of T2 DM, and assessed the major biological processes and the signaling pathways involved. Our lncRNA sequencing results provide a solid basis for understanding the mechanism involved in the effect of liraglutide treatment of diabetes from the perspective of epigenetics. We plan to further explore the function and upstream and downstream regulatory mechanisms of these therapeutic targets.

Conclusions

We explored the potential therapeutic targets and pathways of liraglutide against T2DM in streptozotocin-induced diabetic rats based on lncRNA sequencing. We found 104 lncRNA targets of liraglutide that affect T2DM, with 27 upregulated and 77 downregulated, including NONRATT030354.2, MSTRG.1456.6, and NONRATT011758.2. The major biological processes involved are glucose, lipid, and amino acid metabolism (P value <0.05). Liraglutide alleviates T2DM mainly through the following pathways: Wnt, PPAR, amino acid metabolism signaling, mTOR, and lipid metabolism-related pathways. Our results have clinical significance and may assist in future treatment of T2DM.