Competing endogenous RNA network analysis explores the key lncRNAs, miRNAs, and mRNAs in type 1 diabetes

Background Type 1 diabetes (T1D, named insulin-dependent diabetes) has a relatively rapid onset and significantly decreases life expectancy. This study is conducted to reveal the long non-coding RNA (lncRNA)-microRNA (miRNA)-mRNA regulatory axises implicated in T1D. Methods The gene expression profile under GSE55100 (GPL570 and GPL8786 datasets; including 12 T1D samples and 10 normal samples for each dataset) was extracted from Gene Expression Omnibus database. Using limma package, the differentially expressed mRNAs (DE-mRNAs), miRNAs (DE-miRNAs), and lncRNAs (DE-lncRNAs) between T1D and normal samples were analyzed. For the DE-mRNAs, the functional terms were enriched by DAVID tool, and the significant pathways were enriched using gene set enrichment analysis. The interactions among DE-lncRNAs, DE-miRNAs and DE-mRNAs were predicted using mirwalk and starbase. The lncRNA-miRNA-mRNA interaction network analysis was visualized by Cytoscape. The key genes in the interaction network were verified by quantitatively real-time PCR. Results In comparison to normal samples, 236 DE-mRNAs, 184 DE-lncRNAs, and 45 DE-miRNAs in T1D samples were identified. For the 236 DE-mRNAs, 16 Gene Ontology (GO)_biological process (BP) terms, four GO_cellular component (CC) terms, and 57 significant pathways were enriched. A network involving 36 DE-mRNAs, 8 DE- lncRNAs, and 15 DE-miRNAs was built, such as TRG-AS1—miR-23b/miR-423—PPM1L and GAS5—miR-320a/miR-23b/miR-423—SERPINA1 regulatory axises. Quantitatively real-time PCR successfully validated the expression levels of TRG-AS1- miR-23b -PPM1L and GAS5-miR-320a- SERPINA1. Conclusion TRG-AS1—miR-23b—PPM1L and GAS5—miR-320a—SERPINA1 regulatory axises might impact the pathogenesis of T1D.

other specific types, such as monogenic diabetes, unclassified diabetes and hyperglyacemia first detected during pregnancy [6,7]. The symptoms of T1D mainly include thirst, polydipsia, polyuria, polyphagia, fatigue, and rapid weight loss [8].An annual 3%-4% increase in the incidence of T1D in childhood was estimated in several developed countries [9]. Besides, T1D decreases approximately 11-13 years of life expectancy in developed countries, and even more time in developing countries [10]. Quality of life of T1D patients was significantly decreased than general population. Therefore, investigation in the pathogenesis of T1D is warranted to improve the prognosis of T1D patients.
Considerable efforts have been made on the molecular mechanism of T1D during the past decades. The roles of non-coding RNAs, including microRNA (miRNA), long non-coding RNAs (lncRNAs) in the regulation of T1D have only recently been recognized [11]. Through reducing C-C motif chemokine receptor 2 expression, miR-125a-5p limits the migration and function of regulatory T cells (Tregs) in the pancreas of T1D patients [12]. Activity enhancement of miR-181a promotes the expression of nuclear factor of activated T cells 5 (NFAT5) while suppresses the induction of forkhead box P3 + Tregs; therefore, miR-181a/NFAT5 axis may provide targets for limiting islet autoimmunity in T1D patients [13,14]. The sera level of miR-375 is decreased in T1D patients, indicating that the changed circulating level of miR-375 may be a valuable marker of inflammation and metabolic alterations in T1D [15]. The crosstalk between p38 mitogen-activated protein kinase signaling pathway and the lncRNA metastasis associated lung adenocarcinoma transcript 1 (MALAT1) is correlated with the endothelial cell function in diabetic rats, and MALAT1 inhibition may be a promising approach for anti-angiogenic treatment of diabetic microvascular complications [16]. Besides, lncRNAs MEG3, TUG1 and PVT1 were also being reported to contribute in the pathophysiology of T1D and T1D-associated complications [17][18][19]. However, only a small part of non-coding RNAs in T1D have been revealed, and more comprehensive understanding of T1D should be developed to facilitate design of preventive therapeutic modalities.
MiiRNAs can silence gene expression by binding mRNAs, while competing endogenous RNAs (ceRNAs), such as lncRNAs can bind miRNAs via miRNA response elements competitively to regulate gene expression [20,21]. In this study, the microarray dataset of T1D was searched from Gene Expression Omnibus (GEO) database. Afterwards, differential expression analysis, enrichment analysis, and lncRNA-miRNA-mRNA network (ceRNA network) analysis were executed to investigate the important regulatory axises involved in T1D.

Data source and data preprocessing
From GEO database (http://www.ncbi.nlm.nih.gov/geo/), the gene expression profile of T1D (accession number: GSE55100, including GPL570 and GPL8786 datasets) was downloaded. This dataset was deposited by Gu et al. [22]. There were 12 peripheral blood mononuclear cell (PBMC) samples from T1D patients and 10 PBMC samples from normal controls in both of GPL570 (including mRNA and lncRNA expression data) and GPL8786 (including miRNA expression data) datasets. The clinical characteristics of the 22 samples were listed in Table 1.
Using the Release 26 (GRCh38.p10) reference genome in GENCODE database (https ://www.genco degen es.org/ human /) [23], sequence alignment was conducted. Only the unique alignment sequences were remained. Based on corresponding GTF annotation files, mRNAs (with annotation information of "protein coding") and lncRNAs (with annotation information of "antisense", "sense_intronic", "lincRNA", "sense_overlapping", "pro-cessed_transcript", "3prime_overlapping_ncRNA", and "non_coding") were identified, respectively. Subsequently, the probes without relevant gene symbols were filtered out. For multiple probes that mapped to one gene symbol, their average value was used as the final expression value of this gene.

Differential expression analysis
For GPL570 and GPL8786 datasets, differential expression analysis between T1D and normal samples were performed using limma package (http://www.bioco nduct or.org/packa ges/2.9/bioc/html/limma .html) [24] in R. The corresponding p values of all genes were obtained through statistical test and were conducted with multiple test correction using Benjamini & Hochberg method [25]. The |log fold change (FC)|≥ 0.5 and adjusted p value < 0.05 were considered as the thresholds for obtaining differentially expressed mRNAs (DE-mRNAs). Meanwhile, the p value ≤ 0.05 was regarded as the cutoff criterion for selecting differentially expressed miR-NAs (DE-miRNAs) and differentially expressed lncRNAs (DE-lncRNAs).

Enrichment analysis for the DE-mRNAs
The Gene Ontology (GO)_biological process (BP), GO_molecular function (MF), and GO_cellular component (CC) functional terms involving the DE-mRNAs were analyzed using DAVID tool (version 6.7, https :// david -d.ncifc rf.gov/) [26], and the corresponding results were visualized by the R package GOplot (http://cran.rproje ct.org/web/packa ges/GOplo t) [27]. Using gene set enrichment analysis (GSEA) [28], Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment for the DE-mRNAs was carried out. The adjusted p value < 0.05 was used for screening the significant results.

Construction of lncRNA-miRNA-mRNA interaction network
Using corr.test (parameters: ci = F, adjust = "BH") in R package psych [29], the Pearson correlation coefficients [30] of the expression values of the DE-lncRNAs and the DE-mRNAs were calculated. The lncRNA-mRNA pairs with |r|≥ 0.7 and adjusted p value < 0.05 were selected and a lncRNA-mRNA co-expression network was developed By using Cytoscape software (http://www.cytos cape.org) [31]. Using starbase database (version 3.0, http://starb ase. sysu.edu.cn/) [32], lncRNA-miRNA interactions were predicted for the DE-lncRNAs in the lncRNA-mRNA co-expression network. Simultaneously, miRNA-mRNA interactions were predicted for the DE-mRNAs in the lncRNA-mRNA co-expression network by using with mirwalk database (version 3.0, http://mirwa lk.umm.uniheide lberg .de/) [33]. Finally, a lncRNA-miRNA-mRNA interaction network was visualized by Cytoscape software [31] by integrating the DE-lncRNAs and DE-mRNA regulated by the same DE-miRNA. Furthermore, GO and KEGG enrichment analysis for the mRNAs involved in the interaction network were conducted using DAVID tool [26].

Validation of key genes using quantitatively real-time PCR (qRT-PCR)
T1D patients were recruited from the department of Endocrinology, China-Japan Union Hospital of Jilin University. Age and sex matched healthy controls were recruited from volunteers of physical examination. The experiments were approved by Ethics Committee of China-Japan Union Hospital of Jilin University [No.(2020)linshen (20,201,127)]. Written informed consent was received from all participants. Total RNA were extracted from the PBMCs of T1D patients and normal control samples using RNAiso Plus (TaKaRa, Shiga, Japan). The total RNA was reversed into cDNA with pri-meScript RT Master Mix (TaKaRa) and amplified on ABI ViiA7 (ThermoFisher, USA). GAPDH was used as the internal reference of DE-lncRNAs and DE-mRNAs and U6 was regarded as the internal reference of DE-miR-NAs. The primers were listed in Table 2.

Statistical analysis
The data of qRT-PCR was processed by GraphPad Prism 5.0 (San Diego, CA, USA). Data were presented as mean ± standard deviation. Differences between two groups were determined by t-test. P < 0.05 was regarded as statistical significance level.

Differential expression analysis
A total of 1671 mRNAs, 1990 lncRNAs, and 533 miR-NAs were identified from the gene expression profile of GSE55100.
TRG-AS1 has been reported to be oncogenic in glioblastoma [34], hepatocellular carcinoma [35] and tongue squamous cell carcinoma [36]. It might act as a ceRNA to regulate miRNA-543-Yes-associated protein 1 in tongue squamous cell carcinoma [36], miR-4500-BACH1 in hepatocellular carcinoma [35] and miR-877-5p-SUZ12 in glioblastoma [34]. It also identified to be involved in repeated implantation failure [37]. However, its role in diabetes was not characterized    TUBA1C, TUBA1A, THBS1, C3, CORO1A, ATP6V1C1, RAB5C, ITGB1, HGS, TUBB1,  ATP6V1B2, CD36, THBS3, ATP6V0B, FCGR2B, MRC1, TUBB2A, CD14, ATP6V0C,  THBS4, CYBB, RAB7A, TLR2, CLEC7A, ATP6V1A, ATP6V0D1, NCF2, TLR6, NCF1 CCR10, IL3RA, IL17RA, TNFSF13, IL19, CCL8, EGFR, IL11RA, OSM, CXCL12, CXCR5,  IFNLR1, TNFRSF18, IL17RC, PRL, IFNA17, PLEKHO2, IL9R, IL17RB, IL1RAP, IL1A,  IL18, HGF, VEGFA, CCR6, IFNGR2, TNFSF10, PF4V1, CXCL8, CCL20, TNFSF8,  IL1R1, CSF2RA, CCL2, CCR9, CSF2RB, TNF, TNFRSF10C, CXCR2, CXCR1, CSF3R,  CXCL1, IL1R2, TNFRSF17,  previously. In our study, we predicted the miRNAs and mRNAs that might be regulated by TRG-AS1. The lncRNA-miRNA-mRNA network displayed that it could directly regulate 5 miRNAs including miR-423 and miR-23b, and multiple mRNAs, including PPM1L indirectly. MiR-23b/27b expression is decreased in the muscle stem cells of T2D patients, which exerts a promyogenic function via the p53 pathway [38]. Nuclear factor, erythroid 2/miR-423-5p axis can induce gluconeogenesis and hyperglycemia through inhibiting the family with sequence similarity 3 member A -adenosine triphosphate (ATP)-serine/threonine kinase Akt pathway, which is involved in the progression of T2D and nonalcoholic fatty liver disease [39]. Increased miR-320 impairs lipid metabolism and gluconeogenesis by targeting adiponectin receptor 1 (AdipoR1), and thus miR-320 may be taken as a possible target for T2D therapy [40]. The elevated levels of the T lymphocytes expressing gamma-delta T cell receptor are implicated in the islet autoimmune process in individuals at high risk of T1D, which may serve as a promising indicator for the development of T1D [41,42]. PPM1L mediates inositol-requiring protein-1 phosphorylation and endoplasmic reticulum stress signaling, which is considered as a causal gene for metabolic abnormalities [43]. The macrophage-enriched network has a causal correlation with metabolic disease traits, which involved three obesity genes (including PPM1L) [44]. Therefore, we hypothesized that TRG-AS1 and PPM1L might also play roles in the mechanisms of T1D through the TRG-AS1-miR-23b/miR-423-PPM1L regulatory axis. GAS5 is a member of 5′ terminal oligopyrimidine class which could regulate cell growth, proliferation, and survival [45]. Reduced serum levels of GAS5 are related to diabetes; therefore, serum GAS5 levels combined with other parameters may be used for identifying people at high risk of diabetes more accurately [46]. Shi et al. demonstrated that GAS5 regulates insulin signaling in adipocytes, and suggested it might be a potential target for T2DM [47]. GAS5 knockdown causes cell cycle arrest and impairs insulin synthesis and secretion in Min6 pancreatic β-cells, and thus GAS5 may function in maintaining the identity and function of β cells [48]. In this study, GAS5 was found decreased in T1D patients in both RNA-seq results and the qRT-PCR results, which is consistent with previous studies. We found GAS5 could regulate 7 miR-NAs, including miR-320a, miR-23b and miR-423, and tens of mRNAs, including SERPINA1 in T1D. MiR-320 negatively mediates the expression of fibronectin, endothelin 1, and vascular endothelial growth factor via extracellular signal-regulated kinases 1 and 2 in Fig. 2 The long non-coding RNA (lncRNA)-mRNA co-expression network. Red and yellow circles separately represent mRNAs and lncRNAs high glucose-treated human umbilical vein endothelial cells, which may provide a novel approach for treating diabetic complications [49]. Wei et al. found miR-320 mimic correlated with impaired gluconeogenesis and lipid metabolism by regulating adipoR1 [40]. SERPINA1 is a serine protease inhibitor that could target elastase, plasmin, thrombin, trypsin, chymotrypsin, and plasminogen activator. It was shown to be decreased in serum of obese mice and human subjects and the imbalance between SERPINA1 and neutrophil elastase contributed to insulin resistance [50]. Lower levels of SERPINA1 selectively impaired the ATP-binding cassette transporter A1 cholesterol efflux capacity in T2D [51]. Other members in this family were also reported to participate in development of diabetes. Anti-SERPINB13 antibody contributes to Reg gene expression and beta cell proliferation, and the immunological response may hinder the progression of T1D [52]. The serum concentrations of SERPINA12 (vaspin) is increased in T2D patients, which may be a candidate marker for evaluating the risk of severe macrovascular complications and the status of old T2D patients [53]. Thus, GAS5-miR-320a/miR-23b/miR-423-SERPINA1 regulatory axis might also function in the pathogenesis of T1D. Fig. 3 The long non-coding RNA (lncRNA)-miRNA-mRNA interaction network. Red, yellow, and blue circles represent mRNAs, lncRNAs, and miRNAs, respectively Table 5 The Gene Ontology (GO) functional terms and pathways enriched for the mRNAs involved in the long non-coding RNA (lncRNA)-miRNA-mRNA interaction network BP