Abstract
Keywords
Background
Coronary heart disease is caused by coronary artery atherosclerosis and dyslipidemia, which results in the narrowing or occlusion of the lumen, leading to myocardial ischemia or necrosis. In the United States, coronary heart disease accounts for more than half of all deaths from heart disease, so the disease represents a serious threat to human health. 1 Since the introduction of early primary and secondary prevention measures, such as drug intervention and cardiac vascular transplantation, the global death rate from coronary heart disease has decreased. 1 However, coronary heart disease continues to seriously endanger the health of both middle-aged and older people, and its related complications can greatly affect their quality of life. As an example, vein graft restenosis (VGR), a re-narrowing of the vessels after bypass, can lead to ischemia of organs and tissues. 2
Microarray analysis can simultaneously capture the expression of tens of thousands of genes and explore genomic changes associated with disease initiation and progression. 3 As early as 2010, Kullo applied bioinformatics to analyze early risk assessment in patients with coronary heart disease. 4 Since the beginning of the 21st century, researchers have increasingly applied bioinformatics technology to study differentially expressed genes (DEGs) associated with disease progression, and to explain their role in biological processes, molecular functions, and signal pathways. Therefore, in this study we adopted a bioinformatics approach to screen genetic data from patients with coronary heart disease who underwent vascular transplantation, with or without VGR.
DEGs were identified, and gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to construct a protein–protein interaction (PPI) network. Finally, the predicted role of DEGs in restenosis following vascular transplantation was subjected to verification by molecular experiments. The identification of genetic and molecular variants associated with restenosis following vascular transplantation will improve our understanding of the underlying mechanism that leads to this condition and provide the basis for novel, targeted therapies.
Methods
Access to public data
One gene expression profile, GSE110398,
5
was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).
6
GSE110398 array data consist of an mRNA expression profile of 16 tissue samples, including six samples from vein grafts removed 1 day post-surgery (1 Day), six from vein grafts removed 7 days post-surgery (7 Day), and four from vein grafts obtained from a sham surgery group (Sham). Gene expression profiles were generated using Agilent-020908
Repeatability test for the GSE110398 dataset
Intra-group data repeatability was verified by performing Pearson’s correlation test. The R programming language (version: 3.6.2, https://www.r-project.org/) was used to provide the software and operating environment for all statistical analyses and to draw the graphs. Correlations between all samples from GSE110398 were visualized using heat maps, which were also drawn using R (version: 3.6.2, https://www.r-project.org/).
Identification of DEGs
GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) is an interactive online tool that was used to identify DEGs from the GEO series. 7 GEO2R was applied to identify DEGs between the 1 Day group and the Sham group, between the 7 Day group and the Sham group, and between the 7 Day group and the 1 Day group. If one probe set lacked a homologous gene, or if one gene had numerous probe sets, the data were removed. The cut-off criteria were log (FC)≥1 or log (FC) ≤–1 and P < 0.05. Venn diagrams were delineated using an online Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/), which could then be used to visualize common DEGs shared between different groups.
Enrichment analysis via Metascape software
Metascape 8 integrates several authoritative data resources such as GO, KEGG, UniProt, and DrugBank to provide pathway enrichment and biological process annotation, and gene-related protein network and drug analyses.
Gene lists
User-provided gene identifiers were first converted into their corresponding
The overlaps between these lists were shown in Circos 9 plots. Overlaps were considered based on gene functions or shared pathways, or between genes that share the same enriched ontology terms. Only ontology terms that contained fewer than 100 genes were used to calculate functional overlaps, to avoid linking genes using very general annotation.
Pathway and process enrichment analysis
To further capture the relationships between the terms, a subset of enriched terms was selected and rendered as a network plot, where terms with a similarity >0.3 were connected by edges. Terms were selected with the most significant p-value from each of the 20 clusters, with the constraint that there were no more than 15 terms per cluster and no more than 250 terms in total. The network was visualized using Cytoscape, 10 where each node represented an enriched term and was colored first by its cluster ID and then by its p-value.
If multiple gene lists were provided, the nodes were represented as pie charts, where the size of a pie was proportional to the total number of hits that fell into that specific term. The pie charts were color coded based on gene list identities, and the size of a slice represented the percentage of genes under the term that originated from the corresponding gene list. This type of plot is particularly useful for visualizing whether terms are shared by multiple lists or unique to a specific list, as well as for understanding how these terms associate with each other within the biological context of a meta-study.
PPI enrichment analysis
For each given gene list, PPI enrichment analysis was carried out using the following databases: BioGrid, 11 InWeb_IM, 12 and OmniPath. 12 The resulting network contained a subset of proteins that formed physical interactions with at least one other member in the list. If the network contained between three and 500 proteins, the Molecular Complex Detection (MCODE) algorithm 13 was applied to identify densely connected network components. MCODE networks identified for individual gene lists were then gathered.
Pathway and process enrichment analyses were independently applied to each MCODE component, and the three best-scoring terms by p-value were retained as functional descriptions of the corresponding components.
Construction and analysis of PPI and hub gene networks
The Search Tool for the Retrieval of Interacting Genes (STRING) online database (http://string-db.org) predicts and traces PPI networks once common DEGs have been imported into it. Here, the STRING database was used to construct the DEG PPI network. Once the degrees were set (degrees ≥10), hub genes were identified using the free visualization software tool, Cytoscape (version 2.8). 14 The most significant hub genes were identified by Venn diagram analysis. A hierarchical clustering heat map of significant hub gene expression was visualized using R. Correlation analysis between significant hub genes was also carried out.
Pearson’s correlation coefficient was used to analyze the correlation between VGR and expression of the most significant hub gene. A multivariable linear regression model was used to identify genes that were independently predictive of VGR. Finally, receiver operator characteristic (ROC) curve analysis was performed to determine the ability of the most significant hub genes to predict VGR. All statistical analyses were conducted using SPSS software (version 21.0; IBM). A p-value < 0.05 was considered statistically significant.
Coronary heart disease patient samples
Pathological sections from vessels were obtained from three patients who received vein graft treatment for coronary heart disease at Tianjin Chest Hospital between 2018 and 2019. Coronary heart disease was diagnosed by narrowing or occlusion of a coronary artery on a coronary angiogram. Samples from control vessels without restenosis and from vessels that underwent restenosis following a vein graft were obtained from each patient. This research conformed to the Declaration of Helsinki and was authorized by the Human Ethics and Research Ethics Committees of Tianjin Chest Hospital. Written informed consent was obtained from all participants.
The aortic valve of samples was dissected and cut into separate parts. Part of the artery was fixed in 4% poly-formaldehyde solution and processed for paraffin embedding and sectioning. The aortic valve was cut into 6-µm transverse sections and stained with hematoxylin–eosin (HE) after fixation in 4% poly-formaldehyde solution embedded in paraffin wax. Immumofluorescence techniques were then used to detect the expression levels of proteins encoded by the most significant hub genes. Briefly, the rabbit anti‑human BPI polyclonal antibody (dilution 1:3000; Abcam, Cambridge, UK) or rabbit anti‑human KIR6.1 polyclonal antibody (dilution 1:3000; Abcam) was added dropwise and incubated overnight at 4°C, washed three times with phosphate-buffered saline for 5 minutes each time, then incubated with a goat anti-rabbit monoclonal antibody conjugated with a fluorescent dye (dilution rate = 1:5000, ab205718, Abcam) at 37°C for 1 hour. 4′,6-diamidino-2-phenylindole was then added and incubated for 10 minutes at room temperature in the dark. Fluorescence was examined using a Nikon Eclipse C1 microscope (Nikon, Tokyo, Japan).
Quantitative reverse transcription (qRT)-PCR assay
Total RNA was extracted using an RNAiso Plus kit (Thermo Fisher Scientific, Waltham, MA, USA) and reverse-transcribed into cDNA using the Servicebio® RT First Strand cDNA Synthesis Kit (Servicebio, Wuhan, China) with 2 × SYBR Green qPCR Master Mix (High ROX) (Servicebio). The reaction was incubated at 55°C for 5 minutes. qRT-PCR was performed using a Light Cycler® 4800 System with specific primers for the most significant hub genes. Primers sequences were: BPI, forward: 5′-CCTTCTCAGAGCCTTACAT-3′, reverse: 5′-TCTCCCTCATCACTTTCC-3′; GAPDH, forward: 5′-ATCCGATTACCGATACCTAGACC-3′, reverse: 5′-ATGGACTATATCCGACGACGA-3′; and KIR6.1, forward: 5′-TATTATCCAGCCTACCTC-3′, reverse: 5′-ATTGCACTAACTACCCAC-3′. PCR cycling conditions were 95°C for 10 minutes then 40 cycles of 95°C for 15 s and 60°C for 1 minute. A melt curve was generated from 60°C to 95°C, with a temperature rise of 0.3°C per 15 seconds. Actin was amplified as an endogenous control.
Results
Validation of public data
Pearson’s correlation test was used to validate the repeatability of intra-group data. Strong relationships were identified in the GSE110398 dataset among individuals in the 1 Day group, the 7 Day group, and among different tissues in the Sham group (Figure 1a).

(a) Pearson’s correlation coefficient between samples. The color reflects the intensity of the correlation. 0< correlation <1 represents a positive correlation; –1< correlation <0 represents a negative correlation. A larger absolute value of a number represents a stronger correlation. (b) Volcano plot of DEGs between the 1 Day group and the Sham group. (c) Volcano plot of DEGs between the 7 Day group and the Sham group. (d) Volcano plot of DEGs between the 7 Day group and the 1 Day group.
Identification of DEGs between different groups
A total of 771 DEGs were identified between the 1 Day group and the Sham group (Figure 1b), 24 between the 7 Day group and the Sham group (Figure 1c), and 126 between the 7 Day group and the 1 Day group (Figure 1d). By comparing gene differences between samples from vein grafts removed 1 day post-surgery and those from vein grafts removed 7 days post-surgery, the 126 more significant genes that could affect the degree of restenosis were identified.
Enrichment summary using Metascape software
The DEGs were mainly enriched in the response to wounding, the response to peptides, inorganic ion homeostasis, blood circulation, extracellular structure organization, reactive oxygen species metabolic processes, the calcium signaling pathway, the response to lipopolysaccharide, anion transport, positive regulation of cell adhesion, the AGE-RAGE signaling pathway in diabetic complications, the transport of small molecules, hematopoietic cell lineage, the cellular response to organic cyclic compounds, and the apoptotic signaling pathway (Figure 2). There were several overlaps and interactions between these DEGs, as shown by the Circos plot (Figure 3a, b). Through the network of enriched terms, DEG lists were shown to be mainly enriched in blood circulation, inorganic ion homeostasis, the response to inorganic substances, the apoptotic signaling pathway, the response to steroid hormones, the response to peptides, hematopoietic cell lineage, the regulation of secretion, extracellular structure organization, anion transport, and the response to wounding (Figure 4a); all p-values were < 0.05 (Figure 4b). Representing the network of enriched terms as pie charts showed that there were strong interactions between all DEG lists (Figure 5).

Heat map of enriched terms across input gene lists, colored by p-value.

Overlap between gene lists: (a) only at the gene level, where purple curves link identical genes; (b) including the shared term level, where blue curves link genes that belong to the same enriched ontology term. The inner circle represents gene lists with hits arranged along the arc. Genes that hit multiple lists are shown in dark orange, and genes unique to a list are shown in light orange.

Network of enriched terms: (a) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (b) colored by p-value, where terms containing more genes tend to have a more significant p-value.

Network of enriched terms represented as pie charts, where pies are color coded based on gene list identity.
Following PPI enrichment analysis, four networks were identified: all lists merged colored by counts (full connection) (Figure 6a), all lists merged colored by counts (keep MCODE nodes only) (Figure 6b), all lists merged colored by cluster (full connection) (Figure 6c), and all lists merged colored by cluster (keep MCODE nodes only) (Figure 6d). Some gene nodes were shared in all DEG lists, and significant modules existed in the network.

Protein–protein interaction network and MCODE components identified in the gene lists. (a) All lists merged Colored by Counts (Full Connection). (b) All lists merged Colored by Counts (Keep MCODE Nodes Only). (c) All lists merged Colored by Cluster (Full Connection). (d) All lists merged Colored by Cluster (Keep MCODE Nodes Only).
Identification of the most significant hub genes
A Venn diagram was used to show the common genes between 1 Day–Sham, 7 Day–Sham, and 7 Day–1 Day. Twenty-four genes were common to both 1 Day–Sham and 7 Day–Sham, nine to both 7 Day–Sham and 7 Day–1 Day, 126 to both 1 Day–Sham and 7 Day–1 Day, and nine were common among all three groups (Figure 7a). The PPI network showed the interactions of common genes between 1 Day–Sham and 7 Day–Sham (Figure 7b), and 10 hub genes were screened from this network:

(a) Venn diagram showing genes common to 1 Day–Sham, 7 Day–Sham, and 7 Day–1 Day. (b) Protein–protein interaction network showing the interactions of genes common to 1 Day–Sham and 7 Day–Sham. (c) The 10 hub genes were screened from the protein–protein interaction network. (d) Venn diagram showing the four most significant hub genes.
Functions of the four key hub genes.
Expression analysis of the most significant hub genes
Compared with the Sham group, the expression of

(a) Heat map of the expression of the four most significant hub genes. (b) Correlation analysis between the four most significant hub genes.
Pearson’s correlation coefficient showed that
Correlation and linear regression analyses between VGR and relevant gene expression.
Pearson’s correlation coefficient between VGR and relevant characteristics; ρ: Pearson’s correlation coefficient.
Multiple linear regression analysis, β: parameter estimate.
Significant variables:
VGR: vein graft restenosis.
Receiver operator characteristic curve analysis of key gene expression for VGR.
Significant variables.
AUC: area under curve; ODT: optimal diagnostic threshold; VGR: vein graft restenosis.
Histological patterns of vessels in the control and VGR groups
HE staining revealed that the lumens of vessels in the VGR group were significantly smaller than those in the control group (

Histological patterns of control and VGR vessels by HE staining.
Expression levels of the most significant genes by histology
The immunofluorescence assay showed that the expression of

Immunofluorescence assay of KIR6.1, PCLP1, EDNRB, and BPI expression in the aortic valve.
Verification of mRNA expression of the most significant hub genes
The qRT-PCR assay showed that the mRNA expression of

Verification of the expression of
Discussion
This study used a bioassay to compare gene expression between restenosis and non-restenosis vessels from patients with coronary heart disease who underwent vascular transplantation, and found that
KIR6.1 is an isomer of the KIR family that can regulate the activity of inward potassium channels and maintain the constriction of vascular smooth muscle. Du et al. found that KIR6.1/K-ATP inhibits the p38 mitogen-activated protein kinase–nuclear factor–kappa-B signaling pathway, transforming the microglial cell phenotype from harmful M1 to beneficial M2; therefore, this may offer a therapeutic target for Parkinson’s disease.
16
Moreover, Diehlmann et al. showed that KIR6.1 inhibited the adipogenic differentiation of undifferentiated mesenchymal stem cells.
17
However, there have been few studies on KIR6.1 in restenosis after vascular transplantation. Ribalet et al. reported that KIR6.1 affects K-ATP channels and indirectly participates in purine metabolism. When KIR6.1 expression was reduced, K-ATP channels were shown to become hyposensitive while purine production declined, which affected the regulation of vascular remodeling. Low KIR6.1 expression also causes the dysfunction of vascular smooth muscle contraction. For example, blood vessel over-dilation can break and degrade the middle artery, resulting in atherosclerotic changes in healthy arteries.18,19 In the present study, our bioinformatics analysis detected low
Liu et al. previously identified potential hub genes associated with vein graft restenosis, but did not validate their findings by experimental methods. 5 Our study not only identified the more significant genes affecting the extent of restenosis by identifying expression differences between samples from vein grafts removed 1 day post-surgery with those from vein grafts removed 7 days post-surgery, but also verified the role of these genes using experimental methods. However, our study was limited by its use of tissue samples. Although KIR6.1 might be a tissue biomarker for vein graft restenosis, it might not be a representative blood biomarker. Further study is needed to confirm our findings, and clinical trials are required to verify the role of the identified hub genes in vascular transplantation restenosis.
Conclusions
Bioinformatics analysis provides evidence for investigating the mechanism of restenosis following vascular transplantation. We identified BPI and KIR6.1 as potential key factors involved in protecting against or reducing restenosis. This study provides the basis for the prevention and targeted treatment of restenosis, as well as new ideas for follow-up research.
