Prioritization of genes involved in endothelial cell apoptosis by their implication in lymphedema using an analysis of associative gene networks with ANDSystem

Background Currently, more than 150 million people worldwide suffer from lymphedema. It is a chronic progressive disease characterized by high-protein edema of various parts of the body due to defects in lymphatic drainage. Molecular-genetic mechanisms of the disease are still poorly understood. Beginning of a clinical manifestation of primary lymphedema in middle age and the development of secondary lymphedema after treatment of breast cancer can be genetically determined. Disruption of endothelial cell apoptosis can be considered as one of the factors contributing to the development of lymphedema. However, a study of the relationship between genes associated with lymphedema and genes involved in endothelial apoptosis, in the associative gene network was not previously conducted. Methods In the current work, we used well-known methods (ToppGene and Endeavour), as well as methods previously developed by us, to prioritize genes involved in endothelial apoptosis and to find potential participants of molecular-genetic mechanisms of lymphedema among them. Original methods of prioritization took into account the overrepresented Gene Ontology biological processes, the centrality of vertices in the associative gene network, describing the interactions of endothelial apoptosis genes with genes associated with lymphedema, and the association of the analyzed genes with diseases that are comorbid to lymphedema. Results An assessment of the quality of prioritization was performed using criteria, which involved an analysis of the enrichment of the top-most priority genes by genes, which are known to have simultaneous interactions with lymphedema and endothelial cell apoptosis, as well as by genes differentially expressed in murine model of lymphedema. In particular, among genes involved in endothelial apoptosis, KDR, TNF, TEK, BMPR2, SERPINE1, IL10, CD40LG, CCL2, FASLG and ABL1 had the highest priority. The identified priority genes can be considered as candidates for genotyping in the studies involving the search for associations with lymphedema. Conclusions Analysis of interactions of these genes in the associative gene network of lymphedema can improve understanding of mechanisms of interaction between endothelial apoptosis and lymphangiogenesis, and shed light on the role of disturbance of these processes in the development of edema, chronic inflammation and connective tissue transformation during the progression of the disease. Electronic supplementary material The online version of this article (10.1186/s12920-019-0492-9) contains supplementary material, which is available to authorized users.


Background
Lymphedema is a chronic progressive disease, resulting in a significant loss of productivity, which affects more than 150 million people worldwide [1]. At the heart of the pathogenesis of lymphedema lies the appearance of chronic high-protein edema, characterized by the pathological accumulation of intercellular fluid with largemolecule proteins in the interstitial space due to a defect in lymphatic drainage caused by congenital malformation (primary lymphedema), lymphatic obstruction or by the destruction of lymphatic vessels (secondary lymphedema). Long existence of the high-protein edema causes chronic inflammation leading to the replacement of adipose tissue with connective tissue, increasing the volume of connective tissue matrix, which subsequently leads not only to an increase in size of body parts, but also to a secondary disturbance of lymphatic transport and drainage [2]. Most often it affects the lower extremities, but can also affect the upper limbs, face, trunk, external genitalia, etc. [3]. Primary lymphedema is a disease caused by dysfunction of lymphatic vessels, their aplasia, dysplasia and hypoplasia. Disturbance of drainage function of lymphatic vessels leads to accumulation of fluid rich in proteins under the epidermis. Most often, edema is localized on the lower limbs. Family forms of primary lymphedema are most often an autosomal dominant disease, but autosomal recessive types of inheritance of this pathology are also found [4]. It is believed that developmental defects and the resulting dysfunction of the lymphatic system are the cause of primary lymphedema, as well as associated syndromes, and can be genetically determined [5]. Secondary lymphedema is most often the result of a complex treatment of breast cancer, including removal of the axillary lymph nodes and/or radiation therapy [6]. Recently a number of associations between gene polymorphisms, including cytokine genes [7], and the development of breast cancer-related lymphedema were identified [8,9]. One of the possible causes of clinical manifestation of primary lymphedema in middle age and the development of secondary lymphedema after treatment of breast cancer may be apoptosis of the endothelium, whose role is discussed in publications on the pathogenesis of lymphedema [10,11].
Apoptosis is a form of cell death, characterized by a number of morphological and molecular features, including exposure of phosphatidylserine on the cell membrane, blebbing of the plasma membrane, cell shrinkage, cytoskeleton rearrangement, nuclear collapse, chromatin condensation, DNA fragmentation and formation of "apoptotic bodies" [12]. A number of works discuss the role of endothelial cell apoptosis in the pathogenesis of lymphedema [10,11,13,14]. However, the association between genes involved in endothelial apoptosis and genes associated with lymphedema has not been yet analyzed.
One of the most effective and frequently used approaches for the identification of the potential associations between genes and diseases is a prioritization [15]. Nowadays, there are many known methods for the prioritization of genes. A first, and the largest group of such methods is based on the analysis of genomic and transcriptomic data, as well as data on homologues [16]. The second class of prioritization methods involves approaches based on the analysis of gene networks and protein-protein interactions networks. An example of such methods is GUILD [17]. Integrated methods that combine approaches from the first and the second groups represent the third group. The well-known examples of such systems from the third group are Topp-Gene [15] and Endeavour [18].
Previously, we developed criteria for the prioritization of genes, by using the well-known systems and also considering the structure of associative gene networks from ANDSystem [19,20]. ANDSystem is a computer tool, designed for the automated extraction of knowledge from the texts of scientific publications and automatic reconstruction of the associative gene networks by using the retrieved information, describing the mechanisms of diseases, as well as other complex phenotypic traits. The knowledge base of ANDSystem contains over 30 million facts describing genetic regulation, gene associations with diseases, protein-protein interactions, catalytic reactions, transport pathways, etc., extracted from more than 25 million PubMed abstracts [21,22]. In particular, ANDSystem was used for the identification of candidate genes associated with comorbidity of preeclampsia, diabetes and obesity [23], asthma and tuberculosis [24], as well as asthma and hypertension [19].
In this study, our prioritization criteria were applied to identify the endothelial apoptosis-related genes, potentially involved in lymphedema. Three approaches were used to assess the quality of prioritization, based on the evaluation of the enrichment of the list of top priority genes by genes for which there was indirect evidence confirming their association with lymphedema: (1) enrichment with genes shared between lymphedema and endothelial apoptosis; (2) enrichment with genes, the names of which are significantly co-occurring in full-text articles with the key word «lymphedema»; (3) enrichment with genes differentially expressed in the murine model of lymphedema. All these quality criteria showed significant enrichment.
Among the genes with the highest priority TNF, TEK, BMPR2, SERPINE1, IL10, CD40LG, CCL2, FASLG and ABL1 can be distinguished. These genes can be used to plan experiments confirming their association with lymphedema, and to be further considered as new candidates for genotyping. The found genes can become the basis for a better understanding of the mechanisms of interaction between endothelial apoptosis and lymphedema.
The over-represented Gene Ontology (GO) biological processes were identified using the DAVID 6.8 tool [31] with the following parameters: the organism -«Homo sapiens», Gene_Ontology -«GOTERM_BP_DIRECT».
Reconstruction of associative gene networks was carried out by using the ANDSystem [21,22].
The betweenness centrality of a node in a gene network was estimated using the networkx package of the Python programming language [32]. This indicator characterizes the number of shortest pathways between all pairs of vertices of analyzed graph passing through a given vertex and reflects the functional significance of gene in gene network.
The Mann-Whitney criterion was calculated using the mannwhitneyu function of the scipy.stats package of the Python programming language [33].
The list of human genes involved in the Gene Ontology biological process «apoptotic process» was obtained using the AmiGO system [34] by the «GO:0006915» query and Organism filter set to «Homo sapiens».
The list of human genes involved in the Gene Ontology biological process «endothelial cell apoptotic process» was obtained with AmiGO using the «GO:0072577» query and with Organism filter set to «Homo sapiens».
For the gene prioritization six criteria were used ( Fig. 1), which were discussed in our previous studies devoted to the analysis of asthma, hypertension and Parkinson's disease [19,20].
Criteria 1 and 2 were the well-known prioritization methods: ToppGene [15] and Endeavour version 3.71 [18], respectively. These resources allow to perform ranking of a test set of genes by a training set of genes according to specific criteria characterizing the proximity of genes from the test set to the genes from the training sample. The methods of these resources use the genetic information (co-localization in the genome), the functional properties of genes (participation in the same GO categories), etc., as well as properties of the vertices of the graph of protein-protein networks. A list of genes associated with lymphedema, described above, was used as a training set for each of these methods. As a test set, the list of genes involved in endothelial cell apoptosis, described above, was used. For the ToppGene the "all Feature" parameter was selected in the "Training parameters" section, and the ranking of genes was based on the value of the "Rank" output parameter. In case of the Endeavour system all settings were set to default, the gene ranking was based on "P-value". Thus, the lowest ranks had genes with the lowest "P-value", while genes with the highest "P-value" obtained the highest ranks.
Criterion 3 was calculated as the proportion of the biological processes, over-represented for a set of genes associated with lymphedema (Additional file 2: Table S2), to all Gene Ontology biological processes where the analyzed gene was involved. Information on the involvement of gene in the Gene Ontology biological process was obtained from the AmiGO system [34]. Ranks for genes were determined with the sorting of the list of genes by descending the proportion of over-represented Gene Ontology biological processes. Thus, the lowest ranks were assigned to genes with the largest proportion of over-represented Gene Ontology biological processes.
Criterion 4 was based on the use of the cross-talk centrality (CTC), calculated using the «Intelligent Filtration» function of ANDSystem. Within this criterion, CTCs were calculated separately for genes from the gene network (criterion 4A) and for their products (criterion 4B). Thus, two indices (CTC gene and CTC protein ) were determined for each gene. The centrality of CTC gene was calculated using the following formula: where N iis a number of interactions of the i-th gene with members of the associative gene network of lymphedema; Mis a total number of nodes from the associative gene network of lymphedema. CTC protein was determined for the i-th protein or miRNA in a same way. CTC gene and CTC protein ranks were determined by sorting the values in descending order the same way as it was done for the criterion 3.
Criterion 5 represents the number of diseases comorbid to lymphedema, which are associated with the analyzed gene. Comorbid diseases are considered to be simultaneously present in one patient more often than can be expected for accidental reasons [35,36]. The list of diseases comorbid to lymphedema was manually created by analyzing the publications corresponding to the following query to the PubMed database: "lymphedema and (comorbid or comorbidity)". A total of 80 publications were manually analyzed and six comorbid diseases were found (Additional file 3: Table S3). All interactions between the analyzed genes and these six comorbid diseases were established using ANDSystem. In addition, the rank by criterion 5 was determined by sorting the list of genes in descending order of the number of comorbid diseases associated with these genes.
In the cases when several genes had the same value of the criterion by which the ranking was performed, the rank of these genes was calculated similarly to the Spearman rank correlation coefficient [37]. I.e., such genes were assigned with the same rank, equal to the average arithmetic ranks of these genes, according to their position, in the sorted list of genes.
Criterion 6 was calculated as the average value of ranks obtained from the criteria 1-5.
The statistical significance of the co-occurrence of gene names from the test set and lymphedema in full-text articles was estimated as follows: (1) The number of full-text publications presented in the PubMed Central system (https://www.ncbi.nlm.nih.gov/pmc/) was calculated according to the following parameters: (a) the official name of the analyzed gene is co-occurring with the "lymphedema" word, (b) the official name of the analyzed gene is mentioned in the article, Fig. 1 The prioritization pipeline overview (c) the word "lymphedema" is mentioned in the article.
(2) Using the hypergeometric distribution, implemented in the hypergeom.sf function of the scipy.stats package of the Python programming language [33], the statistical significance of the co-occurrence of names of the analyzed genes and lymphedema was assessed.
The correlation between the ranks of criteria 1-6 and the significance of co-occurrence of the analyzed genes names and lymphedema in the full-text articles was carried out using the Spearman's rank correlation coefficients. The calculations were performed using the Social Science Statistics system (http://www.socscistatistics.com/tests/spearman/default2.aspx).

Results and discussion
Associative gene network of lymphedema Analysis of information from the CTD [25], Malacards [26], KEGG [27], HPO [28] and DisGeNET [29] databases and from ANDSystem, revealed 69 genes associated with lymphedema (Additional file 1: Table S1). The obtained list of genes was applied to ANDSystem for the reconstruction of the associative gene network of lymphedema. The reconstructed network was automatically expanded in ANDSystem with the products of these genes. The obtained gene network contained 69 genes and 78 proteins, as well as 709 interactions between them (Fig. 2), including the following interaction types: association (460 interactions), protein-protein interaction (128 interactions), gene expression (78 interactions), regulation of gene expression (30 interactions), co-expression (6 interactions), regulation of transport (2 interactions), regulation of degradation (2 interactions), catalysis (2 interactions) and regulation of activity (1 interaction). Twenty one genes from the network were interacting only with their own products, while the remaining 48 genes and 57 proteins were highly connected with each other (Fig. 2). Analysis of the overrepresentation of Gene Ontology (GO) biological processes for a set of genes associated with lymphedema revealed 40 processes significantly enriched with these genes (Additional file 2: Table S2). Among the most over-represented processes were lymphangiogenesis, endothelial cell proliferation, ERK1/ ERK2 cascade and VEGF signaling pathway. The role of these biological processes in the pathogenesis of lymphedema is actively discussed in the literature [8,[40][41][42]. For example, Saito et al., 2013 studied the possibility of therapeutic lymphangiogenesis, using hepatocyte growth factor for the treatment of lymphedema [40]. Miaskowski et al., 2013, described the important role of VEGF signaling pathway in lymphangiogenesis associated with inflammation in lymphedema [8].
It should be noted that among the significantly over-represented GO biological processes were also processes of positive and negative regulation of the apoptosis (Additional file 2: Table S2); this is consistent with recent studies suggesting the important role of apoptosis in the pathogenesis of lymphedema [14]. The associative gene network of lymphedema involved 24 genes/proteins from the GO biological process entitled "apoptotic process" (Fig. 2). Analysis of the centrality of nodes from the gene network of lymphedema (Additional file 4: Table S4) showed that the average value of the betweenness centrality for apoptosis genes (equal to 303) statistically significantly exceeds the value calculated for all the nodes of the gene network (equal to 128), with p-value< 10 − 6 , according to the Mann-Whitney criterion. The high centrality of the apoptotic genes in the lymphedema gene network can be indirect evidence of the key role of the apoptosis process in the pathogenesis of lymphedema.

Prioritization of endothelial apoptosis genes by their potential association with lymphedema
Considering the importance of apoptosis of endothelial cells in the pathogenesis of lymphedema [10,11,13,14], we performed a prioritization of genes involved in the «endothelial cell apoptotic process» (GO: 0072577) using a training set that included genes associated with lymphedema. The prioritization process was based on the use of well-known methods (ToppGene and Endeavour), as well as on our original approaches, allowing to consider the structure of the gene network of lymphedema, over-represented biological processes and associations of genes with diseases comorbid to lymphedema. It appeared that two genes GATA2 and KDR were simultaneously presented in the list of genes from the «endothelial cell apoptotic process», as well as in the list of genes associated with lymphedema [43,44]. During the prioritization, these two genes were excluded from the training set and were used as control genes. Thus, the prioritization was performed for 64 genes from the testing set and 67 genes from the training set, associated with lymphedema.
According to the criterion 6 (Additional file 5: Table  S5), calculated as an average value of criteria 1-5, it turned out that the KDR control gene was on the first place, while the GATA2 gene was on twelfth. Hypergeometric distribution showed that the enrichment of the top 12 genes by these two control genes is statistically significant (p-value = 0.03). It should be noted that a similar analysis, carried out by individual criteria 1-5, showed the absence of any statistically significant enrichment. Thus, the criterion 6, which takes into account five previous criteria, can be considered as the best approach to the genes prioritization.
In the Table 1 are shown top 10 of the highest priority genes according to the criterion 6. Among these top 10 genes SNPs associated with lymphedema are known only for KDR gene [44], which had the first place in the table.
For the other genes, to our knowledge, there was no information about any SNPs associated with this disease and these genes can be considered as new candidates for lymphedema susceptibility.
Also, it was interesting to consider each remaining criterion separately, even though none of them showed any statistical significance during the evaluation of the enrichment of the top of genes, ranked by these criteria, by the control genes. Thus, both methods of prioritization (ToppGene and Endeavour) determined the KRD control gene in the first place in the resulting list. The top 5 of the most priority genes according to both ToppGene (Additional file 6: Table S6) and Endeavour (Additional file 7: Table S7) were KDR, ABL1 and TEK, presented in Table 1. Using the criterion 3 the top 10 included the control gene KDR and the TEK gene (Additional file 8: Table S8). At the same time, according to criterion 4A (Additional file 9: Table S9), the top 10 included eight genes (KDR, TNF, IL10, CCL2, CD40LG, ABL1, BMPR2 and SERPINE1) from the Table 1. While for the criterion 4B (Additional file 10: Table S10) there were seven such genes in a list (TNF, IL-10, KDR, CD40LG, CCL2, ABL1 and BMPR2). According to the criterion 5, which considers the comorbidity of diseases, the TNF gene, associated with six comorbid diseases, was in the first place (Additional file 11: Table S11), the second priority was given to CCL2, IL10 and SERPINE1 genes, associated with five comorbid diseases. These genes also appeared to be in the top 10 list ( Table 1).
Verification of the prioritization results using the full-text articles Jenssen et al., 2001 proposed an approach for automated identification of potential interactions between biological objects, based on an assessment of the statistical significance of the co-occurrence of terms in scientific publications [45]. Because this approach is not implemented in ANDSystem the estimation of the co-occurrence between the analyzed genes and lymphedema can serve as an indirect evidence of the correctness of our prioritization results. It can be expected that the most priority genes would have the highest frequency of mentioning in articles together with lymphedema. As the ANDSystem knowledge base was created by automated analysis of PubMed, to obtain results that could be even more free of the ANDSystem data, we performed a manual analysis of the full-text articles from PubMed Central (Additional file 12: Table S12). It was found that of the 64 of analyzed genes, the 21 genes were statistically significantly more often mentioned in publications together with the "lymphedema" word (p-value< 0,05 with FDR correction). 13 out of these 21 genes appeared to be presented among the top 20 genes, ranked according to criterion 6 ( Table 2), while the most significant association with lymphedema was observed for the GATA2 and KDR control genes (p-value = 10 − 195 and p-value = 10 − 136 , respectively).
According to the hypergeometric distribution, the statistical significance of such enrichment of the top 20 genes had a p-value< 0.0004. It should be noted that only criterion 5, which considers the comorbidity, also had such a low p-value, while the p-values of other criteria were higher. At the same time, the highest Spearman's rank correlation coefficient between the genes ranks and the significance of co-occurrence of the analyzed genes names and lymphedema in the full-text articles (r = 0,529), was observed for the criterion 6 ( Table 2).

Agreement between the prioritization results and gene expression data from the mouse lymphedema model
Lymphedema is characterized by a chronic stasis of lymph in the tissues. In Tabibiazar Table S13) were detected by GEO2R (http://www.ncbi.nlm.nih.gov/geo/ geo2r/). After the mapping of these differentially expressed probes on orthologues human genes from the 10,201 probes 3289 differentially expressed human genes were remained, while of the 2219 probes remained 735 genes. The remaining probes were removed from the analysis due to the lack of any intersections with human genes. After combining these two sets 3494 differentially expressed genes were obtained. It appeared that 13 of 64 analyzed genes were in the combined set. According to criterion 6, among the top 20 of the highest priority genes there were 8 of 13 differentially expressed genes (KDR, PLCG1, SERPINE1, CD40LG, IL10, CCL2, PDPK1 and THBS1). It was shown by hypergeometric distribution that such enrichment is statistically significant (p-value = 0.012, Table 3). However, the greatest enrichment value was observed for criterion 4B. This criterion is based on an assessment of the CTC centrality of the nodes of the gene network, corresponding to proteins.

Interactions of the top 10 candidates with participants of the lymphedema gene network
The associative gene networks reconstructed by ANDSystem were used for the analysis of interactions of the top 10 candidate genes, obtained by criterion 6, with the genes/proteins associated with lymphedema (Additional file 14: Table S14, Fig. 3).
The highest priority, according to Table 1, had the control gene KDR, which encodes Vascular Endothelial Growth Factor Receptor 2. KDR is one of the main regulators of endothelial cell growth [46], it provides endothelial survival [47], proliferation, migration [48], sprouting [49] and tubular morphogenesis [50,51]. It is known that a number of polymorphisms in this gene are associated with lymphedema with extreme risk estimates OR > 2, p-value< 0.05 [44]. In the mouse tail lymphedema model, using the quantitative real-time polymerase chain reaction, it was shown that the messenger RNA (mRNA) expression of KDR was increased in the group treated by the low-level laser therapy compared to the control [52]. The associative gene network describing the interactions of the KDR gene consisted of 14 nodes (Fig. 3a). In particular, the KDR protein had protein-protein interactions with VEGF-C, VEGF-D, VEGFR-3 and CBL [53][54][55][56]. It is known that mutations in VEGFC, VEGFR-3 and CBL genes are associated with lymphedema [4,[57][58][59][60][61][62][63]. Serum level of VEGF-D was shown to be significantly higher in the group of patients with primary lymphedema compared with controls [64]. It is also known from the literature that the KDR and SHP2 genes have a positive correlation of their expression levels [65]. Mutations in the SHP2 gene are associated with lymphedema [66]. HGF in combination with VEGF-A can activate VEGFR-2 [67]. Mutations in the HGF gene are also associated with lymphedema [68,69], and increased expression of HGF improves lymphedema [70,71].
The second place in Table 1 belonged to the Tumor Necrosis Factor alpha (TNF) gene, which encodes cytokine involved in systemic inflammation [72]. Figure 3B shows the 20 nodes associated with this gene, in particular, TNF was associated with IL-6. It is known from the literature that this protein is capable to activate the expression of the IL-6 gene [73], increase stability and secretion of IL-6 [74,75]. In turn, IL-6 has been shown to be increased in models of lymphedema [76]. Another example is the HMGB1 protein, which was also connected to TNF in the associative gene network. It is known that HMGB1 can induce expression of TNF [77]. Also studies show the 2.4-fold increased level of HMGB1 in secondary lymphedema [78]. Besides, TNF was found to be associated with adiponectin in the network, a level of which can be reduced in response to the TNF [79]. Shimizu et al., 2013 showed that adiponectin can promote a lymphatic vessel formation resulting in amelioration of lymphedema [80]. The potential role of TNF in the molecular mechanisms of lymphedema is discussed in the literature. It is known that TNF can induce endothelial cells apoptosis [81][82][83]. Földi et al., 2000 showed a decrease in the level of expression of TNF after the complex decongestive physiotherapy in patients with peripheral leg lymphedema [84]. Ji, 2007, reports the importance of the role of TNF in the pathophysiological changes of the lymphatic endothelial cells (LECs) and inflammatory lymphangiogenesis [85]. Anuradha et al., 2012, shows the association of inflammatory cytokines (IL-1β, IL-12, and TNF) with pathogenesis of lymphatic filarial infection, which can lead to lymphedema [86]. Jeong et al., 2013 by using the western blot analysis, showed that hyaluronidase treatment of lymphedema in a mouse tail model declined expression of TNF [87].
The third place in Table 1 had the TEK genea receptor tyrosine kinase, involved in a signaling pathway related to the embryonic vascular development [88]. It is known that TEK can reduce proliferation of Fig. 3 Gene networks describing the interactions of the top 10 of the most promising candidate genes a KDR, b TNF, c TEK, d BMPR2, e SERPINE1, f IL10, g CD40LG, h CCL2, i FASLG and j ABL1 with genes/proteins associated with lymphedema, reconstructed by ANDSystem. Candidate genes and their products are illustrated with large icons endothelial cell and apoptosis [89], while mutations in this gene can cause autosomal dominantly inherited forms of venous malformations [90,91]. The reconstructed gene network contained four interactions of the TEK gene with genes/proteins (FOXC2, NEMO, Shp2 and CBL) associated with lymphedema (Fig. 3c). From the literature it is known that the loss of FOXC2 can result in TEK downregulation [92] meanwhile, a complete loss or a significant gain of FOXC2 function can lead to perturbation of lymphatic vessel formation and lymphedema [93]. According to information from the Innatedb [94], HPRD [95] and BioGrid [96] databases, the TEK protein can interact with the NEMO, Shp2 and CBL proteins, mutations in which are associated with lymphedema [63,66,97].
The bone morphogenetic protein receptor type II (BMPR2) genea serine/threonine receptor kinase involved in cell growth and differentiation, osteogenesis and adipogenesis [98] was on the fourth place. A reduced BMPR2 expression can induce mitochondrial dysfunction of endothelial cells, promoting a pro-inflammatory and pro-apoptotic state [99,100]. Kim, Kim, 2014 showed that knockdown of Bmpr2a and Bmpr2b result in lymphatic defects in developing zebrafish [101]. In the constructed associative gene network, BMPR2 has four interactions with the VEGFC, IL6, HGF and FOXL1 genes/proteins (Fig. 3d). It is known that silencing of BMPR2 results in down-regulation of VEGFC expression [102]. Also, VEGFC plays an important role in the functioning of lymphatic vessels while mutations in this gene are associated with lymphedema [57][58][59][60]. Soon et al., 2015 showed that mutations in BMPR2 lead to higher levels of IL6 [103], which level is increased in lymphedema models [76]. The expression of BMPR2 was up-regulated by HGF [104]. The increase of expression of HGF improves lymphedema [71,105], and mutations in this gene are associated with the disease [68,69]. According to the HPRD database, the FOXL1 protein interacts with BMPR2, and mutations in the FOXL1 gene are associated with lymphedema [106].
Serpin Family E Member 1 (SERPINE1), which is the inhibitor of tissue plasminogen activator and urokinase [107] that regulates fibrinolysis [108], appeared to be in the fifth place of the table. The SERPINE1 gene mediates the anti-apoptotic effect in the endothelial cell [109]. A number of studies have shown the association of SER-PINE1 with the metastasis of tumors in the lymph nodes [110][111][112][113]. In the associative gene network, SERPINE1 interacts with four genes/proteins: IL6, HGF, HMGB1 and APN (Fig. 3e). It is known that SERPINE1 expression decreases with the addition of IL-6 [114] that is increased in lymphedema models [76]. HGF increases SERPINE1 expression [115], and increase of HGF expression improves lymphedema [70,71]. It was shown that levels of HMGB1 and SERPINE1 had positive correlation [116], while correlation of APN and SERPINE1 was negative [117]. The administration of APN improved the edema of injured tails in the mouse model of lymphedema [80].
The sixth in the Table 1 was an anti-inflammatory cytokine IL10, playing an important role in the immunoregulation and inflammation [118]. It was shown that IL-10 significantly blocked endothelial apoptosis [119]. However, other study showed that IL-10 had the capacity to induce macrophage apoptosis [120]. An increased gene expression of IL10 was found in keratinocytes derived from limb affected by lymphedema [121] and in wounded lymphedematous skin [122]. In the associative network IL-10 appeared to be connected with the APN and HGF genes (Fig. 3f ). It is known that APN can promote amelioration of lymphedema [80] and is able to significantly increase IL-10 gene expression and protein secretion [123]. HGF can increase plasma IL-10 concentration [124,125].
The seventh line in the Table 1 was taken by the CD40LG gene, which is expressed in T cells and its product is presented on their surface [126]. CD40LG can bind CD40 on the B cell surface and is involved in T cell proliferation and cytokine production [127]. It is known from the literature that the product of the CD40LG gene can increased apoptosis of endothelial cells [128]. This gene was found to be interacting only with IL-6 in the associative gene network (Fig. 3g). Sommer et al., 2009 showed that CD40LG can up-regulate expression of IL-6 [129].
On the eighth line was the CCL2 (C-C Motif Chemokine Ligand 2) gene, which encodes the cytokine possessing chemotactic activity for monocytes [130] and basophils [131]. Down-regulation of CCL2 by miR-495 resulted in inhibited apoptosis of human umbilical vein endothelial cells [132]. In the associative gene network this gene was interacting with 5 genes/proteins (Fig. 2h). Among the genes involved in interactions with CCL2 were APN, CDC42, HGF, IL6 and TSC2, which are associated with lymphedema [68-71, 76, 80, 133-135]. It is known that APN elevates mRNA and protein level of the CCL2 and stimulates release of CCL2 in primary human monocytes [136]. Ablation of CDC42 induced an overexpression of CCL2 [137]. Müller et al., 2012 showed that CCL2 can induce HGF [138]. It was shown that IL6 induced CCL2 [139] and the IL6 level positively correlated with the CCL2 [140]. Loss of TSC2 function can result in overexpression of CCL2 [141].
The FASLG gene, located on the ninth row of the Table 1, is a member of tumor necrosis factor superfamily [142,143]. FASLG is encoding transmembrane protein Fas ligand that is responsible for the induction of apoptosis triggered by binding to FAS [144][145][146]. Abnormal FASLG expression is associated with a metastasis of tumors in the lymph nodes [147][148][149]. In the reconstructed associative gene network FASLG appeared to be associated only with IL6 (Fig. 2i). It is known that FASLG stimulation can enhanced IL-6 release [150].

Conclusion
In this work, we performed a search for new potential participants of lymphedema molecular-genetic mechanisms based on the prioritization of genes involved in endothelial cell apoptosis. Six criteria, including the well-known ToppGene and Endeavour methods, as well as our original approaches [19,20] were used. The use of original methods allowed taking into account the overrepresented Gene Ontology biological processes, structural features of the associative gene network of lymphedema and endothelial apoptosis, as well as diseases comorbid to lymphedema. Verification of the prioritization quality using three different criteria showed significant enrichment of the most priority genes with known genes, which have simultaneous interactions with lymphedema and endothelial apoptosis, as well as with genes differentially expressed in the murine model of lymphedema. Among genes, involved in endothelial cell apoptosis, TNF, TEK, BMPR2, SERPINE1, IL10, CD40LG, CCL2, FASLG and ABL1 were identified as the most promising candidates that can be used for planning the experiments concerning their possible associations with lymphedema. Besides, analysis of the function of these genes can help in understanding the moleculargenetic role of endothelial cell apoptosis in lymphedema.

Additional files
Additional file 1:

Funding
The publication cost was covered by the Russian Science Foundation grant «Programmed cell death induced via death receptors: Delineating molecular mechanisms of apoptosis initiation via molecular modeling» 14-44-00011.

Availability of data and materials
Results are shared in the additional files.

About this supplement
This article has been published as part of BMC Medical Genomics Volume 12 Supplement 2, 2019: Selected articles from BGRS\SB-2018: medical genomics. The full contents of the supplement are available online at https:// bmcmedgenomics.biomedcentral.com/articles/supplements/volume-12supplement-2.
Author's contributions VAI and OVS development of criteria for prioritization and gene selection. OVS, DBU, VVN, PSD and TVI performed the reconstruction of gene networks. OVS, DBU, VVN and VAI produced results on analysis of lymphedema gene network. OVS, VAI and INL performed analysis of apoptosis gene network. OVS and VAI produced results on analysis of full-text articles and differentially expressed genes. OVS, VVN and VAI drafted the manuscript. VAI, OVS, VVN, INL and TVI revised the manuscript. VAI supervised the whole study. All authors have read the manuscript and approved the final version.