BPI and KIR6.1 as significant hub genes for vein graft restenosis.

BACKGROUND: Vein graft restenosis (VGR), which appears to be caused by dyslipidemia following vascular transplantation, seriously affects the prognosis and long-term quality of life of patients.
METHODS: This study analyzed the genetic data of restenosis (VGR group) and non-stenosis (control group) vessels from patients with coronary heart disease post-vascular transplantation and identified hub genes that might be responsible for its occurrence. GSE110398 was downloaded from the Gene Expression Omnibus database. A repeatability test for the GSE110398 dataset was performed using R language. This included the identification of differentially expressed genes (DEGs), enrichment analysis via Metascape software, pathway enrichment analysis, and construction of a protein-protein interaction network and a hub gene network.
RESULTS: Twenty-four DEGs were identified between VGR and control groups. The four most important hub genes (KIR6.1, PCLP1, EDNRB, and BPI) were identified, and Pearson's correlation coefficient showed that KIR6.1 and BPI were significantly correlated with VGR. KIR6.1 could also sensitively predict VGR (0.9 < area under the curve ≤1).
CONCLUSION: BPI and KIR6.1 were differentially expressed in vessels with and without stenosis after vascular transplantation, suggesting that these genes or their encoded proteins may be involved in the occurrence of VGR.

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 Oryctolagus cuniculus (Rabbit) Oligo Microarray (Platform GPL13288; Agilent, Santa Clara, CA, USA).

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 Homo sapiens Entrez gene IDs using the latest version of the database (last updated on 2019-03-12). If multiple identifiers corresponded to the same Entrez gene ID, they were considered to be a single Entrez gene ID in the downstream analyses. Each gene list was assigned a unique color, which was used throughout the subsequent analysis.

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).

Figure 1.

(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.

DEG, differentially expressed gene.

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).

Figure 2.

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

Figure 3.

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.

Figure 4.

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.

Figure 5.

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.

Figure 6.

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).

MCODE, Molecular Complex Detection algorithm.

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: IL1A, BPI, PCLP1, SAA3, KIR6.1, KCNS3, TNF, EDNRB, EZR, and ABCG2 (Figure 7c). The four most significant hub genes (KIR6.1, PCLP1, EDNRB, and BPI; Figure 7d) were present in the nine common genes in Figure 7a, and the hub genes identified in Figure 7c.). The functions of the four most significant hub genes are shown in Table 1.

Figure 7.

(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.

Table 1.

Functions of the four key hub genes.

Gene symbol	Full name	Function
KIR6.1	Potassium voltage-gated channel subfamily J member 8 (KCNJ8)	Vascular smooth muscle contraction and inwardly rectifying K+ channels. Related gene ontology annotations include inward rectifier potassium channel activity and ATP-activated inward rectifier potassium channel activity.
PCLP1	Podocalyxin-like (PODXL)	Regulation of adhesion, cell morphology, and cancer progression. Functions in cancer development and aggressiveness by inducing cell migration and invasion through its interaction with the actin-binding protein EZR. Affects EZR-dependent signaling events, leading to increased activities of mitogen-activated protein kinase and phosphoinositide 3-kinase pathways in cancer cells.
EDNRB	Endothelin receptor type B	Member of the endothelin receptor group of G-protein-coupled receptors that also includes ETA. Located primarily in vascular endothelial cells where it functions in vasoconstriction, vasodilation, and cell proliferation.
BPI	Bactericidal/permeability-increasing protein	Lipopolysaccharide binding protein associated with human neutrophil granules. Has antimicrobial activity against Gram-negative organisms.

Expression analysis of the most significant hub genes

Compared with the Sham group, the expression of KIR6.1, PCLP1, and EDNRB was significantly lower, and the expression of BPI was significantly higher in the 7 Day group (p < 0.05; Figure 8a). The correlations between the four most significant hub genes are shown in Figure 8b.

Figure 8.

(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 KIR6.1 (ρ = −0.630, p = 0.009), PCLP1 (ρ = −0.988, p < 0.001), EDNRB (ρ = −0.895, p < 0.001), and BPI (ρ = −0.081, p = 0.765) were all significantly correlated with VGR (Table 2). After holding all other variables at a fixed value, the multivariate linear regression model showed that VGR remained significantly associated with the expression of PCLP1, EDNRB, and BPI (p < 0.001) (Table 2). To identify an accurate threshold for the hub genes most sensitively able to predict VGR, we constructed ROC curves. The expression of KIR6.1, PCLP1, and EDNRB was associated with a diagnosis of VGR (0.9 < area under the curve ≤1, p ≤ 0.0001) (Table 3).

Table 2.

Correlation and linear regression analyses between VGR and relevant gene expression.

Gene symbol	VGR
	Pearson’s correlation coefficient		Multiple linear regression
	ρa	p-value	βb	p-value
KIR6.1	−0.630	0.009*	−0.039	0.195
PCLP1	−0.988	<0.001*	−0.098	<0.001*
EDNRB	−0.895	<0.001*	−0.133	<0.001*
BPI	−0.081	0.765	−0.076	<0.001*

Pearson’s correlation coefficient between VGR and relevant characteristics; ρ: Pearson’s correlation coefficient.

Multiple linear regression analysis, β: parameter estimate.

Significant variables: p < 0.05.

VGR: vein graft restenosis.

Table 3.

Receiver operator characteristic curve analysis of key gene expression for VGR.

Gene symbol	VGR
	AUC	p-value	ODT
KIR6.1	1.000	0.001*	8.079
PCLP1	1.000	0.001*	12.361
EDNRB	1.000	0.001*	8.791
BPI	0.600	0.515	13.125

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 (p < 0.05; Figure 9).

Figure 9.

Histological patterns of control and VGR vessels by HE staining.

VGR, vein graft restenosis, HE, hematoxylin–eosin.

Expression levels of the most significant genes by histology

The immunofluorescence assay showed that the expression of KIR6.1, PCLP1, and EDNRB was significantly downregulated in VGR tissues compared with the control group (p < 0.05; Figure 10).

Figure 10.

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 KIR6.1, PCLP1, and EDNRB was significantly lower in the VGR group than in the control group (p < 0.05, Figure 11a–c). However, BPI expression was significantly upregulated in VGR tissue compared with the control group (Figure 11d).

Figure 11.

Verification of the expression of KIR6.1, PCLP1, EDNRB, and BPI at the mRNA level.

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, PCLP1, EDNRB, and BPI were differentially expressed between them. Of these genes, BPI was highly expressed in restenosis vessels, while KIR6.1, PCLP1, and EDNRB were expressed at low levels. Compared with a previous report,[5] this study offers the advantages of advanced bioinformatics analysis and molecular investigation of patients with VGR.

BPI encodes bactericidal/permeability enhancing protein which is a bidirectional, bar-shaped molecule rich in lysine N-terminals. It is the main component of neutrophils and is a member of the lipid transfer/lipopolysaccharide-binding protein gene family, which is involved in the transport of lipids and lipoproteins. Eduardo et al.[15] previously showed that high BPI expression was an independent risk factor for high levels of low-density lipoprotein and cholesterol after adjusting for age and body mass index. Consistent with this, our bioinformatics analysis showed that BPI was highly expressed in patients with restenosis. Based on the biological role of BPI, we hypothesize that BPI molecules contribute to the formation of restenosis following vascular transplantation by influencing atherosclerotic changes. As well as an indicator for blood lipid monitoring, BPI could also be used as an indicator of the risk of restenosis after vascular transplantation and a targeted therapy molecule.

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 KIR6.1 expression in restenosis patients, which is consistent with our current findings. We speculate that KIR6.1 may inhibit restenosis following vascular transplantation through these mechanisms.

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.