In silico analysis of the molecular regulatory networks in peripheral arterial occlusive disease.

BACKGROUND: Peripheral arterial occlusive disease (PAOD) is a global public health concern that decreases the quality of life of the patients and can lead to disabilities and death. The aim of this study was to identify the genes and pathways associated with PAOD pathogenesis, and the potential therapeutic targets.
METHODS: Differentially expressed genes (DEGs) and miRNAs related to PAOD were extracted from the GSE57691 dataset and through text mining. Additionally, bioinformatics analysis was applied to explore gene ontology, pathways and protein-protein interaction of those DEGs. The potential miRNAs targeting the DEGs and the transcription factors (TFs) regulating miRNAs were predicted by multiple different databases.
RESULTS: A total of 59 DEGs were identified, which were significantly enriched in the inflammatory response, immune response, chemokine-mediated signaling pathway and JAK-STAT signaling pathway. Thirteen genes including IL6, CXCL12, IL1B, and STAT3 were hub genes in protein-protein interaction network. In addition, 513 miRNA-target gene pairs were identified, of which CXCL12 and PTPN11 were the potential targets of miRNA-143, and IL1B of miRNA-21. STAT3 was differentially expressed and regulated 27 potential target miRNAs including miRNA-143 and miRNA-21 in TF-miRNA regulatory network.
CONCLUSION: In summary, inflammation, immune response and STAT3-mediated miRNA-target genes axis play an important role in PAOD development and progression.

Introduction

Peripheral artery occlusion disease (PAOD) is an atherosclerotic condition involving non-cardiac and non-cerebral arteries. Approximately 202 million cases of PAOD were diagnosed worldwide in 2010,[ and is estimated to affect almost 41.13 million people in China by 2020 due to a rapidly ageing population.[ PAOD results in severe arterial stenosis and insufficient blood supply to the distal limb causing intermittent cicada, resting pain, disability and even death. PAOD-related mortality rate increases from 0.07/105 in 40 to 44-year age group to 28.71/105 among those older than 80 years. The highest proportion of patients with disability-adjusted life years in 2010 was 47.88/105 in Western Europe.[ PAOD is closely associated with coronary events and cerebrovascular disease,[ and other risk factors include smoking, diabetes, hypertension, hypercholesterolaemia, and aging. Although anti-inflammatory, anticoagulant and antiplatelet drugs, surgical intervention and endovascular treatments[ can improve blood supply and reduce PAOD complications, the atherosclerotic lesions and arterial stenosis are unaffected. Gene and cell therapies have been developed in recent years, but the therapeutic outcomes have been less than satisfactory.[ Therefore, it is essential to dissect the molecular basis of PAOD pathogenesis in order to identify novel therapeutic targets.

MicroRNAs (miRNAs) are a class of endogenous single-chain non-coding RNAs consisting of 18 to 22 nucleotides that regulate gene expression at the post-transcriptional level by binding to a seed sequence.[ miRNAs are associated with a wide range of physiological and pathological processes, and have been identified as therapeutic targets for various conditions, including cardiovascular diseases.[ For example, circulating levels of miR-21, miR-218, and miR-211 are potential diagnostic fingerprints of diabetic atherosclerosis.[ In addition, miR-172 was identified as a target of anti-atherogenic drugs in the ApoE-/- atherosclerotic mouse model.[ Recently, miRNA-21 and miRNA-143 were identified as biomarkers of peripheral arterial disease and vascular restenosis after endovascular treatment,[ indicating their involvement in the development and progression of PAOD.

Transcription factors (TF) regulate gene expression by binding to the promoter of target genes, and are often dysregulated during pathological conditions. Studies show that some recombinant TFs can promote angiogenic growth factor expression and improve perfusion in limb ischemia and angiogenesis.[ STAT3 is closely related to atherosclerosis, and its inhibition reduced atheromatous plaque formation both in rabbit and murine atherosclerosis models.[ MiRNA-21, the confirmed biomarker of PAOD, is also a target of STAT3.[ Taken together, a hitherto unknown TF/miRNA axis likely plays an important role in driving PAOD. In this study, we mined gene expression data of PAOD and control specimens from published studies and transcriptomic datasets, and identified crucial signaling pathways, genes and miRNAs involved in PAOD. We also constructed miRNA/target gene and TF/miRNA/target gene networks relevant to PAOD. Our findings provide new insights into the pathogenesis of PAOD and identify potential therapeutic targets.

Methods

Data collection

The GSE57691 dataset based on GPL10558 Illumina HumanHT-12 V4.0 expression beadchip, including gene expression data of aortic specimens from 9 AOD and 10 control donors, was downloaded from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Differentially expressed genes (DEGs) between the AOD and control aortic specimens were identified using the GEO2R online tool, with |LogFC| > 1 and P < .05 as the criteria, and presented in the form of Volcano plots using the ggplot2 package in R software. Genes related to PAOD were mined from published genomic data and MEDLINE literature using pubmed2ensembl (http://pubmed2ensembl.ls.manchester.ac.uk/),[ with “Peripheral arterial occlusive disease,” “Homo sapiens genes (GRCh37) database,” and “filter on MEDLINE” as the queries. The genes common to the DEGs in GSE57691 and identified from text mining were analyzed further, and henceforth referred to as PAOD-related genes. The miRNAs associated with PAOD were similarly mined from the PubMed database using search terms “Arterial occlusive disease” [Title/Abstract], “Peripheral arterial occlusive disease” [Title/Abstract], “Lower extremity arterial occlusive disease” [Title/Abstract], “peripheral arterial disease” [Title/Abstract] AND “miRNA” OR “microRNA”. Additional approval by an ethics committee was not necessary because the datasets included in the current study were downloaded from public databases, and data acquisition and application were performed according to GEO publication guidelines and data access policies.

Gene ontology (GO) and pathway enrichment analysis

The database for annotation, visualization and integrated discovery (DAVID) version 6.8 (https://david.ncifcrf.gov/home.jsp) was used for the GO functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the PAOD-related genes identified above. P value < .05 was set as the threshold for the enriched GO terms and significant pathways.

Protein–protein interaction (PPI) network and modular analysis

The PPI network of the relevant genes was constructed using STRING (https://string-db.org/) version 11,[ and analyzed by the Cytoscape software. The nodes with more than 10 degrees in the PPI were considered hub genes. Modular analysis of PPI was also performed with MCODE app in Cytoscape software, and the subsequent enrichment analysis in DAVID database.

Prediction of cognate miRNAs of PAOD-related genes and construction of the miRNA-target gene regulatory network

The miRNAs targeting the PAOD-related were predicted using the miRanda (http://www.microrna.org/), miRDB (http://www.mirdb.org/), miRWalk (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk/), RNA22 (https://cm.jefferson.edu/rna22/Interactive/), and TargetScan (http://www.targetscan.org/vert_72/) databases, and those predicted by at least 4 databases were selected for constructing the miRNA-target gene regulatory network by Cytoscape software.

Prediction of transcription factors of miRNA and construction of the TF-miRNA-target gene regulatory network

The TFs regulating miRNAs in the miRNA-target gene regulatory network were predicted using literature-curated TF-miRNA regulation data in the TransmiR (http://www.cuilab.cn/transmir) database.[ The genes common to the predicted TFs and PAOD-related genes were considered differentially expressed TFs (DETFs). The TF-miRNA-target gene, DETF- miRNA-target and PPI hub gene-miRNA-target gene regulatory networks were constructed using Cytoscape software. The TF binding sites on miRNAs and the miRNA binding sites on target genes were predicted by the JASPAR (http://jaspar2018.genereg.net) and microRNA.org (http://www.microrna.org/microrna/) databases respectively.

Results

Identification of PAOD-related genes

A total of 1705 DEGs were identified in the GSE57691 dataset, including 75 up-regulated and 1630 down-regulated genes (Fig. 1 and Supplementary data-1). In addition, 591 genes related to PAOD were identified after text mining (Supplementary data-2). The overlap of both gene sets revealed 59 genes, including 15 up-regulated and 44 down-regulated genes (Fig. 2 and Table 1), that were henceforth analyzed as PAOD-related DEGs. Sixteen miRNAs, including miRNA-320a, miRNA-572, miRNA-21 and miRNA-143, were among the PADO-related DEGs (Supplementary data-3).

Figure 1

Volcano plot showing DEGs between PAOD and control samples in the GSE57691 dataset. Green nodes indicate down-regulated genes with logFC < -1 and P < .05. Red nodes indicate up-regulated genes with logFC > 1 and P < .05. FC = fold change.

Figure 2

Venn diagram showing DEGs in GSE57691 and from text mining. (A) GSE57691 and text mining respectively identified 60 and 576 up-regulated genes, with 15 overlapping genes, and (B) 1586 and 547 down-regulated genes, with 44 common genes. DEGs = differentially expressed genes.

Table 1

A total of 59 shared DGEs between GSE57691 dataset and gene sets obtained from text mining were found, including 15 up-regulated and 44 down-regulated genes.

GO annotation and KEGG pathway enrichment analyses

To determine the biological relevance of the overlapping PAOD-related DEGs, we next determined the biological process (BP), cellular component (CC) and molecular function (MF) annotations of GO analysis. The significantly enriched BP terms included inflammatory response (GO:0006954, GO:0002675), immune response (GO:0006955), cell adhesion (GO:0007155), regulation of cell proliferation (GO:0008284, GO:0008285), regulation of apoptotic process (GO:0043066), regulation of ERK1 and ERK2 cascade (GO:0070374), regulation of MAPK cascade (GO:0000187, GO:0043410), and cytokine and chemokine-mediated signaling pathway (GO:0019221, GO:0070098) (Fig. 3A and Supplementary data-4). Protein binding (GO:0005515), growth factor activity (GO:0008083), chemokine activity (GO:0008009), cytokine activity (GO:0005125) and cell adhesion molecule binding (GO:0050839) were the significantly enriched MF terms among the DEGs (Fig. 3B and Supplementary data-4). Finally, extracellular space (GO:0005615), extracellular region (GO:0005576), extracellular exosome (GO:0070062), plasma membrane (GO:0005886, GO:0005887), cell surface (GO:0009986) and focal adhesion (GO:0005925) were the significantly enriched CC terms (Fig. 3C and Supplementary data-4). The KEGG pathway analysis showed that DEGs were mainly associated with hsa04510: Focal adhesion, hsa04060: Cytokine-cytokine receptor interaction, hsa04062: Chemokine signaling pathway, hsa04670: Leukocyte transendothelial migration, hsa04151: PI3K-Akt signaling pathway, and hsa04630: Jak-STAT signaling pathway (Fig. 3D and Supplementary data-4).

Figure 3

GO and KEGG pathway enrichment of PAOD-related DEGs. The significantly enriched (A) BP, (B) MF and (C) CC terms, and (D) pathways. BP = biological process, CC = cellular component, DEGs = differentially expressed genes, MF = molecular function.

PPI network and modular analysis

To determine the potential interactions between the DEGs and identify the hub genes, a PPI network was established using 47 genes (12 up- and 35 down-regulated genes) and 180 edges (Fig. 4). IL6, CXCL12, IL1B, ITGB1, APP, CAV1, STAT3, CTGF, SPARC, CDC42, SOD1, PPARG, and TIMP2 had more than 10 degrees, and were defined as hub genes. In addition, module1 and module2 were clustered from the PPI network. The former had 8 genes and 21 edges that were enriched in regulation of cell proliferation (GO:0008285, GO:0042102), regulation of MAPK cascade (GO:0043410), extracellular region (GO:0005576), positive regulation of JNK cascade (GO:0046330), and hsa04621: NOD-like receptor signaling pathways (Fig. 5A and Supplementary data-5). The latter consisted of 14 genes and 31 edges, and were enriched in chemotaxis (GO:0006935, GO:0070098, GO:0060326, GO:0030593, GO:0008009), immune response (GO:0006955), cell adhesion (GO:0007155, GO:0005925, GO:0050839), inflammatory response (GO:0002523), and hsa04062: Chemokine signaling pathway and hsa05205: Proteoglycans in cancer pathways (Fig. 5B and Supplementary data-5).

Figure 4

PPI network of PAOD-related DEGs showing 47 nodes and 180 edges. Red and blue nodes represent the up- and down-regulated genes respectively. Greater node size corresponds to higher degree. Greater edge size and darker color correspond to higher combined score. DEGs = differentially expressed genes, PPI = protein–protein interaction.

Figure 5

Module analysis of PPI network. (A) module1 and (B) module2.PPI = protein–protein network.

The miRNA-target gene regulatory network

To identify the potential miRNAs involved in PAOD, those targeting the PAOD-related DEGs were predicted using the miRanda, miRDB, miRWalk, RNA22, and TargetScan databases, which revealed 311 miRNAs (Supplementary-data-6) including hsa-miR-143 and hsa-miR-21. A total of 513 miRNA-target gene pairs were identified, including miRNA-143-PTPN11, miR-143-CXCL12 and miRNA-21-IL1B. A PPI network was constructed with 311 miRNAs, 40 target genes and 513 edges (Fig. 6).

Figure 6

miRNA-target gene regulatory network. Green nodes represent miRNAs. Red and blue nodes represent up- and down-regulated target genes respectively. Greater node size corresponds to higher degree. Arrows indicate miRNA-target gene relationship.

TF-miRNA-target gene regulatory network analysis

The TFs regulating the PAOD-related miRNAs were predicted using TransMir, which revealed 295 TFs, of which STAT3, PPARG, USF1, PDGFA, IL1B, and IL6 were DETFs. The TF-miRNA regulatory network was constructed, and consisted of 434 nodes and 1261 edges including 295 TFs and 139 miRNAs (Fig. 7A and Supplementary data-7). The hub genes identified in the PPI network were also selected to construct the TF-miRNA-target gene regulatory network (Fig. 7B), which showed 251 nodes, 191 TFs, 53 miRNAs, 10 hub genes, and 699 edges (Supplementary data-8). Finally, the DETF-miRNA-target gene network was also constructed, including 39 miRNAs, 22 target genes and 6 DETFs, of which the DETF STAT3 had the highest degree (Fig. 7C and Supplementary data-9). Hsa-miR-143 and hsa-miR-21 were predicted as the targets of STAT3 (Fig. 8A and B), PTPN11 as a target DEG of hsa-miR-143, and STAT3 and IL1-B as the target DEGs of hsa-miR-21 (Fig. 8C).

Figure 7

Molecular regulatory networks in PAOD. (A) TF-target miRNA network. (B) TF-target miRNA-target hub gene network. Blue nodes represent hub genes in PPI network. (C) DETF-target miRNA-target gene network. Red, green and blue nodes represent miRNAs, TFs and target genes respectively. Greater node size corresponds to higher degree. Arrows indicate TF-target miRNA relationship. “T” shape indicates miRNA-target gene relationship. DETF = differential expression transcription factor, PPI = protein-protein network, TF = transcription factor.

Figure 8

STAT3-target miRNA-target gene axis. (A) Binding site consensus sequence of STAT3. (B) Binding site between STAT3 and miR-143 and between STAT3 and miR-21. (C) Binding site between miR-143 and PTPN11 and between miR-21 and STAT3, IL1B.

Discussion

Peripheral arterial occlusive disease (PAOD) has become more frequent due to an ageing population, and often leads to distal limb ischemia that results in reduced quality of life and death. Furthermore, PAOD is an occult condition that goes undetected in the early stages. There is no effective cure at present, and traditional open surgery and endovascular treatment has no retarding effect on the lesion. Therefore, the underlying pathological mechanism of PAOD, and the potential pathways and core genes, have gained considerable attention in recent years with the aim of identifying potential therapeutic targets. In our study, we discovered several putative PAOD-associated genes and pathways, and built miRNA-target gene and TF-miRNA-target gene regulatory networks. The inflammatory response, immune response, cell adhesion change, JAK-STAT signaling pathway and chemokine-mediated signaling pathway were significantly associated with PAOD progression. Furthermore, STAT3, IL6, IL1B and CXCL12, and the STAT3-miRNA-143/PTPN11, miR-143/CXCL12 and STAT3-miRNA-21/IL1B regulatory pairs were likely involved in PAOD development and progression.

The major risk factors of PAOD include arteriosclerosis occlusion, thromboangiitis obliterans and multiple arteritis.[ The immune-inflammatory responses and cell adhesion play important roles in all stages of arteriosclerosis,[ and targeting either can prevent and attenuate disease development and progression.[ In this study, inflammation, cell adhesion, and immune-related BP terms and pathways were significantly enriched in the PDAO-related DEGs, indicating a similar pathological mechanism as arteriosclerosis and the possibility of an anti-inflammatory therapeutic strategy. Consistent with our findings, Tan et al reported that the DEGs between advanced and early carotid atherosclerotic plaques were also associated with the immune system, chemokine signaling pathway and focal adhesion.[ Zhang et al found that the chemokine signaling pathway, focal adhesion and JAK-STAT signaling pathway were enriched in the DEGs between coronary artery disease (CAD) and healthy tissue samples.[ In an animal model of atherosclerosis also, these pathways were significantly enriched,[ further underscoring their therapeutic potential.

Important hub genes in the PPI network of the PAOD-related DEGs included those encoding for the pro-inflammatory cytokines IL-6 and IL-1B, the chemokine CXCL12, and the transcription factor STAT3. Polymorphisms in the IL6 promoter region are correlated to the plasma levels of atherogenic markers like fibrinogen, high sensitivity C-reactive protein, apolipoprotein A1 and high density lipoprotein (HDL) cholesterol, all of which are risk factors of CAD and atherosclerosis.[ IL-6 gene polymorphisms are independent risk factors in PAOD as well, with the GG genotype of the 174 locus present at a higher frequency in PAOD patients compared to the controls, and associated with greater disease severity.[ Similarly, the 511C/T polymorphism in IL-1B gene is associated with a greater risk of CAD among the Chinese.[

CXCL12 regulates the function of several immune and inflammatory cells by interacting with its cognate receptor, and its levels are increased in the sera and atherosclerotic lesions of patients with coronary artery occlusion.[ Döring et al reported an athero-protective role of CXCL12, and found that CXCL12/ CXCR4 signaling recruited the endothelial progenitor cells and plaque-stabilizing vascular smooth muscle progenitor cells to the lesions. However, it also activates inflammatory and immune cells, which are known to promote atherosclerosis.[ In this study, CXCL12 was up-regulated in PAOD, and its exact role likely depended on the lesion area and the severity of atherosclerosis. In addition, CXCL12 was predicted as a target of miR143, which is correlated to postoperative stent restenosis in lower extremity AOD.[ Therefore, the role and potential mechanism of miR-143/ CXCL12 axis in PAOD needs further study.

STAT3 is a pleiotropic TF that can be activated by various cytokines. It is downregulated in the atherosclerotic lesions of the ApoE-/- mice, and its forced overexpression decrease inflammation, lipid accumulation and vascular smooth muscle cell proliferation, indicating an athero-protective function.[ However, phosphorylated STAT3 and JAK2 aggravated IL-1B-induced adventitial inflammation in a rat model of atherosclerosis, resulting in intimal proliferation.[ Pravastatin treatment decreased the levels of phosphorylated STAT3 in the ApoE-/- knockout mice, suggesting that STAT3 activation promotes atherosclerotic progression.[ Taken together, STAT3 has a complex role in atherosclerosis development and progression. In our study, STAT3 expression was down-regulated and the JAK-STAT signaling pathway was also enriched, indicating a protective role in PAOD. Another study showed that STAT3 activation and the vascular endothelial growth factor receptor 1 (VEGFR1)-STAT3 signaling pathway were inhibited both in human and murine experimental peripheral arterial disease. In contrast, activation of STAT3 and the endothelial growth factor receptor 1 (EGFR1)-STAT3 pathway increased perfusion in these mice, indicating that STAT3 and its signaling pathways are potential therapeutic targets in PAOD.[ STAT3 was the target of 15 putative miRNAs, including miRNA-21 which is increased during vascular restenosis in lower extremity AOD, and is a diagnostic marker of the same.[ Interestingly, in the DETF-miRNA-target gene network, miRNA-21 and miRNA-143 were also predicted as target miRNAs of STAT3 and IL-1B, while PTPN11 was the target gene of both miRNAs. Chen et al reported that STAT3 upregulated miRNA-21 by binding to its promoter region, and promoted AngII-induced angiogenesis in human microvascular endothelial cells.[ In addition, serum levels of miRNA-143 are down-regulated in patients with atherosclerosis, CAD and PAOD, and predictive of in-stent restenosis for PAOD patients.[ Furthermore, Bai et al reported that STAT3 up-regulation induced miRNA-143 while its silencing had the opposite effect.[ PTPN11 is involved in low-density lipoprotein cholesterol metabolism and CAD progression,[ although its role in PAOD or atherosclerosis is unknown. Taken together, the STAT3/miRNA-21/IL1B and STAT3/miRNA-143/PTPN11 axes play important roles in PAOD progression.

There are several limitations in our study that should be emphasized. First, our conclusions were based on in silico data, and need to be further verified in experimental studies. In addition, the DEGs were obtained from text mining and gene expression profiles, which are derived from small samples, resulting in possible bias.

Conclusions

Our findings indicate that inflammation and immune responses, and chemokine-mediated signaling are pivotal in PAOD development and progression. The miR-143/CXCL12 axis and the STAT3/miR-21/miR-143/target gene axis also appear to be pathologically relevant, and should be studied further to determine their therapeutic utility.

Ethics approval and consent to participate

Approval by an ethics committee was not necessary since the datasets included in the current study were downloaded from public databases.

Acknowledgments

The authors thanked the contributors of GSE57691 dataset for sharing the gene expression data of aortic specimens from AOD and control donors on open access.

Author contributions

Renbing Bi conceived and instructed the work. Renbing Bi, Xuwen Guan, Xiaoyan Yang checked the associated database and analyze raw data. Xuwen Guan wrote and revised the manuscript. Chunming Wang provided fund support. All of the authors read and approved the final manuscript.