Novel immune cell infiltration-related biomarkers in atherosclerosis diagnosis

Data source

In the current study, gene expression data from 32 AS plaque samples at stage IV and/or V lesions including core and shoulders of the plaque and 32 distant macroscopically intact control samples in GSE43292 (Ayari & Bricca, 2013) were used as the testing set to find diagnostic AS biomarkers. The demographic data of GSE43292 cohort have been reported by Ayari & Bricca (2013) in a previous literature. In addition, 29 atherosclerotic carotid artery samples and 12 healthy artery samples in GSE100927 (Steenman et al., 2018) were used as an external set to validate the expression patterns and diagnostic value of biomarkers identified in GSE43292. GSE43292 and GSE100927 were sourced from GEO database (http://www.ncbi.nlm.nih.gov/geo). The workflow of the current study was presented in Fig. 1.

Exploration of differentially expressed genes (DEGs) in GSE43292

Limma program in R was applied to mine DEGs using |logFC (fold change) | ≥ 1 and adjusted p-value < 0.05. The “ggplot2” package of R was used to create volcano plot. The heatmap was produced using “pheatmap” in R.

Identification of differentially infiltrated immune cells (DIICs) between AS and control

The proportions of IICs in AS and control samples was evaluated using the ssGSEA algorithm. This process was completed by using “GSVA” R package. The immune cells exhibiting significant differences between AS and control were examined using the Wilcoxon method with Benjamin & Hochberg adjusted p-value < 0.05.

Weighted co-expression network analysis

WGCNA was performed (Horvath & Dong, 2008) on all samples in the training set in order to screen the gene modules most associated with DIICs. To remove outliers, the hierarchical clustering trees for all samples was constructed, followed by the selection of optimal β value to build the scale-free network. Next, Pearson’s correlations between DIICs and gene modules were determined and presented in the heatmap. Finally, the most positive module and the most negative module correlated with DIICs were selected as key modules, and key moduler genes were used for the following analysis.

Mining and analysis of DEGs related to DIICs

To get DEGs related to DIICs, we overlapped DEGs with key modular genes identified in WGCNA and presented by a Venn diagram. Furthermore, KEGG pathway and GO analysis, comprising molecular function (MF), cellular composition (CC) and biological process (BP) were performed by “clusterProfiler” in R (Yu et al., 2012).

Identification of hub DIIC-related DEGs in AS

To identify hub DIIC-related DEGs in AS, a protein-protein interaction (PPI) network was first developed by uploading DIIC-related DEGs into STRING database (http://string-db.org) (Szklarczyk et al., 2015). Then Cytoscape software was applied for identifying core network by MCODE plug-in (Shannon et al., 2003). Genes in the core network were defined as hub DIIC-related DEGs in AS. Next, receiver operating curves (ROC) were plotted to assess the role of hub genes in AS diagnosis. If a hub gene has an area under ROC (AUC) > 0.7, then it was considered as a potential marker in AS diagnosis. In addition, their expressions and diagnostic potential of hub genes were validated in GSE100927. Then, the hub genes with consistent expression patterns and diagnostic potential were taken into the next analyses.

Characteristics and regulatory network of hub genes

To investigate the characteristics of hub genes and their relationships, we (1) analyzed their functional similarity by R package “GOSemsim”, (2) performed Spearman correlation analysis to determine whether their expressions were correlated, and (3) calculated angiogenesis and immune checkpoint scores using ssGSEA algorithm, followed by the calculation of correlations between hub gene expressions and angiogenesis/immune checkpoint scores. Furthermore, hub genes were imported into the miRNet database (https://www.mirnet.ca/) to search for miRNAs and transcription factors (TFs) that target and regulate the expressions of hub genes. Finally, the regulatory networks of miRNAs-hub genes and TFs-hub genes were developed using Cytoscape software.

Cell culture and RT-qPCR

To validate hub genes’ expressions, we purchased THP-1 cell line from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China, CAS No: KG224). Cell culture were performed in RPMI-1640 medium (Sigma-Aldrich, St. Louis, MO, USA) containing 10% fetal bovine serum (FBS) in T25 flasks. At 37 °C and 5% CO₂ THP-1 monocytes were induced to differentiate into macrophage by adding PMA at 100 ng/ml to the media for 48 h. Subsequently, the macrophages were incubated with 80 ug/ml ox-LDL for 24 h to transform into foam cells. Total RNA were extracted from foam cells and macrophages, respectively, using TRIzol Reagent (Thermal Fisher Scientific, Waltham, MA, USA). The purity and concentration of extracted RNA were detected for the following cDNA synthesis by SureScript-First-strand-cDNA-synthesis-kit (Servicebio, Guangzhou, China). Next, cDNA was applied for qPCR under 40 cycles at 95 °C for 60 s, 95 °C for 20 s, 55 °C for 20 s and 72 °C for 30 s using qPCR reagent from Servicebio (Wuhan, China). Gene expressions were analyzed by 2^−ΔΔCt method. Primers of HMOX1, CHI3L1 and MMP9 were listed in Table 1.

102 DEGs and 27 DIICs were identified between AS and control

We discovered that 49 and 53 genes were significantly elevated and reduced, respectively, in AS than in controls in GSE43292 (Fig. 2A). The expressions of top 40 DEGs were presented in the heatmap (Fig. 2B). Meanwhile, ssGSEA algorithm was applied to estimate the infiltrations of 28 ICs in controls and AS as shown in Fig. 2C. After the Wilcoxon test, we observed that the average proportions of 28 ICs were much higher in AS group, and there were significant differences of ICs between AS and control except for type 2 T helper cells (Fig. 2D).

99 DIIC-related DEGs were screened by WGCNA

Next, we performed WGCNA to obtain DIIC-correlated modules. The sample clustering tree identified no outlier in the training set (Fig. 3A), followed by the generation of the sample dendrogram and trait heatmap (Fig. 3B). Afterwards, 10 is selected as the best soft threshold to develop the scale-free network (Fig. 3C). Finally, eight gene modules, including yellow, blue, gold, black, pink, blueviolet, turquoise and brown, were obtained (Fig. 3D). The correlations between DIICs and modules were calculated and presented in the heatmap (Fig. 3E). Accordingly, we selected most negatively (brown) and positively (pink) correlated modules with DIICs. Finally, a total of 4,568 genes in the pink and brown modules were screened as DIIC-related genes in AS. Thereafter, by overlapping 102 DEGs with 4,568 DIIC-related genes, we identified 99 DIIC-related DEGs (Fig. 3F).

By GO and KEGG pathway analyses, we identified 24 significantly enriched KEGG pathways, the top 10 of which were tryptophan metabolism, cAMP signaling pathway, long-term potentiation, PPAR signaling pathway, diabetic cardiomyopathy, complement and coagulation cascades, circadian entrainment, ECM-receptor interaction, Malaria and long-term depression (Figs. 4A and 4B). As for GO analysis, 52 BP, 11 MF and 21 CC were significantly enriched. As shown in Figs. 4C and 4D, the top 10 GO terms, such as tryptophan metabolic process, extracellular matrix disassembly, aromatic amino acid family metabolic process, and the corresponding DIIC-related DEGs involved in those GO terms, such as FABP4, MMP9, CXCL10 were presented.

10 hub genes could distinguish AS from control samples

The interactions among DIIC-related DEGs were explored in the STRING database (Fig. 5A) and were visualized by Cytoscape software (Fig. 5B). Then we screened core network from PPI network by MCODE plug-in (Fig. 5C). A total of 11 DIIC-related DEGs, including CASQ2, NEXN, MMP12, CXCL10, PLN, HMOX1, SELE, RYR2, CHI3L1, MMP9, ACTC1 were extracted from the core network (Fig. 5D). Among them, the expressions of MMP9, MMP12, HMOX1, CHI3L1, SELE and CXCL10 were upregulated, while the expressions of CASQ2, NEXN, RYR2, ACTC1 and PLN were down-regulated in AS samples compared to controls from training set (Fig. 6A). The ROC curves in the training set showed that those 11 genes had the potential to classify AS from controls (AUC > 0.7, Fig. 6B). Moreover, in the external dataset, the expression patterns and diagnostic value of CASQ2, NEXN, MMP12, CXCL10, PLN, HMOX1, RYR2, CHI3L1, MMP9, ACTC1 were validated (Figs. 6C and 6D). Thus, CASQ2, NEXN, MMP12, CXCL10, PLN, HMOX1, RYR2, CHI3L1, MMP9, ACTC1 were identified as hub DIIC-related DEGs and may act as AS diagnostic biomarkers.

Hub genes had close relationship with angiogenesis and immune checkpoints

Next, we explored the relationship among the hub genes. Correlation analysis showed that the expressions of them were highly correlated (Fig. 7A). Among them, MMP9 and CHI3L1 were most positively correlated (Fig. 7B), and NEXN and HMOX1 were most negatively correlated (Fig. 7C). Interestingly, we also found that there was a medium to high degree of similarity in their function (Fig. 7D). Moreover, in consideration the role of angiogenesis (Moulton, 2006; Perrotta et al., 2019) and immune checkpoints (Vuong et al., 2022; Alyagoob, Lahmann & Joner, 2020) in AS development, progression and treatment, we calculated the angiogenesis and immune checkpoints scores based on corresponding gene sets by ssGSEA algorithm. Further analyses showed that 10 hub genes were significantly correlated with angiogenesis score (Fig. 8A) and with immune checkpoints score (Fig. 8B).

Hub genes were regulated by multiple miRNAs and TFs

Next, we investigated the miRNAs and TFs targeting and regulating the expressions of 10 hub genes. A total of 313 miRNA-hub gene pairs were searched to establish the miRNA-hub gene network, which includes 239 miRNAs and 10 hub genes (Fig. 8C). In the network, we found some hub genes were regulated by common miRNAs, for example that MMP9 and RYR2 were both regulated by mir-9-3p and mir-29b-3p. Meanwhile, 53 TF-hub gene pairs were obtained in miRNet database and a TF-hub gene network composed of six hub genes and 47 TFs was constructed (Fig. 8D). Also, we observed that two hub genes were regulated by the same TF, such as CHI3L1-STAT3-HMOX1, CHI3L1-SP1-MMP9, and CXCL10-RELA-MMP12.

RT-qPCR validation

Foam cells derived from macrophages may indicate the initial stages of AS (Chistiakov, Bobryshev & Orekhov, 2016; Maguire, Pearce & Xiao, 2019). Thus, in the current study, the foam cell model established from THP-1 monocytes was used, which was widely used in studying AS and AS-related disease (Yin et al., 2019; Huwait, Al-Saedi & Mirza, 2022; Mehta & Dhawan, 2020). After treatment with PMA, THP-1 cells were differentiated into macrophages and were defined as control group. Then THP-1 derived macrophages were treated with ox-LDL to form foam cells, which mimic the early stages of AS. Then we examined the expressions of HMOX1, CHI3L1 and MMP9 by RT-qPCR. The results showed that the expressions of HMOX1, CHI3L1 and MMP9 were significantly up-regulated in foam cells, which were consistent with the sequencing results (Fig. 9).