Biomarkers Associated with Immune Checkpoint, N6-Methyladenosine, and Ferroptosis in Patients with Restenosis

Introduction

Coronary heart disease (CHD), characterized by unstable atherosclerotic lesions, coronary spasms, and coronary thrombosis, can trigger vessel narrowing and occlusion, eventually resulting in insufficient blood perfusion.¹ Methods such as percutaneous coronary intervention (PCI) and dual antiplatelet therapy have significantly improved the treatment outcome of patients with myocardial ischemia.^2–4 However, approximately 20–40% of patients still experience in-stent restenosis between 3 months and 4 years after implantation.^2,5 Restenosis is a complex and multifactorial disease with genetic and environmental causes.⁶ Research has shown that male gender, stent length, elevated levels of Monocyte to HDL ratio (MHR) and Neutrophil to lymphocyte ratio (NLR) were independent predictors of stent restenosis (SR) in STEMI patients treated with primary PCI with bare-metal stent (BMS) implantation.⁷ In addition, risk factors for restenosis, such as smoking, high serum cholesterol, and poor dietary habits, can only predict the likelihood of restenosis but not further restenosis of intravascular stents.⁸ Additionally, early prediction of restenosis is the first step to improving the clinical management of patients after PCI; however, there are no effective biomarkers for this. Therefore, a better understanding of the mechanisms underlying restenosis is required.

In recent years, research on microarray technology and bioinformatics analysis has increased. By analyzing differential gene expressions in patients with a given disease, biomarkers for diagnosis and prognosis can be identified.^9–12 However, very little research has been conducted on the identification of restenosis-related genes and the prediction of diagnostic markers in patients with restenosis using bioinformatics analysis.¹³ Some studies have found that high-sensitivity C-reactive proteins can predict increased risk of restenosis in patients with coronary artery disease;¹⁴ others have shown that preoperative inflammatory factors such as interleukin (IL)-6, IL-8, tumor necrosis factor-alpha, and IL-23 predict restenosis risk.^15,16 The chief pathological features of restenosis are injury to arterial endothelial cells and proliferation of new tissue, which re-shrinks the coronary artery after stenting.^17,18 The condition occurs in 30% of patients with bare-metal stents.¹⁹ While the development of drug-coated stents in recent years has lowered restenosis incidence to ~25% within one year of implantation, this rate is still relatively high.²⁰ Therefore, the prevention of restenosis after stenting remains a major problem that needs resolving.

In this study, we demonstrated the potential diagnostic markers of restenosis after stent implantation by differential expression analysis of the microarray data from the National Center for Biotechnology Information (NCBI).²¹ Further GO and KEGG analysis of these DEGs determined the signal pathways relevant to restenosis.^22,23 In addition, the protein interactions of these DEGs were studied using the Search Tool for the Retrieval of Interacting Genes (STRING).²⁴ Three genes of diagnostic significance, EEF1D, RPL36, and RPSA, were screened using MCC and MCODE algorithms. Finally, we demonstrated that restenosis was associated with immune checkpoints, ferroptosis, and N6 methyladenosine (m6A) modification.^25–29 This work provides new insights into the molecular mechanisms underlying the development of in-stent restenosis.

Materials and Methods

Microarray Data

Microarray data (GSE46560 and GSE48060) were downloaded from the GEO database³⁰ as MINIML files. GSE46560 was based on the GPL15207 platform of the Affymetrix Human Gene Expression Array. GSE48060 data was based on the GPL570 platform of [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Arrawhile. The GSE46560 dataset included five restenosis and six non-restenosis samples collected from the peripheral blood. In contrast, the GSE48060 dataset included five restenosis and 26 non-restenosis collected from peripheral blood. The extracted data were normalized by log2 transformation. The microarray data were normalized by the normalized quantiles function of the preprocessCore package in R software (version 3.4.1). Probes were converted to gene symbols according to the annotation information of the normalized data in the platform. Probes matching multiple genes were removed from these datasets, The average expression value of the gene measured by multiple probes was calculated as the final expression value. As in different datasets, we extracted multiple data sets with common gene symbols, and marking different datasets as different batches, used the removeBatchEffect function of the limma package in the R software to remove batch effects. The result of the data preprocessing was assessed by boxplot. The PCA plot was drawn to illustrate the samples before and after the batch effect.^21,30,31

Filtering DEGs

The limma package in R was used to identify differential mRNA expression between 10 restenosis and 32 non-restenosis samples from the GSE46560 and the GSE48060 datasets. The P-value was adjusted to correct for false-positive results in the GEO datasets. Adjusted P < 0.05 and |log2 fold-change (FC)| > 1.3 were defined as the statistical criteria for differential expression of mRNAs. Heat maps were implemented using the R package function pheatmap.

Functional Enrichment Analysis

Gene Ontology (GO) is a widely used tool for the functional annotation of genes, especially those involved in molecular function (MF), biological pathways (BPs), and cellular components (CCs). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis is valuable for determining high-level gene functions, such as involvement in biological pathways. The ClusterProfiler package in Bioconductor was used to obtain the functional profiles of the genes, and the ClusterProfiler package (version 3.18.0) in R was employed to visualize the results of GO and KEGG functional analyses.³⁰

Candidate Diagnostic Biomarker Screening

Search Tool for the Retrieval of Interacting Genes (STRING) is an online biological database that facilitates gene analysis and the construction of networks of gene interactions at the protein level.²⁴ The protein-protein interaction (PPI) network was constructed with STRING and the network was visualized with Cytoscape. Then, we used both MCC and MCODE algorithms to calculate the key genes.²⁵

Diagnostic Value of Characteristic Biomarkers in Restenosis

To test the predictive value of the identified biomarkers, we used the GLM function in the R (version 3.6.3) package to build a logistics model and used the GGploT2 package to visualize the results. Receiver operating characteristic (ROC) curves were generated using the mRNA expression data from the GSE46560 and GSE48060 datasets. There were ten patients with restenosis and thirty-two patients without restenosis. The diagnostic value of the identified hub genes was evaluated using the area under the ROC curve (AUC), which was between 0.5 and 1. The closer the AUC value is to 1, the better the diagnostic effect. The AUC between 0.5 and 0.7 represents predictive significance, while AUC ranging from 0.7 to 0.9 means more accuracy. When the AUC value is larger than 0.9, the accuracy is the highest.

Ferroptosis, m6A and Immune-Checkpoint Related Gene Expression in Restenosis

On the basis of previous studies, the ferroptosis-, immune-checkpoint-, and m6A-related genes were verified. We downloaded data from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) in a MINIML format. To obtain the ferroptosis-related genes, we used the systematic analysis of the aberrances and functional implications of ferroptosis in cancer published by Liu et al³² Moreover, m6A-related genes and immune checkpoint-related genes were obtained via the molecular characterization and analysis of the clinical significance of m6A modulators across 33 cancer types published by Xu et al and Zeng et al, respectively.³³ Extracting the expression of immune checkpoint genes, observing the expression value of the immune checkpoint-related genes.^34,35 We used the R package ggord to draw principal component analysis plots, and the R package ggplot2 to construct box plots. All the above analysis methods and R packages were implemented using R 4.0.3 (R Foundation for Statistical Computing 2020).

General Characteristics of Patients with Restenosis

We randomly selected 24 patients with restenosis admitted to the Department of Cardiology, the Second Affiliated Hospital of Harbin Medical University as the study subjects, and 14 patients with coronary heart disease with matched gender and age were randomly selected as the non-restenosis group. All patients with restenosis were treated regularly with aspirin. We also detailed the patients’ other concomitant diseases in the Supplementary Material. Patients with other diseases were mainly treated with commonly used drugs such as metformin or nifedipine. The clinical data of the patients was shown in Supplementary Table 1. The diagnostic criterion for restenosis was as follows: during coronary angiography reexamination, the internal diameter of the stent segment, proximal, and distal stents within 5 mm was narrowed again by over 50%. This study was approved by the Ethics Committee of the Second Affiliated Hospital of Harbin Medical University (KY2022-072). The experiment was in accordance with the principles in the Declaration of Helsinki.

Reverse Transcription-Quantitative Polymerase Chain Reaction (qRT-PCR) Verification

We detected the variations of gene expression levels in the blood of patients with and without restenosis with qRT-PCR using specific primers (Supplementary Table 2). This experiment was used to verify the differentially expressed genes. The experiment was performed with biological triplicates. The reaction system was 10 µL. β-actin was used as the internal control. Quantitative analysis of differential expression was assessed using the SYBR green reaction system on a Bio-Rad (USA) PCR machine. The relative expression levels of mRNAs were calculated by the 2-ΔΔCt method. Graphpad prism7 was used for statistical analysis, the t-test was used for the samples conforming to the normal distribution, and the non-parametric test was used for statistical analysis for the samples not conforming to the normal distribution.

Results

Identification of DEGs

We identified DEGs from the GSE46560 and GSE48060 datasets using the limma quartile normalization of quantitative traits approach with background correction. This analysis revealed 305 DEGs in the restenosis group, comprising 243 downregulated and 62 upregulated genes (Figures 1A and B).

Figure 1 Volcano plots and heat maps of differentially expressed genes in restenosis. Notes: (A) Volcano plots were constructed using fold-change and adjusted P values. The red point in the plot represents the significantly overexpressed mRNAs and the blue point indicates the significantly underexpressed mRNAs in restenosis. (B) Hierarchical clustering analysis heatmap of those mRNAs which were differentially expressed between patients with restenosis and those without restenosis. G1 represents the restenosis group and G2 represents the non-restenosis group. The change in gene expression, denoted in blue to red, increased in the restenosis group.

Functional Correlation Analysis

Using the “clusterProfiler” package in Bioconductor to obtain the functional profiles of genes via GO and KEGG pathway enrichment analysis, we found that the DEGs were mainly concentrated in the following functional categories: ribosome, coronavirus disease-COVID-19, base excision repair, RNA catabolic process, heart looping, embryonic heart tube development, and determination of bilateral symmetry (Figures 2A–D).

Figure 2 KEGG and GO enrichment analysis. Notes: The abscissa indicates gene ratios, and the enriched pathways according to the gene ontology (GO) analysis of the potential targets of the mRNAs are presented in the ordinate. (A) The figure shows the enrichment results of the KEGG pathway analysis of the differentially upregulated genes. (B) The figure shows the enrichment results of the GO pathway analysis of the differentially upregulated genes. (C) The figure shows the enrichment results of the KEGG pathway analysis of the differentially downregulated genes. (D) The figure shows the enrichment results of the GO pathway analysis of the differentially downregulated genes. For the enrichment results, a P <0.05 or FDR <0.05 was considered to indicate enrichment in a meaningful pathway.

Identification and Validation of the Diagnostic Biomarkers

To further explore the central genes related to restenosis and their mechanisms, 62 genes with upregulated expression among the 305 DEGs in the restenosis group were determined and uploaded to the STRING online database to build a PPI network; a PPI network with 74 genes as nodes and 18 edges was generated (Supplementary Figure 1). Nodes represent DEGs enriched in the STRING database, whereas edges reflect the interactions between the DEGs. Since genes with a high binding degree and high clustering coefficient are more important in maintaining the stability of the entire network, we searched for genes with a high binding degree and a clustering coefficient greater than 0.4 throughout the PPI network. The average node degree was 0.486, the average local clustering coefficient was 0.288, and the P-value of PPI enrichment was 0.458. Among the 74 nodes, the top three genes with a high binding degree were determined using Cytoscape (version 3.8.2) MCC (Figure 3A) and MCODE (Figure 3B), as follows: EEF1D, RPL36, RPSA (Figure 3C). These three genes were thus identified as central genes that play key roles in restenosis.

Figure 3 The key genes identified by MCC and MCODE. Notes: (A) Figure represents the 10 key genes calculated by MCC algorithm in cytoHubba; the darker the color, the more critical the gene. (B) Figure denotes the gene related to the module with the highest score (3.000) derived using the Cytoscape plugin MCODE. (C) The intersection of the key genes determined using MCC and MCODE is visualized with a Venn diagram (Purple represents the MCC algorithm and red represents the MCODE algorithm).

Diagnostic Effectiveness of the Biomarkers for Restenosis

The three biomarkers used to distinguish restenosis samples from control samples demonstrated high diagnostic values (Figure 4). With respect to the predictive capacity, the AUC value was 0.700 for EEF1D (95% CI = 0.472–0.928), 0.625 for RPL36 (95% CI = 0.384–0.866), 0.695 for RPSA (95% CI = 0.467–0.923).

Figure 4 Receiver operating characteristic (ROC) curve of differentially expressed genes related to restenosis, Independence. Abbreviations: TPR, true positive rate; FPR, false negative rate. Note: Each point on the curve corresponds to FPR and TPR at different thresholds.

Restenosis with m6A, Immune Checkpoint and Ferroptosis

Studies have shown that the mechanisms of ferroptosis are closely linked to the development of cardiovascular disease. Therefore, effective regulation of the ferroptosis process in cardiomyocytes has the potential to provide new action targets for the clinical prevention and treatment of restenosis. We selected 24 genes related to ferroptosis and observed variations in their expression between patients with and without restenosis (Figure 5A). GLS2 expression was closely related to ferroptosis in patients with restenosis (Figure 5B). Many regulators participate in RNA methylation, including methyltransferases (writers), RNA-binding proteins (readers), and demethylases (erasers).³⁵ Therefore, we selected genes related to these three regulator types and investigated their relationship with restenosis (Figure 6A). We found that the expression level of YTH N6-methyladenosine RNA binding protein 2 (YTHDF2) increased significantly in patients with restenosis (P < 0.01), suggesting that it may play a role in the mechanism of restenosis (Figure 6B). In addition, we also found a significant correlation between CTLA4 and restenosis among eight immune checkpoint-related genes evaluated.^34,35 (Figure 7A and B).

Figure 5 The expression distribution of the ferroptosis-related mRNAs in the restenosis and non-restenosis groups. Notes: (A) Ferroptosis-related gene analysis results. (B) The expression distribution of GLS2 gene in restenosis, where the horizontal axis represents restenosis and non-restenosis samples, the vertical axis represents the gene expression distribution. *P < 0.05, the asterisk represents the degree of importance (*P). The significance of the two groups of samples passed the Wilcox test.

Figure 6 The expression distribution of the m6A-related mRNAs in the restenosis and non-restenosis groups. Notes: (A) m6A-related gene analysis results. (B) The expression distribution of YTHDF2 gene in restenosis, where the horizontal axis represents restenosis and non-restenosis samples, the vertical axis represents the gene expression distribution. *P < 0.05, the asterisk represents the degree of importance (*P). The significance of the two groups of samples passed the Wilcox test.

Figure 7 The expression distribution of the immune checkpoint-related mRNAs in the restenosis and non-restenosis groups. Notes: (A) Immune checkpoint-related gene analysis results. (B) The expression distribution of CTLA4 gene in restenosis, where the horizontal axis represents restenosis and non-restenosis samples, the vertical axis represents the gene expression distribution. **P < 0.01, the asterisk represents the degree of importance (*P). The significance of the two groups of samples passed the Wilcox test.

Validation of Diagnostic Genes

We verified the identified DEGs using qRT-PCR in peripheral blood samples of patients with restenosis and controls (Figures 8A–F). Among the evaluated genes, the expression levels of RPSA and RPL36 were up-regulated in restenosis patients (Figures 8E and F). Notably, the expression level of RPSA was significantly increased in the patients (Figure 8F). Finally, we built a flow chart for the whole experiment design (Figure 9).

Figure 8 qRT-PCR validation of gene expression. Abbreviation: ns, not significant. Notes: Quantitative real-time PCR (qRT-PCR) of (A) GLS2, (B) CTLA4, (C) EEF1D, (D) YHDF2, (E) RPL36 and (F) RPSA in the blood of restenosis compared to the control. *P < 0.05, **P < 0.002.

Figure 9 Flow chart of restenosis. Note: Figure was created using the open tool, Figdraw.

Discussion

Restenosis is characterized by a gradual narrowing of coronary artery in-stent lesions,³⁶ and approximately 30% of patients with CAD undergoing PCI experience restenosis.^2,37 Restenosis usually occurs in a few months after PCI, which is often accompanied by arterial damage and new intima formation. In recent years, balloon dilatation and drug-eluting stents have reduced the restenosis incidence and are thus widely used; however, approximately 20% of restenosis cases still occur within 2 years after DES or PCI.³⁸ Additionally, these two methods do not effectively inhibit the proliferation and migration of vascular smooth muscle cells, nor do they induce adequate reendothelialization in blood vessels.^39–41 Thus, restenosis remains a major challenge in clinical practice, and developing potential biomarkers to reduce restenosis incidence is necessary. Fortunately, gene-chip technology has allowed for the increased use of gene expression data in clinical studies, and microarrays are applied extensively in peripheral blood studies and cardiopathy.^42–45 For example, studies have shown that the higher the sST2 plasma levels, the lower the expression of ST2L on the membrane of the cells composing the carotid atherosclerotic plaques.⁴⁶ Although several studies have examined factors related to AMI or CHD occurrence, for instance, some research report that ECG index can predict Diastolic dysfunction (DD) in patients without known coronary artery disease. Such as the electrocardiographic diastolic index (EDI), which is an inexpensive, feasible, and easy to use formula, appears to have a considerable role in predicting diastolic dysfunction (DD) in adult patients,⁴⁷ little is known about the predictive value of in-stent restenosis after AMI.⁴⁸ Therefore, it is necessary to prevent restenosis occurrence and explore the underlying molecular mechanism. Here, we used GEO gene expression datasets to test restenosis-related differential gene expression and identify DEGs via functional analyses. Next, we employed MCC and MCODE calculation methods to screen out genes as potential diagnostic markers. Subsequently, we verified the relationship between restenosis-related genes with ferroptosis, m6A, and the immune checkpoint. Our novel study should contribute to the better diagnosis and treatment of vascular restenosis. Additionally, the findings may provide a definite biomarker to predict restenosis risk in patients undergoing PCI or DES.

In the enrichment analysis of differentially expressed genes, we found fewer enrichment pathways for genes associated with restenosis. We believe that the cause of this problem is due to the pathological mechanism of RNA-regulated restenosis. It is well known that the pathological mechanisms of RNA-regulated restenosis include the following two main aspects: first, neointimal proliferation, including smooth muscle migration and proliferation, and extracellular matrix deposition; second, vascular remodeling. For the former, infiltration of various inflammatory cells is predominant, whereas the current mechanism of vascular remodeling leading to restenosis is not dominated by inflammatory infiltration. Therefore, when performing differential expression analysis, we took into account the diversity of mechanisms by which it occurs, and the absence of KEGG- and GO-related enrichment pathways would not affect the role of these genes in restenosis. In our analysis, this result may be related to the independent role of these DEGs in the pathogenesis of restenosis and may be a novel direction for the mechanism of restenosis occurrence. In the analysis of clinical patient information. We found that patients with restenosis had a higher probability of having hypertension, diabetes mellitus, and a higher history of smoking. However, we believe that this is a common risk factor for coronary artery disease and cannot be a specific risk factor for restenosis.

The eukaryotic translation extension factor-1 (EEF1) protein complex is composed of non-ribosomal protease factors; among these members, EEF1D plays a carcinogenic role by maintaining the Akt-mTOR and Akt-Bad signaling pathways in osteosarcoma.⁴⁹ The RPSA gene encodes SA ribosomal protein, which has a variety of functions, including laminin binding, ribosomal function, nuclear function, extracellular matrix modification, and signal transduction.⁵⁰ These previous studies coincide with the results of our current study. Metabolism is a strong predictor of inflammation and injury. In recent years, there have been increasing reports of a metabolic disorder involving intracellular lipid peroxidation that produces toxic lipids and triggers cell death, which is termed ferroptosis. The ferroptosis-related gene GLS2 was among the restenosis-associated genes identified in this study. Moreover, studies have confirmed that the m6A modification significantly affects the occurrence and development of restenosis.³⁹ But our results revealed that the expression level of the methylation-related gene YTHDF2 was not increased in patients with restenosis. YTHDF2 is an m6A reader that binds to m6A-containing mRNAs to regulate the localization and stability of the bound mRNAs.⁵¹ Moreover, YTHDF2 regulates the occurrence and development of several tumors, including breast and prostate cancers.^52–57 Therefore, methylation of YTHDF2 may play a pathophysiological role in restenosis through structural changes of it. Among our results of immune checkpoint genes, CTLA4 was also not significantly associated with restenosis.

Studies have shown that the development of restenosis is associated with damage to endothelial and smooth muscle cells. Smooth muscle cells are relevant to maintaining the integrity of the heart muscle. PCI damages the endothelial cells, destroying them as a whole and exposing the subendothelial matrix to the blood. In turn, vascular smooth muscle cells continue to proliferate and the endothelium overgrows, thereby increasing the risk of restenosis. Indeed, laminin is a family of extracellular matrix glycoproteins. Involved in cell adhesion, differentiation, migration, signaling, neurite growth and metastasis, RPSA is the cell surface receptor for laminins and plays an important role in the adhesion of cells to the basement membrane and the subsequent activation of signaling pathways.⁵⁸ The gene encoding RPL36 is located on position 19p13.2 at the short arm of human chromosome 19 and consists of 104 amino acids with a molecular weight of 12128 daltons.⁵⁹ It may therefore also play a role in cell fate determination and tissue morphogenesis. Thus, RPSA and RPL36 play a key role in the development of restenosis.

Our findings from public databases and bioinformatics tools were thus verified in clinical samples, demonstrating that during coronary artery stenting, a significant increase in RPSA and RPL36 expression could predict the occurrence of in-stent restenosis. However, it is interesting to note that the significant changes in the expression levels of EEF1D, GLS2, CTLA4, and YTHDF2 in the restenosis patients in our clinical sample did not match the results of bioinformatic analysis. One potential reason for this discrepancy could be that these genes may be associated with a particular phenotype in the restenosis-associated changes, and thus did not cause the change of expression quantity. This study also has some limitations. First, this study is a retrospective study, and some prospective studies are needed in the future to further diagnose the occurrence of restenosis. Second, the total number of restenosis samples analyzed in this work was small. Thus, our results need to be further confirmed in experiments with larger sample sizes. Third, the ferroptosis, m6A, and immune checkpoint-related genes involved in this experiment are only a preliminary screen, and further experiments are still needed for validation.

Conclusion

In conclusion, we identified RPL36 and RPSA as the diagnostic biomarkers for restenosis. This study can contribute to the timely diagnosis and improved treatment of restenosis.