Bioinformatic Reconstruction and Analysis of Gene Networks Related to Glucose Variability in Diabetes and Its Complications

Glucose variability (GV) has been recognized recently as a promoter of complications and therapeutic targets in diabetes. The aim of this study was to reconstruct and analyze gene networks related to GV in diabetes and its complications. For network analysis, we used the ANDSystem that provides automatic network reconstruction and analysis based on text mining. The network of GV consisted of 37 genes/proteins associated with both hyperglycemia and hypoglycemia. Cardiovascular system, pancreas, adipose and muscle tissues, gastrointestinal tract, and kidney were recognized as the loci with the highest expression of GV-related genes. According to Gene Ontology enrichment analysis, these genes are associated with insulin secretion, glucose metabolism, glycogen biosynthesis, gluconeogenesis, MAPK and JAK-STAT cascades, protein kinase B signaling, cell proliferation, nitric oxide biosynthesis, etc. GV-related genes were found to occupy central positions in the networks of diabetes complications (cardiovascular disease, diabetic nephropathy, retinopathy, and neuropathy) and were associated with response to hypoxia. Gene prioritization analysis identified new gene candidates (THBS1, FN1, HSP90AA1, EGFR, MAPK1, STAT3, TP53, EGF, GSK3B, and PTEN) potentially involved in GV. The results expand the understanding of the molecular mechanisms of the GV phenomenon in diabetes and provide molecular markers and therapeutic targets for future research.


Introduction
The concept of blood glucose variability (GV) attracts increasing attention in both basic and clinical science. Maintenance of GV within normal values is provided by a complex regulatory molecular-genetic network, which remains not fully understood. Many patients with diabetes have increased GV, which might be a major therapeutic challenge [1,2]. The highest GV is observed in insulin-treated type 1 and type 2 diabetic subjects. However, the initial increase in GV can be revealed in persons with impaired glucose tolerance [3].
Thus, the network associated with GV includes a variety of molecular components. Hormones constitute a significant part of this network. This is consistent with the well-known fact that maintenance of glycemic fluctuations within the physiological limits is controlled by fine and well-balanced hormonal regulation. Circadian rhythms in secretion of growth hormone, melatonin, prolactin, cortisol, leptin, and ghrelin are highly correlated with sleep and can contribute to the intra-day glycemic fluctuations [35][36][37].
Among non-hormonal factors, interleukin-6 (IL6) showed the highest number of connections in the GV-related network. Circadian rhythm in IL-6 secretion with an elevation in day-time hours was described [38]. As IL-6 modulates glucose disposal and insulin sensitivity [39], one could speculate it can affect intraday GV. On the other hand, excessive glucose fluctuations can increase IL-6 production. In cultured human endothelial cells, an intermittent increase in glucose levels caused a more pronounced increase in IL6 expression than a consistently high glucose [40]. The association between GV and plasma level of IL-6 has been shown in individuals with metabolic syndrome [41].
In patients with type 1 diabetes, an increase in the IL-6 levels was observed after a two-hour episode of hypoglycemia [42]. Hyperglycemia, following hypoglycemia, caused a further increase in IL-6 concentration [18]. The genes associated with insulin secretion and signaling turned out to be important nodes of the GV network. In humans, mutations in these genes (GCK, HNF1A, ABCC8, and KCNJ11) are causes of monogenic forms of diabetes known as maturity-onset diabetes of the young [43]. An increased GV was reported in patients with diabetes caused by mutations in the HNF1A [44] and ABCC8 [45]. Mutations in the insulin receptor gene (INSR) are associated with insulin resistance and hyperglycaemia and can lead to hyperinsulinemic hypoglycaemia in adults and children [46,47]. The GCK, KCNJ11, ABCC8, and INSR were identified as genetic variants associated with type 2 diabetes and/or glycemic traits in genome-wide association studies [48]. However, the contribution of the polymorphisms of these genes to GV in healthy individuals and patients with diabetes has not been studied yet. The role of a number of molecules, identified as the network components, in GV phenomenon remains to be clarified.
At the next step, we performed the comparisons of expression patterns of GV-related genes in different tissues with the use of Bgee [49]. The greatest expression of GV-related genes was observed in the cardiovascular system, followed by pancreas, adipose and muscle tissues, gastrointestinal tract, and kidney; retina and nerves demonstrated a less maximal expression score (Table S4). This fact gives further support to notion that GV could be the unifying link between impairment of insulin synthesis and signaling, glucose metabolism, and diabetic vascular complications. Table 2 demonstrates the results of the Gene Ontology (GO) enrichment analysis by DAVID service [50], aimed at revealing the overrepresented biological processes for GV-related genes. In addition to the regulation of insulin secretion and different aspects of glucose metabolism, this analysis revealed some signaling pathways considered as the cornerstones in regulation of cellular metabolism (mitogen-activated protein kinase (MAPK) cascade, protein kinase B signaling, JAK-STAT cascade, regulation of peptidyl-tyrosine phosphorylation), cell cycle (regulation of mitotic nuclear division and cell proliferation), and cell-cell signaling. Positive regulation of nitric oxide biosynthetic process was identified among overrepresented processes. It turned out that the biological processes identified as overrepresented are closely related to oxidative stress, glycation, endothelial dysfunction, and other processes that are considered as important players in the development of diabetes complications. In particular, oxidative stress can directly alter glucose homeostasis [51] and it itself is triggered by hyperglycemia [52]. The intensity of glycation is proportional to glucose concentration and depends on glucose metabolism [53,54]. MAPK cascade and nitric oxide biosynthesis are known as principal pathways in the biology of endothelial cells [55,56].
These data support the notion that GV-related genes could play an important role in pathogenesis of diabetic complications. We tested this assumption at the next stage of our work.

The Role of GV-Related Genes in the Networks of Diabetes Complications
The ANDSystem was used to construct gene networks associated with diabetic complications: cardiovascular disease ( Figure S3, Table S5), diabetic nephropathy ( Figure S4, Table S6), diabetic retinopathy ( Figure S5, Table S7), and diabetic neuropathy ( Figure S6, Table S8). For each of these gene networks, the role of GV-related genes was assessed by analyzing the betweenness centrality coefficients of the network participants (Table 3). Betweenness centrality reflects the number of shortest pathways between all pairs of genes in the analyzed gene network that go through a given gene. It describes the functional importance of a gene in a gene network.
Analysis shows that GV-related genes are significantly overrepresented in all networks of diabetes complications. In addition, the centrality of the GV-related genes was, on average, significantly higher than the average centrality of all participants of the considered gene networks (Table 3). These data support the key role of GV-related genes in the molecular orchestra of diabetes complications.
In order to find out what biological processes involved in pathogenesis of diabetic complications are triggered by glycemic fluctuations, the following analysis was carried out. Among the participants of the gene networks of diabetes complications, we selected those that were directly regulated by GV-related genes, according to ANDSystem. Thus, we identified 141 genes in the network of cardiovascular diseases, 219 genes in the network of diabetic nephropathy, 132 genes in the network of diabetic retinopathy, and 44 genes in the network of diabetic neuropathy. Further, using the DAVID service [50], the GO biological processes overrepresented for these genes were identified (Table S9). Significance of difference between the coefficient of betweenness centrality of genes associated with GV and all genes in the network It was found that the response to hypoxia was overrepresented in all four sets of genes, and in all cases, it was among the top ten most significantly overrepresented biological processes. This suggests an important role of hypoxia as a mediator between GV and diabetes complications. Diabetes is characterized by altered oxygen deprivation signaling in the targeted organs, including the vasculature and kidney. Cellular response to hypoxia is controlled by hypoxia-inducible factor (HIF), a transcriptional regulator. The HIF family is considered as a master regulator of endothelial biology with rising interest in the field of diabetes-driven atherosclerosis [57]. It is postulated that a shift in the balance of HIF-1α and HIF-2α promotes proinflammatory and profibrotic pathways in glomerular and renal tubular cells [58]. A growing body of research indicates HIF-1α and its target genes as the contributors to the retinal neovascularization in diabetic retinopathy [59]. In obesity and type 2 diabetes, inflammation and hypoxia of adipose tissue may play an important role in the development of insulin resistance [60]. In its turn, intermittent hypoxia, which is used to improve the adaptation to hypoxia, reduced blood glucose in patients with type 2 diabetes [61]. These data are in agreement with our findings indicating the response to hypoxia as an important mechanism linking GV and chronic diabetic complications.
Other 15 biological processes were revealed to be overrepresented simultaneously for three sets of genes important for cardiovascular disease, diabetic nephropathy, and diabetic retinopathy (Table 4). They are involved in the inflammatory response, platelet function, angiogenesis, regulation of blood pressure, and nitric oxide biosynthesis. These processes, along with endothelial dysfunction, are considered as the cornerstones in the pathogenesis of diabetic vascular complications [62]. Table 4. GO biological processes overrepresented (p-values with FDR correction < 0.05) for genes that are regulated by GV-related genes and participate in the gene networks of cardiovascular disease, diabetic nephropathy, and diabetic retinopathy.

Gene Ontology Biological Processes
Gene Ontology ID Molecular pathways, implicated in the regulation of cell cycle, proliferation, and cell-cell signaling, have been recognized among significantly overrepresented biological processes related to GV and complications. Among them, extracellular signal-regulated kinase (ERK) 1 and 2 cascade should be mentioned. ERK cascade plays a crucial role in multiple cellular processes, such as cell proliferation, differentiation, adhesion, migration, and survival [63]. Glucose can activate ERK1/2 in beta-cells, contributing to an increase in mitogenesis [64]. Moreover, activation of ERK is a component of insulin signaling. Accordingly, the targeting of the ERK cascade is considered as a potential treatment for insulin resistance and type 2 diabetes [65]. The activation of ERK1/2 is involved in pathogenesis of cardiac hypertrophy and dysfunction; recent studies showed the role of ERK in development of diabetic cardiomyopathy [66]. High glucose increases ERK phosphorylation in mesangial cells, promoting the synthesis of extracellular matrix, a key event in diabetic kidney disease [67].
Thus, it was revealed that GV-related genes occupy central positions in the gene networks of diabetes complications. The response to hypoxia, low-grade inflammation, abnormalities of angiogenesis and platelets function, disturbances in the regulation of blood pressure and vasodilatation, as well as alterations of cellular signaling pathways, could link GV with the development of diabetic complications.

Identification of New Candidate Genes in the GV Network
At the first step of this study, we identified by the ANDSystem the set of genes that are directly associated with both hyperglycemia and hypoglycemia, attributed them to GV. The Pathway Wizard tool of the ANDSystem allows one to identify indirect interactions between objects in the gene networks that pass through genes-intermediaries. Using this opportunity, it was possible to find genes/proteins for which there is no evidence on a direct connection with hypoglycemia and hyperglycemia, but there are indications of indirect relationships with glucose abnormalities through genes-intermediaries.
As genes-intermediaries, we considered the GV-related genes described in the first section of this article. The genes/proteins that were associated with genes-intermediaries, we considered as candidate genes potentially involved in GV network. The analysis identified 334 such candidate genes. Further, in order to rank candidates and determine the most promising ones, gene prioritization was carried out using the ToppGene system [68]. The GV-related genes were used as a training set.
The THBS1 gene received the highest priority. This gene encodes an adhesive glycoprotein thrombospondin 1 (THSP1), which regulates cell-cell and cell-matrix interactions. THSP1 is an important player in the regulation of adhesion, mobility, and growth of endothelial cells; it has important regulatory function in angiogenesis, inflammation, and tissue remodeling also [69]. Patients with diabetes and coronary artery disease have increased THSP1 levels in plasma [70]. The enhanced expression of THBS1 was observed in mice with streptozotocin-induced diabetes. In cultured fibroblasts, high glucose increased expression of THBS1 [71]. Being involved in atherosclerosis, inflammation, ischemia/reperfusion injury, cardiac hypertrophy, and heart failure, THBS1 plays a significant role in the development of cardiovascular diseases [72]. In diabetic complications, this molecule is considered to be a key inductor of fibrosis [73]. Figure 2 shows the network of interactions between the THBS1 and GV-related genes. It was demonstrated that THBS1 is able to bind to CD36, a transporter and sensor of free fatty acids. Binding of CD36 to THBS1 initiates anti-angiogenic signals, induction of endothelial cell apoptosis, and the inhibition of angiogenesis [74]. THBS1 differentially regulates the release of IL-6 and IL-10 by human monocytes [75]. Both hypoglycemia and hyperglycemia increased serum concentrations of IL-6 [76,77], while hyperglycemia down-regulated IL-10 mRNA levels [78]. THBS1 modulates epidermal growth factor 1 (EGR1) function [79]. At the same time, up-regulation of THBS1 is mediated by EGR1 [80]. Hyperglycemia is responsible for EGR1 activation, acetylation, and its prolonged expression [81].
Fibronectin is a component of the extracellular matrix, accumulating in the retinal and glomerular basement membranes and mesangial matrix in diabetic retinopathy and diabetic kidney disease [82,83]. High ambient glucose increases fibronectin production in human mesangial cells [84] and retinal pigment epithelial cells [85]. In mesangial cells, the effect was mediated via THSP and TGF-beta [84].
Heat shock protein 90α (Hsp90α), encoded by the HSP90AA1 gene, is the stress-inducible isoform of the molecular chaperone Hsp90. The expression of HSP90AA1 increases in some forms of cancer [86]. Previously, HSP90AA1 was recognized by bioinformatic analysis as a node in a protein-protein interaction network in diabetic nephropathy [87]. However, the exact role of Hsp90α in diabetes and its complications remains to be clarified.
Both EGF and EGFR, encoding epidermal growth factor (EGF) and its receptor respectively, were present in the list of gene-intermediates. In diabetes, EGF/EGFR system was discussed primarily as a promoter of wound healing [88]. Recently it was demonstrated that EGF signal is involved in diabetes-induced vascular dysfunction, remodeling, and transcriptome dysregulation associated with renal involvement [89]. In patients with type 2 diabetes, soluble EGFR correlates with fasting blood glucose, insulin resistance, and glycated hemoglobin HbA1c [90]. Interestingly, EGFR can regulate the transcription of proopiomelanocortin [91], a hypothalamic appetite regulator and ACTH precursor modulated by glucose level [92]. Both EGFR and TP53 regulate expression of IGFBP3 [93,94], which level is reduced in response to hyperglycemia [95]. Glycoprotein, a component of extracellular matrix, is involved in the cell adhesion, wound healing, blood coagulation, and tumor metastasis.

HSP90AA1
Heat shock protein 90 alpha family class A member 1 Chaperone promoting proper folding of target proteins during cell stress.

EGFR Epidermal growth factor receptor
Receptor for members of the epidermal growth factor family. It is located on the cell surface and promotes cell proliferation.

MAPK1
Mitogen-activated protein kinase 1 MAP kinase that acts as a trigger for a variety of biochemical signals, being involved in cell proliferation, differentiation, and development.

STAT3
Signal transducer and activator of transcription 3 Signal transducer and activator of transcription in response to cytokines and growth factors. It mediates the expression of many genes in response to cellular stimuli and plays a key role in cell growth, apoptosis, and immune processes.

TP53
Tumor protein p53 Tumor suppressor protein that responds to various cellular stresses and regulates the expression of target genes. It causes a stop of the cell cycle, apoptosis, DNA repair, and metabolic changes. Mutations in this gene are associated with different human cancers.

EGF Epidermal growth factor
Epidermal growth factor, which acts as a mitogenic factor by binding to a cell surface epidermal growth factor receptor. It plays an important role in the growth, proliferation, and differentiation of many cell types.

GSK3B
Glycogen synthase kinase 3 beta Serine-threonine kinase, a member of the glycogen synthase kinase subfamily, regulates glucose homeostasis and is involved in energy metabolism, inflammation, mitochondrial dysfunction, and apoptosis.
[80]. Hyperglycemia is responsible for EGR1 activation, acetylation, and its prolonged expression [81]. Fibronectin is a component of the extracellular matrix, accumulating in the retinal and glomerular basement membranes and mesangial matrix in diabetic retinopathy and diabetic kidney disease [82,83]. High ambient glucose increases fibronectin production in human mesangial cells [84] and retinal pigment epithelial cells [85]. In mesangial cells, the effect was mediated via THSP and TGF-beta [84].
Heat shock protein 90α (Hsp90α), encoded by the HSP90AA1 gene, is the stress-inducible isoform of the molecular chaperone Hsp90. The expression of HSP90AA1 increases in some forms of cancer [86]. Previously, HSP90AA1 was recognized by bioinformatic analysis as a node in a proteinprotein interaction network in diabetic nephropathy [87]. However, the exact role of Hsp90α in diabetes and its complications remains to be clarified.
Both EGF and EGFR, encoding epidermal growth factor (EGF) and its receptor respectively, were present in the list of gene-intermediates. In diabetes, EGF/EGFR system was discussed primarily as a promoter of wound healing [88]. Recently it was demonstrated that EGF signal is involved in diabetes-induced vascular dysfunction, remodeling, and transcriptome dysregulation associated with renal involvement [89]. In patients with type 2 diabetes, soluble EGFR correlates with fasting blood glucose, insulin resistance, and glycated hemoglobin HbA1c [90]. Interestingly, EGFR can regulate the transcription of proopiomelanocortin [91], a hypothalamic appetite regulator and ACTH precursor modulated by glucose level [92]. Both EGFR and TP53 regulate expression of IGFBP3 [93,94], which level is reduced in response to hyperglycemia [95].
Mitogen-activated protein kinase 1 (MAPK1), also known as ERK2, is a molecular target of metformin, a principal antihyperglycemic agent [96]. MAPK1 is one of the essential molecules implicated in canonical insulin signaling, an important modulator of protein synthesis [96]. In the liver, MAPK1 is involved in regulation of glucose and lipid metabolism required for physiological metabolic adaptation [97]. In mice on a high-fat/high-sucrose diet, MAPK1 knockout promoted insulin resistance and impairment of glucose tolerance [98].
Signal transducer and activator of transcription 3 (STAT3) is a principal modulator of β-cell survival and function [99]. It is involved in the regulation of hepatic insulin sensitivity and gluconeogenesis [100]. The excessive STAT3 signaling contributes to skeletal muscle insulin resistance [101]. Recent data indicated the role of STAT3 signaling in high-glucose-induced podocyte injury [102], tubular cell dysfunction [103], and diabetic retinopathy [104].
Tumor protein p53 (TP53), a tumor suppressor, was recently identified as a metabolic regulator. In the liver, TP53 improves insulin sensitivity via inhibition of MAPKs and NF-κB pathways [105]. Mitogen-activated protein kinase 1 (MAPK1), also known as ERK2, is a molecular target of metformin, a principal antihyperglycemic agent [96]. MAPK1 is one of the essential molecules implicated in canonical insulin signaling, an important modulator of protein synthesis [96]. In the liver, MAPK1 is involved in regulation of glucose and lipid metabolism required for physiological metabolic adaptation [97]. In mice on a high-fat/high-sucrose diet, MAPK1 knockout promoted insulin resistance and impairment of glucose tolerance [98].
Signal transducer and activator of transcription 3 (STAT3) is a principal modulator of β-cell survival and function [99]. It is involved in the regulation of hepatic insulin sensitivity and gluconeogenesis [100]. The excessive STAT3 signaling contributes to skeletal muscle insulin resistance [101]. Recent data indicated the role of STAT3 signaling in high-glucose-induced podocyte injury [102], tubular cell dysfunction [103], and diabetic retinopathy [104].
Tumor protein p53 (TP53), a tumor suppressor, was recently identified as a metabolic regulator. In the liver, TP53 improves insulin sensitivity via inhibition of MAPKs and NF-κB pathways [105]. The protein regulates the function of glucose transporters affecting their transcription and translocation. It negatively regulates glycolysis and positively regulates gluconeogenesis [106]. In the beta cells, TP53 mediates gluco-and lipotoxicity and induces apoptosis [106]. In patients with type 2 diabetes, the level of TP53 increased with the age, duration of diabetes, and waist-to-hip ratio [107].
Glycogen synthase kinase 3 beta (GSK3B) is a serine/threonine protein kinase that is involved in glycogen metabolism, regulation of gene transcription, cytoskeleton organization, cell cycle, and apoptosis. The enzyme inhibits glycogen synthase and decreases glycogen storage in skeletal muscles and, to a lesser extent, in the liver; it also deteriorates glucose tolerance and insulin sensitivity. Insulin inactivates GSK3B via protein kinase B (Akt) [108]. It was demonstrated that GSK3B modulates inflammatory response in the islets [109]. Women with gestational diabetes had increased GSK3B activity in skeletal muscle and adipose tissue, where the enzyme promoted the inflammatory response [110]. Blood glucose fluctuations suppressed the expression of phosphorylated GSK3B and contributed to cardiac fibrosis in a model of diabetic cardiomyopathy [111].
Phosphatase and tensin homolog (PTEN) is essential for both cellular growth and metabolism. PTEN counteracts insulin signaling via PI3K/Akt pathway and reduces insulin sensitivity in mouse adipose tissue, liver, and β-cells [112]. In humans, mutations in PTEN gene increase the risks of obesity and cancer, but decrease the risk of type 2 diabetes [113]. Recent experimental studies indicate the role of PTEN in β-cell function and apoptosis [99], hyperglycemia-induced angiogenesis suppression [114], and development of diabetic nephropathy [115].
As an indirect confirmation of the relevance of performed gene prioritization, we checked the co-occurrence of the candidate gene together with the terms denoting GV in full-text PubMed Central articles. We used information from full-text articles to test our prioritization because ANDSystem stores the information gathered from the abstracts, but not full-text articles. For each candidate gene, using the hypergeometric distribution, the likelihood of meeting the name of this gene, together with the terms denoting GV in full-text of PubMed Central articles, was estimated (Table S11). At the next step, the Spearman's rank correlation coefficient was calculated between the prioritization rank of a gene and the probability to discover the gene in articles together with the terms describing GV by chance. The correlation coefficient was found to be 0.3 and the statistical significance of the correlation was 1.6 × 10 −8 . This indirectly confirms the adequacy of the gene prioritization. It means that the top candidate genes are promising for experimental verification of their relationship with GV.
For further verification of obtained results, we used transcriptomic data stored in GEO (https://www.ncbi.nlm.nih.gov/geo/info/), the largest open repository of experiments in the field of gene expression. We were able to find only one study (GEO ID: GSE40779) related to glucose fluctuations in diabetes. In this study, carried out in Rattus norvegicus, non-coding RNA profile was assessed in the heart under glucose fluctuations compared to sustained hyperglycemia [116]. It was found that the expression of 20 of 350 scanned miRNAs significantly differed between diabetic hearts and diabetic hearts with glucose fluctuations [116]. We have conducted the search for the target genes of these 20 differentially expressed miRNAs in the miRTarBase database [117]. For 16 differentially expressed miRNAs, 69 target genes in Rattus norvegicus have been found. Four miRNAs (rno-miR-136, rno-miR-296*, rno-miR-326, and rno-miR-30c-2*) had no target genes in miRTarBase. At the next step, human orthologous genes were found for these 69 target genes using the NCBI Orthologs service (Table S12). Next, we matched these genes with those we identified as the candidate genes most likely associated with GV.
It turned out that five genes, namely PTEN, TNFRSF1A, CCND1, TLR2, and RELA, which we proposed as the candidate genes for experimental verification, are targets of microRNAs differentially expressed in the GSE40779 experiment. According to the hypergeometric distribution, the probability of observing this for random reasons is low (p-value = 0.006). In addition, all of these five genes were among the top 100 candidate genes with highest priority predicted by us, which, according to the hypergeometric distribution, is unlikely by chance (p-value = 0.002). Moreover, three of five discussed genes were among the top 50 candidate genes with highest priority, which is also unlikely to be observed for random reasons (p-value = 0.03), and the PTEN was among the top ten most-priority candidates. This indirectly confirms the results of our bioinformatics analysis.

Study Limitations and Future Remarks
Our study is not without limitations. Taking into account the lack of the data on the molecular aspects of GV itself, we identified GV-related genes as those associated with both hyperglycemia and hypoglycemia. Consequently, some molecular aspects of physiological glucose fluctuations might not be recognized. Future reconstruction and analysis of tissue-specific gene networks may also be important. As the ANDSystem operates with PubMed data only, some relevant information could have been missed. Our study is a hypothesis-generating one. The role of newly identified candidate genes, as well as some molecular pathways and biological processes, needs further experimental verification.

Materials and Methods
Associative gene networks were built using the ANDSystem version: 20.0413b646 (ICG SB RAS, Novosibirsk, Russia) [33,34,118]. The computer tool ANDSystem was developed for automatic reconstruction of associative molecular-genetic networks describing the mechanisms of diseases and other complex phenotypic traits. It utilizes text mining, an automatic knowledge extraction from the texts of scientific publications. Previously, the ANDSystem was used to identify candidate genes associated with comorbidity of asthma and tuberculosis [119], preeclampsia, diabetes and obesity [120], primary open-angle glaucoma [121], asthma and hypertension [122], and lymphedema [123].
Comparative analysis of the expression of GV-related genes in the tissues was carried out by the "Expression comparison" function (https://bgee.org/?page=expression_comparison) provided on the website of Bgee resource [49].
Gene Ontology enrichment analysis was performed using the DAVID service version: 6.8 (LHRI, Frederick, MD, USA) [50] with the following parameters: organism-"Homo sapiens", Gene_Ontology-"GOTERM_BP_DIRECT". Gene Ontology biological processes were considered as overrepresented if the p-values with false discovery rate (FDR) correction were less than 0.05.
The calculation of betweenness centrality coefficients for the gene network participants was made by the Statistics function in the Analysis section of the ANDSystem [35].
The identification of candidate genes potentially involved in the regulation of GV was carried out using the Pathway Wizard tool of ANDSystem [118]. The parameter "Pathway length" was set equal to 4. "Object 1" was hypoglycemia, "Object 2", "Object 3" and "Object 4"-any human gene/protein, "Object 5"-hyperglycemia. For the interactions "Object 1-"Object 2" and "Object 4"-"Object 5" the following "Interaction property" was set: "Direction"-any; all databases in the "Database filter" section; all types of interactions in the section "Interaction type filter". For the interactions "Object 2"-"Object 3" the following "Interaction property" was set: "Direction"-"Right to left"; all databases in the "Database filter" section; all types of interactions in the section "Interaction type filter", except "association", "coexpression", "expression", "interaction", "involvement". For the interactions "Object 3"-"Object 4" the following "Interaction property" was set: "Direction"-"Left to right"; all databases in the "Database filter" section; all types of interactions in the "Interaction type filter" section, except for "association", "coexpression", "expression", "interaction", "involvement". Genes presented at the "Object 3" were considered as the candidate genes.
The prioritization of candidate genes was carried out using the ToppGene system (BMI CCHMC, Cincinnati, OH, USA) [68]. A list of GV-related genes was used as a training set. The candidate gene list was used as a test set. The parameters in the "Training parameters" section have been set: "all Feature".
The information about the function of the candidate gene products was obtained from the NCBI Gene database (https://www.ncbi.nlm.nih.gov/gene/).
The evaluation of the co-occurrence of the candidate gene together with the terms denoting GV in full-text PubMed Central articles was carried out as follows. All articles mentioning GV were found in PubMed Central interface (https://www.ncbi.nlm.nih.gov/pmc/) by the query: "glucose variability" OR "glycemic variability" OR "glycemia variability" OR "glucose fluctuations" OR "glycemic fluctuations" OR "sugar fluctuations" OR "glucose fluctuation" OR "glycemic fluctuation" OR "sugar fluctuation" OR "glucose excursions" OR "glycemic excursions" OR "unstable blood glucose" OR "fluctuating blood glucose". For each candidate gene, all articles in PubMed Central mentioning this gene were counted. All articles in which the candidate gene name and a list of terms denoting GV occurred simultaneously were counted. Further, using these numbers, statistical significance of the co-occurrence of names of genes-candidates and GV was assessed. For this, a hypergeometric distribution, implemented in the hypergeom.sf function of the scipy.stats package [124] of the Python programming language, was used. The false discovery rate (FDR) correction (https://www.sdmproject.com/utilities/?show=FDR) for the multiple comparisons was carried out using the "Signed Differential Mapping" tool [125].
The correlation between the rank of a candidate gene and the probability of meeting the name of this gene together with terms denoting glycemic variability in full-text articles of PubMed Central was estimated using Spearman's rank correlation coefficient. For this, the stats.spearmanr function of the scipy.stats package [124] of the Python programming language was utilized.
Search for Series of experiments of gene expression was performed in GEO (https://www. ncbi.nlm.nih.gov/geo/info/) by query: "a phrase denoting GV" AND (diabetic OR diabetes). The following phrases one by one have been used as "a phrase denoting GV": "glucose variability", "glycemic variability", "glycemia variability", "glucose fluctuations", "glycemic fluctuations", "sugar fluctuations", "glucose fluctuation", "glycemic fluctuation", "sugar fluctuation", "glucose excursions", "glycaemia excursions", "unstable blood glucose", "fluctuating blood glucose". The open experimental data used in our work [116] have the GEO ID: GSE40779 and are available by the link https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE40779. The target genes of differentially expressed miRNAs were extracted from the miRTarBase database [117]. Human genes that are orthologous of Rattus norvegicus genes, targeted by differentially expressed miRNAs, were found using the NCBI Orthologs service (https://www.ncbi.nlm.nih.gov/kis/info/how-are-orthologs-calculated/). Assessment of the statistical significance of the enrichment of candidate-gene list by genes, targeted by differentially expressed miRNAs, was performed using a hypergeometric distribution. For this purpose, we used the hypergeom.sf function of the scipy.stats package [124] of the Python programming language.

Conclusions
The identification of the molecular effects of GV in diabetes and its complications is a new challenge for Endocrinology and Molecular Medicine. In this study, we have performed for the first time the reconstruction and analysis of the GV-related network, using the bioinformatics analysis and text mining of publications indexed in PubMed. Firstly, we had constructed the network of GV, which included 37 genes/proteins related to both hyperglycemia and hypoglycemia. In this network, we have identified the genes of hormones regulating glucose metabolism, proinflammatory cytokines, growth factors, regulators of insulin secretion and signaling, transcriptional factors, and other molecules. The cardiovascular system, pancreas, adipose and muscle tissues, gastrointestinal tract, and kidney were recognized as the loci with the highest expression of GV-related genes. We have identified regulation of insulin secretion, glucose homeostasis, as well as some signaling pathways, which are involved in the regulation of cellular metabolism (MAPK cascade, protein kinase B signaling, JAK-STAT cascade, regulation of peptidyl-tyrosine phosphorylation), cell cycle (regulation of mitotic nuclear division and cell proliferation), and cell-cell signaling among biological processes that are most actively regulated by GV-related genes. We have shown that GV-related genes are central to the networks of diabetic complications, including cardiovascular disease, diabetic nephropathy, retinopathy, and neuropathy. Finally, we have identified new candidate genes (THBS1, FN1, HSP90AA1, EGFR, MAPK1, STAT3, TP53, EGF, GSK3B, and PTEN), promising for experimental verification of their role in the GV phenomenon. The obtained results provide further possibilities for the identification of new molecular markers and therapeutic targets in diabetes and its complications.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/22/ 8691/s1. Figure S1: View of the gene network of hyperglycemia reconstructed using the ANDSystem. Figure  S2: View of the gene network of hypoglycemia reconstructed using the ANDSystem. Figure S3: View of the gene network, associated with cardiovascular disease, reconstructed using the ANDSystem. Figure S4: View of the gene network, associated with diabetic nephropathy, reconstructed using the ANDSystem. Figure S5: View of the gene network, associated with diabetic retinopathy, reconstructed using the ANDSystem. Figure S6: View of the gene network, associated with diabetic neuropathy, reconstructed using the ANDSystem. Table S1: Interactions between participants of the gene network associated with hyperglycemia. Table S2: Interactions between participants of the gene network associated with hypoglycemia. Table S3: gene network associated with glycemic variability. Table S4: Bgee expression analysis. Table S5: Interactions between participants of the gene networks associated with cardiovascular disease. Table S6: Interactions between participants of the gene networks associated with diabetic nephropathy. Table S7: Interactions between participants of the gene networks associated with diabetic retinopathy. Table S8: Interactions between participants of the gene networks associated with diabetic neuropathy. Table S9: GO enrichment analysis for glycemic variability genes found in diabetes complications. Table S10: ToppGene prioritization results. Table S11: co-occurrence of candidates with glycemic variability in PubMed Central. Table S12: analysis of differentially expressed miRNAs from GSE40779 experiment.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

FDR
False discovery rate EGF Epidermal growth factor ERK Extracellular signal-regulated kinase GO Gene ontology GSK3B Glycogen synthase kinase 3 beta GV Glucose variability HIF Hypoxia-inducible factor Hsp90α Heat shock protein 90α ID Identification number MAPK Mitogen-activated protein kinase PTEN Phosphatase and tensin homolog miRNAs microRNAs STAT3 Signal transducer and activator of transcription 3 THSP1 Thrombospondin 1