1.
Introduction
Atherosclerosis (AS) is a progressive, chronic, lipid-driven inflammatory disease of the arterial wall and an underlying cause of most cardiovascular diseases [1]. Many factors contribute to the accumulation of atheromatous plaques in the arterial wall, such as hypercholesterolemia, hypertension, diabetes mellitus, obesity, and smoking [2]. Despite the development of preventive strategies targeting classical risk factors, AS remains a major cause of morbidity and mortality worldwide [3]. At present, the lack of timely detection and effective control has led to ineffective control of the occurrence and development of atherosclerotic stenosis and vulnerable plaques [4]. Therefore, it is imperative to explore the accurate molecular markers and potential pathogenesis to contribute to the diagnosis and treatment of atherosclerotic diseases.
AS is a chronic inflammatory disease that involves infiltration of immune cells such as neutrophils, macrophages, T cells, and B cells, into the inner layer of vessel walls [5]. The role of immune cells in AS has been widely studied. For instance, a previous study reported T cell recruitment at an early stage in the development of atherosclerotic lesions [6]. In the presence of low-density lipoprotein, chemoattractants and adhesion molecules are expressed in the intima, leading to T cell infiltration into the plaques [7]. The function of B cell subsets in atherogenesis is controversial, as it has both atherogenic and atheroprotective properties [8]. Although numerous studies have reported the underlying pathogenesis of immune cell infiltration in AS, only a few studies have analyzed the potential genes associated with immune cell infiltration in AS [9].
In this study, we downloaded the AS-associated gene expression profile dataset (GSE57691) from the NCBI Gene Expression Omnibus (GEO) database and analyzed the proportion of immune cells and immune infiltration-related genes via bioinformatic analysis. A competing endogenous RNA (ceRNA) network was subsequently constructed to explore the potential regulatory mechanisms of immune-related genes in AS. Our results provide evidence of the involvement of the genes related to immune cell infiltration in the occurrence and development of AS.
2.
Materials and method
2.1. Data sources and preprocessing
On June 16, 2020, the gene expression profile data (GSE57691 [10]) of human arterial tissue samples were downloaded from the NCBI GEO [11] (https://www.ncbi.nlm.nih.gov/) database. This dataset contained 68 samples, and only 19 samples related to AS were selected for analysis, including nine samples of arterial occlusion and 10 samples from healthy controls. The clinical information of these 19 AS-related samples is shown in Supplementary Table 1. The test platform was an Illumina HumanHT-12 V4.0 expression beadchip.
2.2. Screening of significantly differentially expressed (DE) RNAs
The detailed platform annotation information (probe corresponding RNA type and symbol) of the Illumina HumanHT-12 V4.0 expression beadchip was downloaded from the ENSEMBL database (https://asia.ensembl.org/index.html). Based on the long-noncoding RNA (lncRNA), microRNA (miRNA), and messenger RNA (mRNA) annotation information included in the HUGO Gene Nomenclature Committee (HGNC) [12] database (http://www.genenames.org/), the lncRNAs, miRNAs, and mRNAs in the expression profile dataset were annotated.
The limma version 3.34.7 [13] in R3.4.1 language was subsequently used to screen significantly DERNAs (lncRNA, miRNA, and mRNA) between the AS and the control comparison groups, and false discovery rate (FDR) < 0.05 and |log2 fold change (FC)| > 1 were selected as the thresholds. According to the expression levels of these identified DERNAs in samples, bidirectional hierarchical clustering based on the Pearson correlation algorithm was conducted using R3.4.1 pheatmap version 1.0.8 (https://cran.r-project.org/web/packages/pheatmap/index.html) [14].
2.3. Quantitative evaluation of the immune cell proportion in the dataset sample
Each gene expression profile dataset was uploaded to estimate the proportion of immune and cancer cells (EPIC, https://gfellerlab.shinyapps.io/EPIC_1-1/) [15], and the proportion of immune-infiltrating cells in each sample was assessed based on genome-wide expression levels. EPIC can analyze the infiltration rate of the following seven immune cells based on the expression data: B cells, CD4+ T cells, CD8+ T cells, endothelial cells, macrophages, natural killer cells, and other immune cells.
2.4. Identification of genes associated with immune infiltration
Using the cor function in R3.4.1, the Pearson's correlation coefficient between the significant DEmRNAs and the immune-infiltrating cell types in the sample was calculated, and significant DEmRNAs closely associated with immune infiltration were identified. In addition, the biological processes related to immunity were searched from the AmiGO 2 database (http://amigo.geneontology.org/amigo) with "immune" as the keyword, and the genes associated with immune-related biological processes were downloaded. Meanwhile, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp/) was used to download all pathways associated with "immune." Next, the intersection sections between DEmRNAs and immune function-related genes were selected for further analyses, including Gene Ontology (GO) and KEGG enrichment analyses (criteria for significance screening: p value < 0.05) using DAVID 6.8 (https://david.ncifcrf.gov/) [16,17].
2.5. ceRNA network construction and functional analyses
The connection relationships between DElncRNA and DEmiRNAs were searched from the DIANA-LncBasev2 database [18] (http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php?r=lncbasev2%2Findex-experimental), and only the connection pairs with opposite DE direction were retained to construct DElncRNA-DEmiRNA connection network. Thereafter, the starBase Version 2.0 database [19] (http://starbase.sysu.edu.cn/) that integrates the target gene prediction information from TargetScan, PicTar, RNA22, PITA, and miRanda databases was used to search for the target genes of DEmiRNAs. We selected the regulatory relationships included in at least one of the databases as DEmiRNA-target gene pairs. Then, the DEmRNAs were mapped to the target genes, and the miRNA and mRNA pairs that were negatively correlated in expression levels were retained to build a DEmiRNA-DEmRNA regulatory relationship.
Finally, DElncRNA-DEmiRNA and DEmiRNA-DEmRNA connection networks were synthesized, and a ceRNA network comprising DElncRNAs, DEmiRNAs, and DEmRNAs was proposed and visualized using Cytoscape version 3.6.1 (https://cytoscape.org/) [20]. DAVID 6.8 was used for GO biological process and KEGG pathway enrichment analyses of these genes in the ceRNA regulatory network. A P value less than 0.05 was selected as the threshold for screening significant enrichment.
3.
Results
3.1. DERNA identification
In the dataset, 17861 mRNAs, 902 lncRNAs, and 46 miRNAs were annotated. For the AS versus control comparison groups in the dataset, a total of 1749 DERNAs meeting the thresholds were screened, including 1673 mRNAs (1259 downregulated and 414 upregulated), 63 lncRNAs (14 downregulated and 49 upregulated), and 13 miRNAs (8 downregulated and 5 upregulated). The volcano plot is shown in Figure 1A, and the distribution of DElncRNAs, DEmiRNAs, and DEmRNAs is displayed in Figure 1B. The bidirectional hierarchical clustering heatmaps based on the expression level of DERNAs are shown in Figure 1C, which infers that the three types of RNAs could effectively differentiate the AS samples from the healthy control samples.
3.2. Quantitative evaluation of the proportion of immune cells
EPIC was used to quantitatively evaluate the types of immune-infiltrating cells in each sample, as shown in Figure 2A. We used the intergroup T test to compare each type of immune-infiltrating cells in the AS and control groups. As shown in Figure 2B, the proportions of B cells, CD4+ T cells, and CD8+ T cells were significantly different between the two groups.
3.3. Screening of genes associated with immune cell infiltrationk
The correlation coefficients between the three types of immune-infiltrating cells and DEmRNAs were calculated, and only those genes with absolute values of correlation coefficients greater than 0.2 were retained. Finally, 1614 genes were identified. In total, 3020 unique genes related to immune function and 817 unique genes related to immune pathways were downloaded from the database. As shown in Figure 3, 582 intersection genes were observed, which were further compared with 1614 genes above to obtain 341 intersection genes, including HBB, FCN1, IL1B, CXCL8, RPS27A, CCN3, CTSZ, and SERPINA3. The overlapping genes were subjected to GO biological process and KEGG signal pathway enrichment analyses, and a total of 70 significantly related GO biological processes (such as immune response, positive regulation of molecular function, and response to wounding) and 15 KEGG pathways (such as chemokine signaling pathway, leukocyte transendothelial migration, and proteasome) were observed (Figure 4).
3.4. ceRNA network construction and function analysis
After searching the DIANA-LncBasev2 database, 89 DElncRNA-DEmiRNA pairs were obtained. By mapping the 314 immune-related genes to the regulatory target genes of the DEmiRNAs, 409 DEmiRNA-DEmRNA pairs were identified. By integrating DElncRNA-DEmiRNA and DEmiRNA-DEmRNA interaction pairs, a ceRNA network was constructed that included 33 lncRNAs (SNHG10, MBNL1-AS1, PKN2-AS1, LINC01165, MIAT, etc.), 11 miRNAs (miR-6764, miR-6514, miR-1226, miR-7-1, etc.), and 216 mRNAs (e.g., FCN1, SERPINA3, and IL1B) (Figure 5). GO and KEGG enrichment analyses showed that these mRNAs in the ceRNA network were significantly enriched in immune response, positive regulation of molecular function, positive regulation of catalytic activity, response to wounding, and chemotaxis, as well as pathways of chemokine signaling pathway, leukocyte transendothelial migration, regulation of actin cytoskeleton, cytokine-cytokine receptor interaction, and T cell receptor signaling pathway (Figure 6).
4.
Discussion
AS is a chronic inflammatory disorder of the arteries, characterized by subendothelial accumulation of immune cells, extracellular matrix deposition, and cholesterol [21]. In this study, we focused on immune cell infiltration in AS and identified 341 immune-related genes. Based on these genes and the identified DEmiRNAs and DElncRNAs, a ceRNA network comprising 33 lncRNAs, 11 miRNAs, and 216 mRNAs was constructed, and the mRNAs in this network were found to be significantly enriched in some immune-associated functions and pathways.
In comparison with the healthy condition, AS can contribute to the change in cell types in the arterial tissue, which expresses different species of RNA, and T and B cells in the media of the vascular wall play important roles during the development of atherosclerotic processes [22]. Our results showed that B cells, CD4+ T cells, and CD8+ T cells were significantly differently expressed between the two groups, which further indicated their important role in AS. During AS, T cells can recognize self-antigens and induce humoral immunity driven by B cells [23], leading to the destruction of active immune privileges in the media. CD4+ T cells are the main T cell type at all stages of lesion development. It has been reported that local immune responses dominated by CD4+ T cells occur at all stages of AS [24]. Additionally, abundant CD8+ T cell infiltrates have been identified in atherosclerotic lesions of the carotid artery [25]. However, in contrast to CD4+ T cells, no AS-related antigens that may activate CD8+ T cells have been found [26]. Based on our results, we speculate that AS progression may be closely associated with changes in B cells, CD4+ T cells, and CD8+ T cells.
To investigate the underlying molecular mechanisms of AS, immune cell infiltration-related genes were analyzed. A total of 341 immune infiltration-related genes were identified. By combining the obtained DEmiRNAs and DElncRNAs, a ceRNA network was established to investigate the potential regulatory mechanism in AS. FCN1 (ficolin 1) had the highest |log2 FC| among the immune-related upregulated genes in the network. FCN1 is involved in the complement lectin pathway and is elevated in patients with Takayasu arteritis [27]. Abnormal expression of FCN1 has been suggested as a potential target for coronary artery diseases [28]. To our knowledge, the role of FCN1 in AS has not been well investigated. In this study, FCN1 was enriched in immune response-associated GO function, which further demonstrated the important function of the immune system in AS. Further, FCN1 is involved in the SNHG10-miR-6764-FCN1 axis. We speculate that SNHG10 may function as a ceRNA to promote AS progression by upregulating FCN1.
IL1B (interleukin 1 beta) also had a higher log2 FC among the immune-related genes in the ceRNA network, and was enriched in the cytokine-cytokine receptor interaction pathway. IL1B, one of the first recognized cytokines, is a critical driving factor of local and systemic immune responses in atherosclerotic cardiovascular disease [29]. A study suggested that IL1B can promote a pro-atherogenic phenotypic shift in vascular smooth muscle cells by downregulating the expression of smooth muscle markers and inducing the expression of pro-inflammatory cytokines [30]. In the constructed ceRNA network, IL1B was found to be regulated by miR-7-1 in the ceRNA relationship of PKN2-AS1-miR-7-1-IL1B. We speculate that PKN2-AS1 may serve as a ceRNA to facilitate AS by regulating IL1B expression.
SERPINA3 (serpin family A member 3) had the highest |log2 FC| among the immune-related downregulated genes in the network, and was regulated by miR-6514. It was involved in the ceRNA axis of SERPINA3-miR-6514-MBNL1-AS1. SERPINA3 belongs to the serpin superfamily and plays an important role in vessel wall biology [31]. The plasma SERPINA3 level is elevated in acute myocardial infarction as compared with healthy controls [32]. Importantly, SERPINA3 is expressed in the medial smooth muscle cells and endothelial cells of human atherosclerotic lesions [33]. Further, it can regulate immune cells via elastase and cathepsin G regulation [31,34], indicating its role in immunity.
However, this study has some limitations. The number of samples is small, necessitating further studies with larger sample sizes to verify this conclusion. In addition, the immune cell-specific or non-specific biomarkers associated with AS should be further explored, and experimental studies on human and mouse models should be carried out to confirm our findings.
5.
Conclusion
In conclusion, our study identified 1673 DEmRNAs, 63 DElncRNAs, and 13 DEmiRNAs between AS patients and healthy controls, and 341 immune-related genes. B cells, CD4+ T cells, and CD8+ T cells were found to be closely related to AS progression. A ceRNA regulatory network comprising 33 lncRNAs, 11 miRNAs, and 216 mRNAs was proposed, and immune-related genes such as FCN1, IL1B, and SERPINA3 were thought to serve as potential risk markers of atherosclerotic immune cell infiltration. These genes may promote AS progression via ceRNA regulation. Our findings will help improve the understanding of the occurrence and development of AS and explain novel underlying pathological molecular mechanisms for AS.
Acknowledgments
This study was supported by the National Natural Science Foundation of China Youth Science Foundation Project (No. 81700304), the Scientific Research Fund Project of Key Laboratory of the Second Hospital of Tianjin Medical University (No. 2019ZDSYS11), the Foundation of Tianjin Union Medical Center (No. 2017YJ006), and Peking Union Medical Foundation - Ruiyi Emergency Medical Research Fund (No. R2018006), the National Natural Science Foundation of China (No. 82072222), and the Fundamental Research Funds for the Central Universities, China (No.3332019127). These funders have no role in this study.
Conflict of interest
The authors declare that they have no conflict of interest.
DownLoad:
