Identification of Potential Key Genes Involved in the Carotid Atherosclerosis

Introduction

Carotid atherosclerosis is the manifestation of systemic atherosclerosis in the carotid arteries.¹ It usually occurs in adolescence and gradually worsens with age. It is currently believed to be closely related to the occurrence of ischemic stroke in the elderly.² In the early stage of carotid atherosclerosis, it first manifests as intima-media thickening, and then gradually forms atherosclerotic plaques.³ This cause corresponding hemodynamic changes, which in turn leads to ischemic cerebrovascular events. Contrary to the steady or declining trend of most disease, the incidence of atherosclerosis has been shown to increase in both men and women.⁴ In the past decade, the number of available treatments for carotid atherosclerosis has increased. In particular, immunotherapy and targeted therapy have improved the prospects of treatment. However, the efficiency of treatment is insufficient.^5,6 Therefore, it is urgent to find therapeutic targets and deepen the understanding of carotid atherosclerosis.

Some studies have shown that the abnormal expression of some genes is closely related to the development of carotid atherosclerosis. S100A12, a calcium binding protein involved in inflammatory signaling, was also shown to be elevated in patients bearing carotid atherosclerotic lesions.⁷ Serum levels of TNFα and L-selectin, a lectin type adhesion molecule expressed on leukocytes surface, were shown to associate with larger plaque size.⁸ A previous study showed that there is an association between histological features of carotid plaque instability and serum levels of circulating matrix metalloproteinase (MMP)-1, MMP-7, tissue inhibitor of matrix protease (TIMP)-1, and IL-8.⁹ But the molecule mechanism of carotid atherosclerosis has not been fully elucidated. Therefore, it is important to deepen the understanding for the relationship between carotid atherosclerosis and abnormal changes in genes.

In this study, we obtained the differentially expressed genes (DEGs) between atheroma plaque and control samples. We further explored the biological functions of DEGs and constructed a protein–protein interaction (PPI) network. In addition, we obtained 16 hub genes. And we established the carotid atherosclerosis development-related network modules. In summary, we identified of hub genes involved in the development of carotid atherosclerosis. The study may help to provide guideline for the treatment of carotid atherosclerosis.

Materials and Methods

Data Collection and Processing

The expression profiling data of hypertensive patients in our study was derived from one published dataset (GSE43292)¹⁰ available in Gene expression omnibus (GEO) database. The GSE43292 dataset based on the GPL6244 platform contains 32 atheroma plaque and 32 macroscopically intact tissue (control samples) from 32 hypertensive patients. The dataset was profiled using the Affymetrix Human Gene 1.0 ST Array (Santa Clara, CA) platform. The justRMA method in the affy package¹¹ was applied to normalize the raw data of the dataset. If one gene corresponded to multiple probes, the average expression value of these probes was considered to be the expression value of the gene. The workflow of this study is shown in Figure 1.

Figure 1 The workflow of this study. Abbreviation: GEO, gene expression omnibus.

Analysis of Differentially Expressed Genes (DEGs) for Atheroma Plaque and Control Samples

The limma package¹² was used to analyze differentially expressed genes (DEGs) between atheroma plaque and control samples. The DEGs with P value (adjusted false positive rate) <0.05 and |log2 fold change (FC)| >0.5 were considered significant.

Enrichment Analysis

The clusterProfiler package¹³ in R was used to functionally analyze the enriched pathways of the DEGs, including Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. P value <0.05 was considered significant. Gene set enrichment analysis (GSEA) was performed using GSEA java software¹³ against the reference gene sets c2.cp.kegg.v6.2. symbols.gmt and c5.bp.v7.0.entrez.gmt, which were downloaded from the Molecular Signature Database (http://software.broadinstitute.org/gsea/index.jsp).¹³ The nominal (NOM) P value <0.05 and false discovery ratio (FDR) <0.25 were considered statistically significant.

Protein–Protein Interaction (PPI)

To obtain insights into the interactions among DEGs associated with carotid atherosclerosis. A PPI network was constructed using the Search Tool for the Retrieval of interacting Genes (STRING),¹⁴ a database of known and predicted protein interactions. An interaction with score>990 was included in our present study. The PPI network was visualized by Cytoscape 3.8.2 software (http://www.cytoscape.org/).¹⁵

Hub Genes Selection and Analysis

CytoHubba is a common tool for analyzing PPI networks.¹⁶ The hub genes were selected using the cytoHubba plug-in of Cytoscape software. Those overlap genes in the top 50 genes, from the four ranking methods, were defined as hub genes. The four ranking methods include Degree, density of maximum neighborhood component (DMNC), maximal clique centrality (MCC), maximum neighborhood component (MNC).

Establishment of the Carotid Atherosclerosis-Related Network Modules

We used the MCODE plug-in in the Cytoscape software¹⁷ to screen the modules concealed in the carotid atherosclerosis occurrence-related network with the following criteria: Max. depth = 100, K-Core = 7, mode score cutoff = 0.2, and degree cutoff = 2. Likewise, the functional specification of the identified module was determined with the clusterProfiler package as mentioned above. A P value (adjusted by Benjamini & Hochberg method) <0.05 was considered statistically significant.

Statistical Analysis

All analyses were performed using R (version 4.0.2, http://www.r-project.org/). We analyzed the expression levels of genes in the atheroma plaque and control samples using unpaired t-tests provided by limma package. Unless otherwise stated, we considered P < 0.05 to be statistically significant.

Results

Identification of DEGs Associated with the Progression of Carotid Atherosclerosis

A total of 1328 DEGs, including 749 up-regulated DEGs and 579 down-regulated DEGs, were identified in atheroma plaque compared with control samples (Figure 2A). In addition, we found that the expression patterns of these DEGs may distinguish atheroma plaque from control samples (Figure 2B)

Figure 2 The differentially expressed genes (DEGs) for carotid atherosclerosis and control samples. (A) A volcano plot of the for carotid atherosclerosis and control samples. Red indicated upregulated genes and blue indicated downregulated genes. The gray area showed the gene expression below the threshold criteria. (B) The heatmap of the DEGs. Abbreviation: DEGs, differentially expressed genes.

KEGG and GO Enrichment Analyses of DEGs

To analyze the biological classification of DEGs, functional and pathway enrichment analyses were performed. As shown in Figure 3A, changes in biological processes (BPs) of DEGs were significantly enriched in the physical activity of immune cells (neutrophil, leukocyte, T cell, lymphocyte and phagocyte). GO analysis result showed that cellar component (CC) and molecular function (MF) were significantly enriched in “extracellular matrix”, “plasma membrane protein complex”, “collagen-containing extracellular matrix”, “cell adhesion molecule binding” “actin binding” and “amide binding” (Figure 3B and C). The KEGG enrichment analysis result showed that some pathways were related to carotid atherosclerosis, such as PI3K-Akt signaling pathway,¹⁸ Chemokine signaling pathway,¹⁹ Cell adhesion molecules,²⁰ MAPK signaling pathway,²¹ Calcium signaling pathway,²² NOD-like receptor signaling pathway²³ and Fluid shear stress and atherosclerosis²⁴ (Figure 3D). Furthermore, we found that the change of BPs was significantly “leukocyte migration”, “leukocyte proliferation”, “myeloid leukocyte migration”, “phagocytosis”, “regulation of lymphocyte activation”, “response to molecule of bacterial origin” and “T cell activation” in atheroma plaque compared with control samples based on GSEA result (Figure 3E). And the KEGG pathways, “B cell receptor signaling pathway”, “chemokine signaling pathway”, “cytokine-cytokine receptor interaction”, “leukocyte trans endothelial migration” “NOD like receptor signaling pathway”, were significantly enriched in atheroma plaque compared with control samples based on GSEA result (Figure 3F).

Figure 3 The biological functions of differentially expressed genes (DEGs). (A) Cellular component, (B) biological process, (C) molecular function and (D) KEGG pathway analysis. The gene set enrichment analysis (GSEA) for (E) GO and (F) KEGG for carotid atherosclerosis and control samples. Abbreviations: GO, gene ontology; KEGG, Kyoto encyclopedia genes and genome.

PPI Network Construction

To obtain the interactions between the 1328 DEGs in the atheroma plaque compared with control samples, a PPI network was constructed using the STRING database and visualized by the Cytoscape software. As shown in Figure 4, the network included 147 nodes and 178 edges.

Figure 4 The PPI network. A graphic representation of the protein–protein network based on differentially expressed genes of carotid atherosclerosis and control samples. Red nodes indicated upregulated genes and light blue nodes indicated downregulated genes in the carotid atherosclerosis group.

Hub Genes Selection and Analysis

Among the lists of the top 50 genes selected, respectively, by the four ranking methods (Figure 5A–D), 16 genes overlapped were identified (Figure 5E). And 16 hub genes (C3AR1, CCR1, CCR2, CD33, CD53, CXCL10, CXCL8, CXCR4, CYBB, FCER1G, FPR2, ITGAL, ITGAM, ITGAX, ITGB2, and LILRB2) were upregulated in carotid atherosclerosis tissue compared with control samples (Figure 5F).

Figure 5 The selection and analysis of hub genes. Top 50 genes selected based the four ranking methods. (A) Degree, (B) DMNC, (C) MCC and (D) MNC. The darker the color of the node, the higher the score. DMNC, density of maximum neighborhood component; MCC, maximal clique centrality. MNC, maximum neighborhood component. (E) Venn plot. Venn plot of the overlapping genes for the top 50 genes selected based on the four ranking methods. Degree, DMNC, MNC, and MCC. (F) The expression levels for 16 hub genes in carotid atherosclerosis samples compared with control samples.

KEGG and GO Enrichment Analyses of Cluster Genes

A total of 5 clusters were obtained after MCODE (Molecular Complex Detection) algorithm. In addition, the enrichment analysis result showed that changes in GO function of cluster 1 were significantly enriched in “neutrophil activation”, “secretory granule membrane” and “G protein-coupled peptide receptor activity” (Figure 6A), while the changes in KEGG pathways were significantly enriched in “chemokine signaling pathway” and “cytokine-cytokine receptor interaction” (Figure 6B). The changes in GO function of cluster 2 were significantly enriched in “response to type I interferon”, “defense response to other organism”, “endocytic vesicle”, “plasma membrane protein complex”, “toll-like receptor activity”, “protein tyrosine kinase activity” and “cytokine receptor binding” (Figure 7A), while the changes in KEGG pathways of that were significantly enriched in “Tuberculosis” (Figure 7B). The changes in GO function of cluster 3 were significantly enriched in “neutrophil activation”, “azurophil granule lumen”, “primary lysosome”, “azurophil granule”, “vacuolar lumen”, “secretory granule lumen”, “vesicle lumen”, “receptor ligand activity”, “enzyme activator activity” and “receptor regulator activity” (Figure 8A), while the changes in KEGG pathways of that were significantly enriched in “lysosome” and “rRheumatoid arthritis” (Figure 8B). The changes in GO function of cluster 4 were significantly enriched in “regulation of leukocyte activation”, “cytoplasmic vesicle lumen”, “vesicle lumen”, “ubiquitin-protein transferase activity” and “ubiquitin-like protein transferase activity” (Figure 9A), while the changes in KEGG pathways of that were significantly enriched in “PI3K-Akt signaling pathway”, “human papillomavirus infection”, “endocytosis”, “focal adhesion” and “cell adhesion molecules” (Figure 9B). The changes in GO function of cluster 5 were significantly enriched in “muscle system process”, “actin cytoskeleton” and “actin binding” (Figure 10A), while the changes in KEGG pathways of that were significantly enriched in “vascular smooth muscle contraction” (Figure 10B). These results indicated that the formation of carotid atherosclerosis is the result of a coordination disorder of multiple gene modules and various biological pathways.

Figure 6 The GO and KEGG enrichment analysis for the genes in cluster 1 obtained based the MCODE method. (A) The biological processes and (B) KEGG pathways for genes in cluster 1.

Figure 7 The GO and KEGG enrichment analysis for the genes in cluster 2 obtained based the MCODE method. (A) The biological processes and (B) KEGG pathways for genes in cluster 2.

Figure 8 The GO and KEGG enrichment analysis for the genes in cluster 3 obtained based the MCODE method. (A) The biological processes and (B) KEGG pathways for genes in cluster 3.

Figure 9 The GO and KEGG enrichment analysis for the genes in cluster 4 obtained based the MCODE method. (A) The biological processes and (B) KEGG pathways for genes in cluster 4.

Figure 10 The GO and KEGG enrichment analysis for the genes in cluster 5 obtained based the MCODE method. (A) The biological processes and (B) KEGG pathways for genes in cluster 5.

Discussion

A study showed that about 57.79 million people of 30~70-year-old in the world have atheroma plaque in 2020, an increase of 59.13% from 2000.²⁵ However, the efficiency of treatment was low. To find potentially effective and therapeutic targets, there is an urgent requirement to explore the genes that lead to development of carotid atherosclerosis.

The present study utilized a relatively large sample dataset obtained from GEO. These data were analyzed to identify DEGs between carotid atherosclerosis and control samples. KEGG pathway analysis revealed with the 1328 DEGs mainly involved with PI3K-Akt signaling pathway, Chemokine signaling pathway, Cell adhesion molecules, MAPK signaling pathway, Calcium signaling pathway, NOD-like receptor signaling pathway and Fluid shear stress and atherosclerosis. It is well known that the inhibition of PI3K-Akt signaling pathway may reduce the vulnerability of atherosclerotic plaques.¹⁸ And chemokine CXCR4 may limit atherosclerosis by maintaining arterial integrity.¹⁹ Furthermore, geniposide may be against atherosclerosis by inhibiting formation of foam cell and lowering reverse lipid transport via MAPK signaling pathway.²¹ Circulating cell adhesion molecules are correlated with ultrasound-based assessment of carotid atherosclerosis.²⁰ GO enrichment analysis revealed that DEGs were mainly associated with physical activity of immune cells (neutrophil, leukocyte, T cell, lymphocyte and phagocyte). The present results provide bioinformatics evidence for further research.

The 16 overlapping genes among the top 50 genes in the PPI network found using four ranking methods were selected. All 16 genes were upregulated in the carotid atherosclerosis. The result reveals that these genes may serve an important role in the progression of carotid atherosclerosis. And prior the present study, few studies have addressed the gaps in the molecular mechanisms that lead to carotid atherosclerosis development. The genes (C3AR1, CCR1, CCR2, CXCL10, CXCL8, CXCR4) enriched for chemokine-related genes. The genes (ITGAL, ITGAM, ITGAX, ITGB2) enriched for cell adhesion-related molecules. The genes (CD33, CD53, CYBB, FCER1G, FPR2, LILRB2) enriched for immunity-related genes. In carotid plaques, C3aR1²⁶ and FCER1G²⁷ is higher compared to control arteries. CCR1 alters the immuno-inflammatory response in atherosclerosis.²⁸ Anti-CXCL10 treatment in atherosclerosis susceptible mice results in a change into a more stable lesion phenotype.²⁹ The progression of carotid atherosclerosis is related to CD53.³⁰ In the present study, the aim was to focus on DEGs between atheroma plaque and control samples. However, due to lack of experimental validation, it is not clear whether these genes are causal or merely markers. Furthermore, we established of the carotid atherosclerosis development-related network modules based on the MCODE method. A total of 5 clusters were obtained. Furthermore, the enrichment analysis results of each cluster genes showed that biological functions obtained were mostly related to immunity, chemokines and cell adhesion molecules. But there are some differences.

A total of 5 clusters were obtained after MCODE algorithm. And the enrichment analysis result showed that the biological function were associated with carotid atherosclerosis. Interestingly, NF-kb, a crucial transcriptional factor, which controls the transcription of many genes with an established role in atherosclerosis, such as cell proliferation, regulators of apoptosis, acute phase proteins, adhesion molecules, chemokines and cytokines.³¹ Toll-like receptors (TLRs), key players in innate immunity, are upregulated in atherosclerotic lesions.³² Atherosclerosis is an inflammatory disease associated with the activation of innate immune TLRs and nucleotide-binding oligomerization domain-containing protein (NOD)-like receptor pathways.³³ In addition, a previous study showed that a large amount of iron is deposited in the lysosomes of foam cells in early atherosclerotic lesions.³⁴ In summary, we thought the 5 clusters were closely related to carotid atherosclerosis.

There are a few important limitations to the present study. Firstly, the study is only based on the bioinformatics analysis, lack of experimental validation. Secondly, the sample size in this study is not very large. We hope that we can expand the sample size included in the further. Finally, our study is based on a retrospective dataset. We hope that we can verify the findings of this study in a prospective cohort.

Conclusion

In conclusion, the present study identified key DEGs in atheroma plaque compared with control samples. The key genes involved in the development of carotid atherosclerosis may provide valuable therapeutic targets for carotid atherosclerosis.