Skip to main content

Identification of key regulatory genes and their working mechanisms in type 1 diabetes



Type 1 diabetes (T1D) is an autoimmune disease characterized by the destruction of beta cells in pancreatic islets. Identification of the key genes involved in T1D progression and their mechanisms of action may contribute to a better understanding of T1D.


The microarray profile of T1D-related gene expression was searched using the Gene Expression Omnibus (GEO) database. Then, the expression data of two messenger RNAs (mRNAs) were integrated for Weighted Gene Co-Expression Network Analysis (WGCNA) to generate candidate genes related to T1D. In parallel, T1D microRNA (miRNA) data were analyzed to screen for possible regulatory miRNAs and their target genes. An miRNA–mRNA regulatory network was then established to predict the key regulatory genes and their mechanisms.


A total of 24 modules (i.e., clusters/communities) were selected using WGCNA analysis, in which three modules were significantly associated with T1D. Further correlation analysis of the gene module revealed 926 differentially expressed genes (DEGs), of which 327 genes were correlated with T1D. Analysis of the miRNA microarray showed that 13 miRNAs had significant expression differences in T1D. An miRNA–mRNA network was established based on the prediction of miRNA target genes and the combined analysis of mRNA, in which the target genes of two miRNAs were found in T1D correlated genes.


An miRNA–mRNA network for T1D was established, based on which 2 miRNAs and 12 mRNAs were screened, suggesting that they may play key regulatory roles in the initiation and development of T1D.

Peer Review reports


The diagnosis of type 1 diabetes (T1D) is based on clinical manifestations of insulin-dependent diabetes, and the majority of cases are found in childhood or puberty with diabetic ketoacidosis [1]. Furthermore, data has shown that the incidence of T1D in children is still increasing [2, 3]. Traditionally, T1D is considered an immune disease caused by the destruction of pancreatic beta cells by T lymphocytes, leading to hyperglycemia and hypoglycemia-related complications [4]. Local islet inflammation can be partially attributed to the regulation ofbeta cells and infiltrating immune cells. In addition, increasing evidence has shown that islet inflammation can be ascribed to T1D candidate genes and environmental factors, including viral infection [5,6,7,8]. Research exploring the regulatory mechanism of T1D showed that miR-92a could regulate KLF2 expression to influence the function of B cells [9]. In addition, genes such as RNASEH1, BANK1, and SLC40A1 have also been reported to be involved in T1D regulation [10,11,12]. Nevertheless, data concerning key regulatory genes in T1D remain to be clarified, and more effort is required before a clearer understanding of T1D pathogenesis can be acquired. Therefore, screening for T1D regulatory genes may provide new directions for a better understanding of T1D pathology.

Weighted correlation network analysis (WGCNA) is a data mining (or computational) method used to describe correlation patterns among genes through gene co-expression networks in complicated diseases [13]. WGCNA has been used successfully to identify disease-related modules (clusters/communities) and genes. For instance, Chen et al. highlighted two modules and genes associated with bone mineral density through WGCNA [14, 15]. Riquelme et al. compared the data of T1D patients and healthy controls in one microarray using WGCNA and identified modules in co-expression gene networks that may be associated with T1D, in which the possible regulatory genes in T1D were identified [16]. The majority of previous studies mainly focused on single susceptibility genes or genome-wide association studies to explore the possible hereditary factors for T1D [17,18,19,20].

The importance of the immune system in the development of diabetes was reported in a previous study [21], and further evidence indicated the involvement of peripheral blood mononuclear cells (PBMCs) in the immunoregulation of T1D [22]. To broaden the data size, we integrated data from the T1D PBMC microarray to enable data excavation on T1D regulatory genes using T1D as a relative character. The integrated data were then subjected to WGCNA analysis to determine the possible regulatory genes. In a study by Roggli et al., some microRNAs (miRNAs), including miR-21 and miR-34a, were possibly involved in beta cell apoptosis in diabetes [23]. In addition, let-7c-5p and miR-25 have been reported to serve as biomarkers for T1D initiation and development [24,25,26,27]. Thisevidence supports the possible regulatory role of miRNAs in T1D. In this study, we used WGCNA analysis to analyze the T1D miRNA microarray to identify the key regulatory genes, and explored their possible working mechanisms in T1D after conjoint analysis on the T1D miRNA microarray. This study aimed to provide a potential therapeutic basis for a better understanding of the initiation and development of T1D.

Materials and methods

Data download and treatment

T1D microarray data were screened in the Gene Expression Omnibus (GEO) database ( based on the criteria that the sample size for both controls and T1D should be greater than 10, and the sample source should be PBMCs. Three independent T1D microarrays were obtained, including two messenger RNA (mRNA) microarrays and one miRNA microarray (Table 1) [28, 29]. Then, the two mRNA microarrays were integrated using the R limma package and sva package [30,31,32] and subjected to Bayesian adjustments. Genes co-expressed in both microarrays were saved during integration.

Table 1 Data downloaded from GEO database for analysis

WGCNA analysis

The adjusted microarray was analyzed using the R WGCNA package [13, 33]. The top 25% of the maximum expression variance was selected for WGCNA analysis, and 4895 genes remained for WGCNA analysis. The hclust function was used for outlier detection using “average”. The soft threshold was set to the value recommendedby the WGCNA analysis. A one-step network construction was used to construct the gene module. Then, the similarity matrix was calculated based on the eigengenes of each module using the R cor function, hclust function, and Pearson method. Clustering analysis was performed using “average” and the modules with high similarity were integrated during which 0.7 was set as the cut height. Then, the correlation analysis between the integrated module and trait was determined based on the eigengenes using the R cor function and corP valueStudent function. Subsequently, the genes in candidate modules were extracted, and based on eigengenes, the correlation between genes and modules was analyzed using the R cor function and corPvalueStudent function. Genes with a correlation coefficient better than 0.5 were selected for further analysis.

Differential analysis

Modules with significant correlation with T1D were selected, and their correlations were analyzed using the cor function with the module eigengene being used for calculation. Genes with a correlation coefficient above 0.5 were selected as candidate genes for subsequent differential analysis. The R limma package was used to analyze the expression of candidate genes in integrated data, and genes with p < 0.05 were defined as differentially expressed genes (DEGs). Differentially expressed miRNAs were searched in the miRNA microarray using the Limma package based on the criteria of |logFC| > 0.7 and p < 0.05.

Gene ontology (GO) and Kyoto Encyclopedia of genes and genomes (KEGG) enrichment analyses

GO and KEGG enrichment analyses were performed to detect biological functions and potential pathways of DEGs using the R clusterprofiler package with p < 0.05 considered as the cutoff value [34, 35].

miRNA–mRNA network

Funrich v3.1.3 software [36] was used to predict miRNA target genes. The data of target genes and T1D-related DEGs were analyzed for overlaps, in which the target genes with reversed miRNA expression patterns were obtained. The miRNA–mRNA network was constructed using cytoscapev3.9.1 [37].


Outlier detection

Two independent T1D microarrays (GSE156035 and GSE55098) were downloaded from the GEO database. The data in these two microarrays were integrated and adjusted, and 62 samples were identified, including 30 cases of healthy controls and 32 cases of T1D. There were 19,579 coexpressed genes in the two microarrays (Additional file 1: Table S1). Then, sample clusters were detected for outliers (Fig. 1), and the detection showed no obvious outliers; therefore, all samples were included for subsequent analysis.

Fig. 1
figure 1

Sample clustering to detect outliers using WGCNA analysis. Note: WGCNA, weighted correlation network analysis.

WGCNA analysis

Genes were analyzed using the WGCNA method, and modules related to T1D were obtained. To ascertain the obtained modules, they were screened using soft thresholding power (Fig. 2A, B). Screening showed that the scale-free parameter was approximately 0.9 in response to a soft thresholding power of 3. Therefore, the soft thresholding power was set to 3 for subsequent analysis. Genes were analyzed using a one-step method to distinguish them into 31 modules. Each module was marked in a different color (Fig. 2C), and the module in gray indicates non-significant genes.

Fig. 2
figure 2

WGCNA analysis. Note: A, B The scale free topology model fit under different soft thresholding power. The number herein is the relative soft thresholding power. When soft thresholding was set at 3, the approximate scale free topology can be obtained, indicating genes can be efficiently distinguished at the soft thresholding of 3. C Each branch of the dendrogram stands for a gene. Based on the topology overlapped cluster, gens were classified into different modules in different color. Each color stands for one module, in which highly connected genes were included. A total of 31 modules were identified. WGCNA, weighted correlation network analysis.

Integration of similar modules

A total of 31 modules were identified. The similarity among the 31 modules was calculated, and those with a similarity coefficient within 0.7 were integrated (Fig. 3A). A total of 24 modules were obtained, as shown in Fig. 3B, C. A total of 500 genes were randomly selected to visualize module membership (Fig. 3D).

Fig. 3
figure 3

Integration of similar modules. Note: A Correlation analysis on modules before integration. The vertical ordinate is the similarity of modules. Lower the ordinate value, higher the similarity. Red line is the cut height. Cut height at 0.7 was set in this study and modules with cut height under 0.7 were selected for integration. B Gene dendrogram. C Module dendrogram after integration. D Visualization on gene networks. About 500 genes in the gene network were randomly selected and the TOM value was described by heat map.

Correlation analysis among modules

The correlations among the 24 modules were analyzed (Fig. 4A), and the results showed no significant correlations among the modules, indicating the independence of the genome in these modules. Subsequently, the features of the T1D cases were compared with those of the 24 modules using correlation analysis (Fig. 4B). The results showed that the cyan module was negatively correlated with T1D, whereas the blue and light cyan modules were positively correlated with T1D. These findings suggest that the genes in these three modules may be of great importance in the initiation and development of T1D.

Fig. 4
figure 4

Correlation analysis among modules. Note: A Heat map describing the correlation among modules. Red indicates for positively correlated and blue for positively correlated. The horizontal and vertical axis indicates for different colored modules. The histogram in the right stands for color level. B The correlation between modules and T1D development. Numbers are the Pearson coefficient and the p value of correlation is in the brackets. The histogram in the right stands for color level.

Differential analysis of the candidate genes

Genes in the above three modules were extracted, and the correlation between clustered genes and modules was calculated (Additional file 2: Table S2, Additional file 3: Table S3, Additional file 4: Table S4). Genes correlated with module eigengenes (correlation coefficient > 0.5) were selected for subsequent analysis and 926 genes were obtained (Additional file 5: Table S5). The expression differences of these 926 genes in the integrated microarray were analyzed, and 327 genes in T1D samples showed significantly different expression compared to their expression in controls (Additional file 6: Table S6 and Fig. 5). These results further support the hypothesis that these 327 genes play key roles in T1D.

Fig. 5
figure 5

Heat map describing the expression difference of candidate genes. The horizontal axis is the sample code and vertical axis is the gene name. Dendrogram in the left is the cluster of gene expression levels. The histogram in the right stands for color level.

Enrichment analysis

The screened 327 genes were subjected to GO and KEGG pathway enrichment analysis [38,39,40]. GO enrichment analysis revealed that 327 genes were mainly distributed in functional categories, including HSV1 infection, cell split, and nuclear split (Fig. 6A). KEGG enrichment analysis showed that the cell cycle and MAPK signaling pathways were enriched (Fig. 6B). Studies have highlighted that T1D can lead to beta cell loss and influence the splitting of beta cells [41,42,43]. Additionally, the MAPK gene is involved in the initiation and development of T1D [44]. HSV1 infection has also been reported to influence T1D progression [45]. This evidence suggests that the cell cycle and MAPK signaling pathways enriched by KEGG may have an important regulatory role in T1D.

Fig. 6
figure 6

Functional enrichment analysis of candidate genes. Note: A GO enrichment analysis. Horizontal axis is GeneRatio and the vertical axis is category. The histogram in the right stands for color level. B KEGG pathway enrichment analysis. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

MiRNA differential analysis

One independent miRNA microarray, GSE55099, was downloaded from the GEO database, which included 10 control samples and 12 T1D samples. Differential analysis on GSE55099 found 13 miRNAs with significant differential expression (Fig. 7A, B), among which some miRNAs including miR-423, let-7i, and miR-25 were reported to regulate T1D [24, 27, 46, 47]. Among the 13 miRNAs, 5 were upregulated in T1D and 8 were downregulated in T1D.

Fig. 7
figure 7

Differential analysis on T1D related miRNA microarray. Note: A Volcano map describing differential expressed miRNA. Horizontal axis is -log10pvalue and the vertical axis is logFC. Red dot is the miRNA that is upregulated in T1D and green dot is the miRNA that is downregulated in T1D. B Heap map for differential expressed miRNA. T1D, type 1 diabetes.

Prediction of miRNA–mRNA network

The target genes of the 13 candidate miRNAs were predicted (Additional file 7: Table S7). The results of the differential expression of miRNAs in T1D and the prediction of their target genes were analyzed for overlaps with the 327 genes to screen the genes that had reversed expression with the candidate miRNAs. Subsequently, an miRNA–mRNA regulatory network was established. Among the 13 miRNAs, 2 miRNAs were identified with their target genes with significantly differential expression in T1D (Fig. 8). These results suggest that these two miRNAs may play key roles in regulating T1D by mediating different target genes.

Fig. 8
figure 8

miRNA–mRNA regulatory network. Triangle indicates the candidate miRNA and circle is the candidate mRNA. The line between triangle and circle indicating the possible regulatory role between miRNA and mRNA.


Multiple T1D microarrays were integrated to obtain the T1D-related gene expression data. The gene expression data were screened using WGCNA and differential analysis, and 327 genes were identified to have a close relationship with T1D. In parallel, one miRNA microarray related to T1D was downloaded for analysis and 13 candidate miRNAs were identified to be associated with T1D. Further prediction of the 13 miRNAs and correlation analysis with the 327 genes contributed to the miRNA–mRNA regulatory network. Two miRNAs and 12 mRNAs had significant differential expression in T1D. These results suggest that miRNAs and mRNAs may play certain roles in the initiation and development of T1D.


Network-based analysis is an effective approach for exploring the connection between genes and pathways in diseases. This study used WGCNA analysis to integrate the clinical features of T1D patients with genome profiles to classify genes into several modules. Next, the correlation of modules with T1D was conducted to screen the key genes involved in T1D initiation and development.

WGCNA analysis generated three modules that were closely associated with T1D, and the genes in the modules were extracted for correlation analysis. Finally, 926 candidate genes were identified, in which differential expression analysis of these candidate genes yielded 327 genes. GO enrichment of these genes showed that they were involved in the cell-splittingprocess. In T1D, the number of beta cells decreases with the progression of T1D [43]. Nakanishi et al. identified that beta cells were decreased by 90% in T1D cases compared with controls [48]. Similarly, failure to detect beta cells in T1D was reported by Butler [49], and Diedisheim et al. detected no beta cell mass in seven T1D cases with disease progression of 5–25 years [50]. A report on 54 cases of T1D by Campbell-Thompson also failed to detect beta cell mass in the majority of T1D samples [42]. GO analysis and the above reports suggest that split beta cells may be of great importance in the development of T1D. KEGG enrichment analysis further identified that the 327 genes were mainly involved in the MAPK signaling pathway. In agreement with our results, the MAPK signaling pathway was shown to have a close relationship with T1D [51, 52].

To locate the key genes involved in T1D progression, we screened an miRNA microarray in the GEO database, based on which differential analysis was conducted, and 13 miRNAs were obtained. Then, the target genes of these 13 miRNAs were predicted, and further screening of the expression pattern of miRNAs and target genes contributed to the establishment of the miRNA–mRNA network. Eventually, 2 miRNAs and 12 target genes were identified. Consistently, let-7i has been implicated in T1D regulation [27], and miR-25 is believed to be associated with the function of beta cells in T1D [24]. The above results suggest that the candidate miRNAs may directly regulate T1D progression; however, how they regulate T1D development has rarely been reported. In this study, a total of 12 genes were screened, including B3GNT7, IRS2, TGFBR3, PRDM1, DUSP1, SMAD7, TOB1, SIK1, CD69, JOSD1, FOSL2, and DDIT4 by WGCNA analysis, and may be regulated by let-7i and miR-25. Among these genes, DUSP1 and SMAD7 have been reported to be involved in the regulation of T1D [53,54,55,56,57]. In addition, IRS2 is believed to be involved in the functional regulation of INS-1 cells to mediate T1D development [58]. In addition, a close relationship between IRS2 and insulin-related pathways has been documented in previous studies [59, 60]. A Chinese study reported that SMAD7 was implicated in the regulatory effect of the compound Coptodis Rhizoma capsule in diabetic rats [61]. Despite these studies, some candidate genes have not been reported in T1D, such as SIK1, FOSL2, and DDIT4. However, these genes have been found to regulate type 2 diabetes [62,63,64,65,66], suggesting that these 12 candidate genes may have an important regulatory role in T1D. Comprehensive research on these genes will facilitate our understanding of the T1D pathogenesis.

In this study, we identified 2 miRNAs and 12 mRNAs that may have certain effects on the regulation of T1D. Based on this analysis, it was noted that most of the mRNAs and miRNAs were not reported in T1D, which provides new research directions for T1D pathogenesis. In addition, one of the highlights of this study is to conduct further analysis of a possible regulatory mechanism based on the candidate genes. Overall, through WGCNA analysis, this study was conducted to facilitate our understanding of genes and mechanisms related to T1D progression.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the Gene Expression Omnibus (GEO) repository (


  1. Powers AC. Type 1 diabetes mellitus: much progress, many opportunities. J Clin Invest. 2021;131(8):e142242.

    Article  Google Scholar 

  2. Patterson CC, Dahlquist GG, Gyurus E, Green A, Soltesz G, Group ES. Incidence trends for childhood type 1 diabetes in Europe during 1989–2003 and predicted new cases 2005–20: a multicentre prospective registration study. Lancet. 2009;373(9680):2027–33.

    Article  Google Scholar 

  3. Vehik K, Hamman RF, Lezotte D, Norris JM, Klingensmith G, Bloch C, Rewers M, Dabelea D. Increasing incidence of type 1 diabetes in 0- to 17-year-old Colorado youth. Diabetes Care. 2007;30(3):503–9.

    Article  Google Scholar 

  4. American Diabetes Association Professional Practice C. 2. Classification and diagnosis of diabetes: standards of medical care in diabetes-2022. Diabetes Care. 2022;45(Suppl 1):S17–38.

    Article  Google Scholar 

  5. Santin I, Eizirik DL. Candidate genes for type 1 diabetes modulate pancreatic islet inflammation and beta-cell apoptosis. Diabetes Obes Metab. 2013;15(Suppl 3):71–81.

    Article  CAS  Google Scholar 

  6. Dotta F, Censini S, van Halteren AG, Marselli L, Masini M, Dionisi S, Mosca F, Boggi U, Muda AO, Del Prato S, et al. Coxsackie B4 virus infection of beta cells and natural killer cell insulitis in recent-onset type 1 diabetic patients. Proc Natl Acad Sci U S A. 2007;104(12):5115–20.

    Article  CAS  Google Scholar 

  7. Gonzalez-Moro I, Olazagoitia-Garmendia A, Colli ML, Cobo-Vuilleumier N, Postler TS, Marselli L, Marchetti P, Ghosh S, Gauthier BR, Eizirik DL, et al. The T1D-associated lncRNA Lnc13 modulates human pancreatic beta cell inflammation by allele-specific stabilization of STAT1 mRNA. Proc Natl Acad Sci U S A. 2020;117(16):9022–31.

    Article  CAS  Google Scholar 

  8. Op de Beeck A, Eizirik DL. Viral infections in type 1 diabetes mellitus: why the beta cells? Nat Rev Endocrinol. 2016;12(5):263–73.

    Article  CAS  Google Scholar 

  9. Wang W, Wang J, Yan M, Jiang J, Bian A. MiRNA-92a protects pancreatic B-cell function by targeting KLF2 in diabetes mellitus. Biochem Biophys Res Commun. 2018;500(3):577–82.

    Article  CAS  Google Scholar 

  10. Hao L, Mi J, Song L, Guo Y, Li Y, Yin Y, Zhang C. SLC40A1 mediates ferroptosis and cognitive dysfunction in type 1 diabetes. Neuroscience. 2021;463:216–26.

    Article  CAS  Google Scholar 

  11. Pineda-Trujillo N, Rodriguez-Acevedo A, Rodriguez A, Ruiz-Linares A, Bedoya G, Rivera A, Alfaro JM. RNASEH1 gene variants are associated with autoimmune type 1 diabetes in Colombia. J Endocrinol Invest. 2018;41(7):755–64.

    Article  CAS  Google Scholar 

  12. Zouidi F, Stayoussef M, Bouzid D, Fourati H, Abida O, Joao C, Ayed MB, Fakhfakh R, Thouraya K, Monjia H, et al. Association of BANK1 and cytokine gene polymorphisms with type 1 diabetes in Tunisia. Gene. 2014;536(2):296–301.

    Article  CAS  Google Scholar 

  13. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.

    Article  Google Scholar 

  14. Chen YC, Guo YF, He H, Lin X, Wang XF, Zhou R, Li WT, Pan DY, Shen J, Deng HW. Integrative analysis of genomics and transcriptome data to identify potential functional genes of BMDs in females. J Bone Miner Res. 2016;31(5):1041–9.

    Article  CAS  Google Scholar 

  15. He H, Zhang L, Li J, Wang YP, Zhang JG, Shen J, Guo YF, Deng HW. Integrative analysis of GWASs, human protein interaction, and gene expression identified gene modules associated with BMDs. J Clin Endocrinol Metab. 2014;99(11):E2392-2399.

    Article  CAS  Google Scholar 

  16. Riquelme Medina I, Lubovac-Pilav Z. Gene Co-expression network analysis for identifying modules and functionally enriched pathways in type 1 diabetes. PLoS ONE. 2016;11(6):e0156006.

    Article  Google Scholar 

  17. Hakonarson H, Grant SF. Genome-wide association studies in type 1 diabetes, inflammatory bowel disease and other immune-mediated disorders. Semin Immunol. 2009;21(6):355–62.

    Article  CAS  Google Scholar 

  18. Cooper JD, Smyth DJ, Smiles AM, Plagnol V, Walker NM, Allen JE, Downes K, Barrett JC, Healy BC, Mychaleckyj JC, et al. Meta-analysis of genome-wide association study data identifies additional type 1 diabetes risk loci. Nat Genet. 2008;40(12):1399–401.

    Article  CAS  Google Scholar 

  19. Luo DF, Buzzetti R, Rotter JI, Maclaren NK, Raffel LJ, Nistico L, Giovannini C, Pozzilli P, Thomson G, She JX. Confirmation of three susceptibility genes to insulin-dependent diabetes mellitus: IDDM4, IDDM5 and IDDM8. Hum Mol Genet. 1996;5(5):693–8.

    Article  CAS  Google Scholar 

  20. Bennett ST, Lucassen AM, Gough SC, Powell EE, Undlien DE, Pritchard LE, Merriman ME, Kawaguchi Y, Dronsfield MJ, Pociot F, et al. Susceptibility to human type 1 diabetes at IDDM2 is determined by tandem repeat variation at the insulin gene minisatellite locus. Nat Genet. 1995;9(3):284–92.

    Article  CAS  Google Scholar 

  21. Shalev A. Minireview: thioredoxin-interacting protein: regulation and function in the pancreatic beta-cell. Mol Endocrinol. 2014;28(8):1211–20.

    Article  Google Scholar 

  22. Axelsson S, Faresjo M, Hedman M, Ludvigsson J, Casas R. Cryopreserved peripheral blood mononuclear cells are suitable for the assessment of immunological markers in type 1 diabetic children. Cryobiology. 2008;57(3):201–8.

    Article  CAS  Google Scholar 

  23. Roggli E, Gattesco S, Caille D, Briet C, Boitard C, Meda P, Regazzi R. Changes in microRNA expression contribute to pancreatic beta-cell dysfunction in prediabetic NOD mice. Diabetes. 2012;61(7):1742–51.

    Article  CAS  Google Scholar 

  24. Nielsen LB, Wang C, Sorensen K, Bang-Berthelsen CH, Hansen L, Andersen ML, Hougaard P, Juul A, Zhang CY, Pociot F, et al. Circulating levels of microRNA from children with newly diagnosed type 1 diabetes and healthy controls: evidence that miR-25 associates to residual beta-cell function and glycaemic control during disease progression. Exp Diabetes Res. 2012;2012:896362.

    Google Scholar 

  25. Dotta F, Ventriglia G, Snowhite IV, Pugliese A. MicroRNAs: markers of beta-cell stress and autoimmunity. Curr Opin Endocrinol Diabetes Obes. 2018;25(4):237–45.

    Article  CAS  Google Scholar 

  26. Salas-Perez F, Codner E, Valencia E, Pizarro C, Carrasco E, Perez-Bravo F. MicroRNAs miR-21a and miR-93 are down regulated in peripheral blood mononuclear cells (PBMCs) from patients with type 1 diabetes. Immunobiology. 2013;218(5):733–7.

    Article  CAS  Google Scholar 

  27. Garavelli S, Bruzzaniti S, Tagliabue E, Prattichizzo F, Di Silvestre D, Perna F, La Sala L, Ceriello A, Mozzillo E, Fattorusso V, et al. Blood co-circulating extracellular microRNAs and immune cell subsets associate with type 1 diabetes severity. Int J Mol Sci. 2020;21(2):477.

    Article  CAS  Google Scholar 

  28. Yang M, Ye L, Wang B, Gao J, Liu R, Hong J, Wang W, Gu W, Ning G. Decreased miR-146 expression in peripheral blood mononuclear cells is correlated with ongoing islet autoimmunity in type 1 diabetes patients 1miR-146. J Diabetes. 2015;7(2):158–65.

    Article  CAS  Google Scholar 

  29. Santos AS, Cunha-Neto E, Gonfinetti NV, Bertonha FB, Brochet P, Bergon A, Moreira-Filho CA, Chevillard C, da Silva MER. Prevalence of inflammatory pathways over immuno-tolerance in peripheral blood mononuclear cells of recent-onset type 1 diabetes. Front Immunol. 2021;12:765264.

    Article  CAS  Google Scholar 

  30. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7): e47.

    Article  Google Scholar 

  31. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–3.

    Article  CAS  Google Scholar 

  32. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118–27.

    Article  Google Scholar 

  33. Langfelder P, Horvath S. Fast R functions for robust correlations and hierarchical clustering. J Stat Softw. 2012;46(11):i11.

    Article  Google Scholar 

  34. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141.

    CAS  Google Scholar 

  35. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  Google Scholar 

  36. Pathan M, Keerthikumar S, Ang CS, Gangoda L, Quek CY, Williamson NA, Mouradov D, Sieber OM, Simpson RJ, Salim A, et al. FunRich: an open access standalone functional enrichment and interaction network analysis tool. Proteomics. 2015;15(15):2597–601.

    Article  CAS  Google Scholar 

  37. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  Google Scholar 

  38. Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  Google Scholar 

  39. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545–51.

    Article  CAS  Google Scholar 

  40. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  Google Scholar 

  41. Lam CJ, Jacobson DR, Rankin MM, Cox AR, Kushner JA. Beta cells persist in T1D pancreata without evidence of ongoing beta-cell turnover or neogenesis. J Clin Endocrinol Metab. 2017;102(8):2647–59.

    Article  Google Scholar 

  42. Campbell-Thompson M, Fu A, Kaddis JS, Wasserfall C, Schatz DA, Pugliese A, Atkinson MA. Insulitis and beta-cell mass in the natural history of type 1 diabetes. Diabetes. 2016;65(3):719–31.

    Article  CAS  Google Scholar 

  43. Butler AE, Galasso R, Meier JJ, Basu R, Rizza RA, Butler PC. Modestly increased beta cell apoptosis but no increased beta cell replication in recent-onset type 1 diabetic patients who died of diabetic ketoacidosis. Diabetologia. 2007;50(11):2323–31.

    Article  CAS  Google Scholar 

  44. Long D, Chen Y, Wu H, Zhao M, Lu Q. Clinical significance and immunobiology of IL-21 in autoimmunity. J Autoimmun. 2019;99:1–14.

    Article  CAS  Google Scholar 

  45. Liao J, Jijon HB, Kim IR, Goel G, Doan A, Sokol H, Bauer H, Herrmann BG, Lassen KG, Xavier RJ. An image-based genetic assay identifies genes in T1D susceptibility loci controlling cellular antiviral immunity in mouse. PLoS ONE. 2014;9(9):e108777.

    Article  Google Scholar 

  46. Li C, Wei B, Zhao J. Competing endogenous RNA network analysis explores the key lncRNAs, miRNAs, and mRNAs in type 1 diabetes. BMC Med Genom. 2021;14(1):35.

    Article  CAS  Google Scholar 

  47. Liu Y, Ma M, Yu J, Ping F, Zhang H, Li W, Xu L, Li Y. Decreased serum microRNA-21, microRNA-25, microRNA-146a, and microRNA-181a in autoimmune diabetes: potential biomarkers for diagnosis and possible involvement in pathogenesis. Int J Endocrinol. 2019;2019:8406438.

    Article  Google Scholar 

  48. Nakanishi K, Kobayashi T, Miyashita H, Okubo M, Sugimoto T, Murase T, Kosaka K, Hara M. Relationships among residual beta cells, exocrine pancreas, and islet cell antibodies in insulin-dependent diabetes mellitus. Metabolism. 1993;42(2):196–203.

    Article  CAS  Google Scholar 

  49. Meier JJ, Bhushan A, Butler AE, Rizza RA, Butler PC. Sustained beta cell apoptosis in patients with long-standing type 1 diabetes: indirect evidence for islet regeneration? Diabetologia. 2005;48(11):2221–8.

    Article  CAS  Google Scholar 

  50. Diedisheim M, Mallone R, Boitard C, Larger E. Beta-cell mass in nondiabetic autoantibody-positive subjects: an analysis based on the network for pancreatic organ donors database. J Clin Endocrinol Metab. 2016;101(4):1390–7.

    Article  CAS  Google Scholar 

  51. Erener S, Marwaha A, Tan R, Panagiotopoulos C, Kieffer TJ. Profiling of circulating microRNAs in children with recent onset of type 1 diabetes. JCI Insight. 2017;2(4):e89656.

    Article  Google Scholar 

  52. Khan S, Jena GB. Protective role of sodium butyrate, a HDAC inhibitor on beta-cell proliferation, function and glucose homeostasis through modulation of p38/ERK MAPK and apoptotic pathways: study in juvenile diabetic rat. Chem Biol Interact. 2014;213:1–12.

    Article  CAS  Google Scholar 

  53. Yuan H, Xu J, Xu X, Gao T, Wang Y, Fan Y, Hu J, Shao Y, Zhao B, Li H, et al. Calhex231 alleviates high glucose-induced myocardial fibrosis via inhibiting itch-ubiquitin proteasome pathway in vitro. Biol Pharm Bull. 2019;42(8):1337–44.

    Article  CAS  Google Scholar 

  54. Dias HF, Kuhtreiber WM, Nelson KJ, Ng NC, Zheng H, Faustman DL. Epigenetic changes related to glucose metabolism in type 1 diabetes after BCG vaccinations: a vital role for KDM2B. Vaccine. 2022;40(11):1540–54.

    Article  CAS  Google Scholar 

  55. Jin Y, Sharma A, Carey C, Hopkins D, Wang X, Robertson DG, Bode B, Anderson SW, Reed JC, Steed RD, et al. The expression of inflammatory genes is upregulated in peripheral blood of patients with type 1 diabetes. Diabetes Care. 2013;36(9):2794–802.

    Article  CAS  Google Scholar 

  56. Barrett JC, Clayton DG, Concannon P, Akolkar B, Cooper JD, Erlich HA, Julier C, Morahan G, Nerup J, Nierras C, et al. Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes. Nat Genet. 2009;41(6):703–7.

    Article  CAS  Google Scholar 

  57. Takahashi P, Xavier DJ, Lima J, Evangelista AF, Collares CVA, Foss-Freitas MC, Rassi DM, Donadi EA, Passos GA, Sakamoto-Hojo ET. Transcript expression profiles and microRNA regulation indicate an upregulation of processes linked to oxidative stress, DNA repair, cell death, and inflammation in type 1 diabetes mellitus patients. J Diabetes Res. 2022;2022:3511329.

    Article  Google Scholar 

  58. Shen H, Sun J, Liu J, Wang L, Dong L. miR-181d promotes pancreatic beta cell dysfunction by targeting IRS2 in gestational diabetes mellitus. Ginekol Pol. 2021;92(8):563–70.

    Article  Google Scholar 

  59. White MF. IRS2 integrates insulin/IGF1 signalling with metabolism, neurodegeneration and longevity. Diabetes Obes Metab. 2014;16(Suppl 1):4–15.

    Article  CAS  Google Scholar 

  60. Valverde AM, Gonzalez-Rodriguez A. IRS2 and PTP1B: two opposite modulators of hepatic insulin signalling. Arch Physiol Biochem. 2011;117(3):105–15.

    Article  CAS  Google Scholar 

  61. Liu S, Chen XQ, Tang LQ, Yu N, Zhang XL, Du HF. Regulatory effect of compound Coptidis Rhizoma capsule on unbalanced expression of renal tissue TGF-beta1/BMP-7 and Smad signaling pathway in rats with early diabetic nephropathy. Zhongguo Zhong Yao Za Zhi. 2015;40(5):938–45.

    Google Scholar 

  62. Liao Q, Xu W, Luo Q, Wen X. Zhenqing recipe relieves diabetic nephropathy through the SIK1/SREBP-1c axis in type 2 diabetic rats. Am J Transl Res. 2021;13(12):13776–83.

    CAS  Google Scholar 

  63. Pan S, Li M, Yu H, Xie Z, Li X, Duan X, Huang G, Zhou Z. microRNA-143-3p contributes to inflammatory reactions by targeting FOSL2 in PBMCs from patients with autoimmune diabetes mellitus. Acta Diabetol. 2021;58(1):63–72.

    Article  CAS  Google Scholar 

  64. Wang C, Song D, Fu J, Wen X. SIK1 regulates CRTC2-mediated gluconeogenesis signaling pathway in human and mouse liver cells. Front Endocrinol (Lausanne). 2020;11:580.

    Article  Google Scholar 

  65. Li J, Li S, Hu Y, Cao G, Wang S, Rai P, Wang X, Sun K. The Expression level of mRNA, Protein, and DNA methylation status of FOSL2 of Uyghur in XinJiang in Type 2 diabetes. J Diabetes Res. 2016;2016:5957404.

    Article  Google Scholar 

  66. Wang H, Wang J, Qu H, Wei H, Ji B, Yang Z, Wu J, He Q, Luo Y, Liu D, et al. In vitro and in vivo inhibition of mTOR by 1,25-dihydroxyvitamin D3 to improve early diabetic nephropathy via the DDIT4/TSC2/mTOR pathway. Endocrine. 2016;54(2):348–59.

    Article  CAS  Google Scholar 

Download references


Not applicable.


Not applicable.

Author information

Authors and Affiliations



HL and XT contributed to the conception of the study; XH, JL and WJ performed the experiment; XH, JL, WJ and LW contributed significantly to analysis and manuscript preparation; HL, XH and XT performed the data analyses and wrote the manuscript; JL, WJ, LW helped perform the analysis with constructive discussions. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Xin Tan.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Two independent microarrays related to T1D (GSE156035 and GSE55098) were integrated for analysis and 19579 co-expressed genes were identified in both microarrays.

Additional file 2.

The genes in the module 1 were extracted and the correlation of clustered genes with modules was calculated.

Additional file 3.

The genes in the module 2 were extracted and the correlation of clustered genes with modules was calculated.

Additional file 4.

The genes in the module 3 were extracted and the correlation of clustered genes with modules was calculated.

Additional file 5.

Genes with a correlation coefficient above 0.5 were selected for subsequent analysis and 925 candidate genes were obtained.

Additional file 6.

The 925 candidate genes were used for differential analysis and 327 genes were found to have significant differences.

Additional file 7.

Prediction on the target genes of 13 candidate miRNAs (hsa-mir-28, hsa-mir-744, hsa-let-7i, hsa-mir-584, hsa-mir-423, hsa-mir-1249, hsa-mir-1275, hsa-mir-25, hsa-mir-146b, hsa-mir-345, hsa-mir-1237 and hsa-mir-1301.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, H., Hu, X., Li, J. et al. Identification of key regulatory genes and their working mechanisms in type 1 diabetes. BMC Med Genomics 16, 8 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Type 1 diabetes
  • GEO
  • Differentially expressed genes