Identification of key genes and pathways in abdominal aortic aneurysm by integrated bioinformatics analysis.

OBJECTIVES: To identify key genes associated with abdominal aortic aneurysm (AAA) by integrating a microarray profile and a single-cell RNA-seq dataset.
METHODS: The microarray profile of GSE7084 and the single-cell RNA-seq dataset were obtained from the Gene Express Omnibus database. Differentially expressed genes (DEGs) were chosen using the R package and annotated by Gene Ontology and Kyoto Encyclopedia of Genes and Genomics analysis. The hub genes were identified based on their degrees of interaction in the protein-protein interaction (PPI) network. Expression of hub genes was determined using single-cell RNA-seq analysis.
RESULTS: In total, 507 upregulated and 842 downregulated DEGs were identified and associated with AAA. The upregulated DEGs were enriched into 9 biological processes and 10 biological pathways, which were closely involved in the pathogenesis and progression of AAA. Based on the PPI network, we focused on six hub genes, four of which were novel target genes compared with the known aneurysm gene database. Using single-cell RNA-seq analysis, we explored the four genes expressed in vascular cells of AAA: CANX, CD44, DAXX, and STAT1.
CONCLUSIONS: We identified key genes that may provide insight into the mechanism of AAA pathogenesis and progression and that have potential to be therapeutic targets.

Background

Abdominal aortic aneurysm (AAA) is a degenerative change in the abdominal aorta medial layer caused by various etiologies; the layer dilates to form an aneurysm under blood pressure, leading to permanent dilatation of the arterial wall to more than 150% of the normal vessel diameter. AAA is the tenth leading cause of death in elderly men.[1] Epidemiological studies have shown that the incidence of AAA in men increases significantly over the age of 55 years, and the mortality rate is as high as 5.9% in the 80 to 85 year age group.[2] Surgical repair, to date, is the main therapeutic regimen, but clinical care is hampered by the lack of etiological treatment (treatment specific to the etiology of an individual case). Therefore, it is of great clinical significance to find new molecular targets involved in initiation and progression of AAA.

AAA is characterized by chronic inflammation, vascular smooth muscle cell apoptosis, and extracellular matrix remodeling and degradation.[3] Inflammation is most associated with the progression of AAA.[4] Macrophages, located in the media and adventitia, promote the formation of AAA by producing matrix metalloproteinases (MMPs), reactive oxygen species, and inflammatory factors.[5] The adventitial mesenchymal cells, including fibroblasts, myofibroblasts, and “synthetic” vascular smooth muscle cells secrete type I collagen and transforming growth factor (TGF)-β, MMPs involved in extracellular matrix (ECM) remodeling.[6] Previous studies have identified individual lifestyle factors and physiological parameters associated with the incidence of AAA, such as smoking history and cholesterol levels.[7] Recently, a large genome-wide association study meta-analysis identified four new AAA risk loci: SMYD2, LINC00540, MMP9, and ERG.[8]

Thus, previous studies have focused on revealing the target genes associated with AAA. For example, Zhang et al. predicted AAA target genes by a novel protein–protein interaction (PPI) method.[9] However, the cell-specific expression of these hub genes remains to be elucidated. The present study aimed to identify the candidate genes involved in AAA using integrated bioinformatics techniques. We identified differentially expressed genes (DEGs) and then performed enrichment analysis to explore their involved functions and signaling pathways. We discovered hub genes by constructing a PPI network and explored the cell-specific expression of these hub genes. The DEGs may provide insight into the pathogenesis of AAA, which could help elucidate biomarkers and therapeutic targets of AAA.

Materials and methods

Microarray data

The National Center for Biotechnology Information Gene Expression Omnibus database (NCBI-GEO) is a public online database containing many microarray and RNA-seq datasets and from which we obtained two datasets. GSE7084[10,11] is an mRNA microarray profile based on the GPL570 platform (HG-U133_Plus_2; Affymetrix Human Genome U133 plus 2.0 Array), which consists of 10 control aortic tissues from autopsy and 9 human AAA tissues from patients undergoing surgical procedures. GSE118237[12] is a single-cell RNA sequencing profile based on the GPL21103 platform (Illumina HiSeq 4000; Illumina Inc., San Diego, CA, USA), consisting of 3 samples of angiotensin II (AngII)-induced murine AAA tissues. The microarray data of GSE7084 and GSE118237 can be obtained from the NCBI-GEO online database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi).

The current study protocol was approved by Huaian No. 1 People’s Hospital of Review Board in Huaian on 26 July 2019.

Identification of DEGs

Following sample quality control, DEG analysis was performed using the limma package in R software (http://www.bioconductor.org/packages/release/bioc/html/limma.html). The mRNAs with a P-value <0.05 and |log2 fold-change (FC)| > 1 were considered differentially expressed. Volcano plots, heatmaps, and principal component analysis (PCA) plots were also generated in R (www.r-project.org).

GO and KEGG pathway enrichment analysis

FunRich version 3 (http://www.funrich.org/) was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomics (KEGG) analysis, a common method to annotate genes and explore their biological attributes. Statistically significant biological process terms and biological pathway terms for DEGs were illustrated. Transcription factors (TFs) that might regulate DEGs were also predicted. Based on the cumulative hypergeometric distribution test, P < 0.05 was considered to indicate a statistically significant difference.

PPI network construction and visualization

PPI networks provide valuable information regarding cellular functions and signaling pathways. The online database Search Tool for the Retrieval of Interacting Genes/Proteins[13] (http://string-db.org/) was used to search for the interaction of proteins encoded by the identified DEGs. Cytoscape[14] (http://cytoscape.org/) was then used to visualize the PPI network, established based on five calculation methods (Degree, EPC, EcCentricity, MCC, and MNC). The intersecting genes calculated from these five algorithms encode core proteins with important biological regulatory functions.

Single-cell analysis of key genes involved in AAA

The single-cell RNA-seq dataset GSE118237 was downloaded from the NCBI-GEO online database. Data analysis was performed using the Loupe Cell Browser software (10x Genomics, Pleasanton, CA, USA) on Cloupe files. The clusters were displayed based on t-distributed stochastic neighbor embedding (t-SNE) projections of the cell transcriptome. Cell types were identified by their expression levels of cell-specific markers.

Results

Using the R limma package, we obtained 507 upregulated DEGs and 842 downregulated DEGs based on the cut-off criteria (P < 0.05 and |log2 FC| > 1). A volcano plot generated in R shows the DEGs with log2 FC scores and –log10 P-value (Figure 1a). The DEGs were clustered between AAA tissues and normal tissues, as shown in the heatmap (Figure 1b). The PCA plot also demonstrated that the DEGs could accurately distinguish AAA samples from non-AAA samples (Figure 1c).

Figure 1.

(a) Volcano plot of DEGs. Blue dots represented genes with |log2FC| > 1, P < 0.01; red dots represented the rest of the genes (with no significant expression change), (b) Heatmap of DEGs; red represents upregulation and blue represents downregulation in control (Ctl) and AAA tissues and (c) PCA plot of control (Ctl) and AAA tissues. The x- and y-axes denote the variance that each PC accounts for. DEG, differentially expressed gene; FC, fold change; AAA, abdominal aortic aneurysm; PCA, principal component analysis; PC, principal component.

GO and KEGG enrichment analysis

The functions and pathway enrichment of DEGs were analyzed using the Funrich software. As shown in Figure 2a, upregulated DEGs were enriched in signal transduction, cell communication, metabolism, energy pathways, transport, protein metabolism, immune response, apoptosis, and regulation of cell cycle in biological process. They were also involved in TRAIL signaling pathway, proteoglycan syndecan-mediated signaling events, glypican pathway, ErbB receptor signaling pathway, VEGF signaling network, thrombin/protease-activated receptor pathway, plasma membrane estrogen receptor signaling, IFN-gamma pathway, GMCSF mediated signaling events, and Alpha9 beta1 integrin signaling events (Figure 2b). The downregulated DEGs were specifically enriched in muscle contraction, cell growth and maintenance, and aldehyde metabolism in biological process, as well as sphingosine 1-phosphate pathway, class I PI3k signaling events, beta1 integrin cell surface interactions, and Arf6 trafficking events. Finally, we predicted the TF that might regulate our DEGs (Table 1).

Figure 2.

Bubble plots representing biological processes (a, c) and biological pathways (b, d) of upregulated (a, b) and downregulated (c, d) DEGs. Rich factor indicates the percentage of DEGs divided by all genes within certain GO term. DEG, differentially expressed gene; GO, Gene Ontology.

Table 1.

The predicted TF regulating DEGs.

DEGs	TF	Number	P-value
Upregulated	SP4	87	4.82E-02
	ASCL2	35	3.48E-02
	SOAT1	30	3.96E-02
	STAT1	30	3.96E-02
	IRF1	21	3.74E-02
	POU6F1	20	8.73E-03
	ISL2	16	1.55E-02
	BARX1	12	3.10E-02
Downregulated	SP1	244	1.78E-05
	SP4	163	8.45E-06
	KLF7	162	4.92E-05
	EGR1	121	1.70E-04
	HNF4A	105	6.46E-03
	NFYA	95	1.28E-03
	CTCF	88	2.60E-02
	ETS1	84	3.18E-02
	MYF5	83	7.47E-06
	POU2F1	83	2.27E-03
	RREB1	80	2.98E-03
	TCF3	79	4.60E-03
	MEF2A	78	8.25E-05
	NFIC	76	4.35E-03
	RFX1	73	4.01E-02
	PPARG	71	1.98E-02
	JUND	64	1.34E-02
	JUNB	64	1.34E-02
	FOSB	64	1.34E-02
	FOS	64	1.34E-02
	JUN	64	1.34E-02
	TFAP4	62	1.17E-05
	ASCL2	57	8.89E-03
	NKX6-1	56	3.87E-03
	LHX3	49	4.83E-03
	ARID3A	49	8.75E-03
	CREB1	47	1.51E-04
	NHLH1	44	1.71E-02
	HENMT1	44	1.71E-02
	PLAU	42	8.77E-03
	ATF1	42	8.77E-03
	RXRA	41	3.79E-02
	NR1H3	38	3.87E-02
	HSF1	37	1.87E-03
	HNF1A	37	3.12E-03
	NR6A1	36	1.54E-03
	RARA	35	5.87E-03
	RAB40B	35	5.87E-03
	GATA1	34	2.14E-03
	ZNF238	34	8.80E-03
	ZNF513	34	8.80E-03
	JDP2	33	7.78E-03
	ATF3	33	1.49E-02
	FOXA1	33	1.59E-02
	FEV	32	2.89E-05
	PBX1	29	3.90E-02
	PAX6	28	1.83E-03
	HOXB4	28	4.35E-02
	INSM1	27	5.25E-04
	MEIS2	27	1.68E-03
	LMO2	27	1.15E-02
	SRF	25	1.17E-03
	CACD	25	4.79E-03
	MAF	24	2.78E-03
	GFI1	23	3.43E-02
	ONECUT1	23	4.20E-02
	PAX7	22	1.11E-02
	DBX2	21	5.00E-03
	GSC	18	1.90E-02
	PRRX1	18	2.39E-02
	CDX2	17	7.16E-04
	EOMES	17	3.19E-03
	FOXC1	15	2.09E-02
	PITX3	14	1.61E-02
CEBPD	9	4.16E-03
MAFK	8	2.60E-02
HOXC10	6	2.47E-02

P-value <0.05 was statistically significant.

TF, transcription factor; DEG, differentially expressed gene.

All DEGs were submitted to STRING, a biological database and web resource for known and predicted PPI. The PPI network was constructed using Cytoscape software using the settings experiments, neighborhood, database, textmining, and coexpression. After removing isolated genes (those having no interaction with other genes), the upregulated DEGs were displayed with 34 nodes and 36 edges (Figure 3a). According to the degree of each gene, we identified six hub genes with a degree ≥4: AR, CANX, CD44, DAXX, STAT1, and TP53 (Figure 3b). Compared with those in the aneurysm gene database (AGD),[15] AR,[16] STAT1,[17,18] and TP53[19] were previously reported to be associated with AAA, whereas CANX, CD44, and DAXX were novel target genes deserving further investigation.

Figure 3.

(a) The PPI network of upregulated DEGs; red nodes represent hub genes with degree ≥4 and (b) Fold change of each hub gene in AAA compared with normal tissues. PPI, protein-protein interaction; DEG, differentially expressed gene; AAA, abdominal aortic aneurysm.

Single-cell analysis of hub genes involved in AAA

To evaluate the cell-specific expression of our hub genes, we analyzed a single-cell sequencing dataset consisting of three murine AAA tissues using the Loupe Cell Browser software. These cells were clustered into five groups based on k-means clustering and are shown in the t-SNE plot (Figure 4a). The heatmap demonstrated that these clusters differed significantly in transcriptome (Figure 4b). These clusters were mapped to different types of vascular cells based on their highly expressed genes (Table 2). As illustrated by t-SNE plots, CANX was obviously expressed in endothelial cells, smooth muscle cells (SMCs), fibroblasts, and macrophages (Figure 5a). CD44 was mainly expressed in SMCs and a subcluster of fibroblasts (Figure 5b). DAXX and STAT1 were both expressed in some clusters of endothelial cells, SMCs, fibroblasts, and macrophages but their expression levels were relatively low (Figure 5c, 5d). However, AR and TP53 were not detected within any cluster of vascular cells.

Figure 4.

(a) t-SNE plot of 5 clusters corresponding to different cells types in AAA and (b) Heatmap of highly expressed genes among each cluster. t-SNE, t-distributed stochastic neighbor embedding; AAA, abdominal aortic aneurysm.

Table 2.

Highly expressed genes within different cell types.

Cell type	Highly expressed genes	Fold change
Fibroblast	ANGPTL1	6.52
	SFRP4	5.78
	DCN	5.73
	MMP3	5.71
	CLEC3B	5.71
	ANGPTL7	5.71
	CILP	5.66
	CTHRC1	5.61
	GDF10	5.34
	FBLN1	5.32
Erythroid	HBA-A2	7.49
	HBA-A1	7.43
	HBB-BT	7.15
	ALAS2	6.82
	HBB-BS	6.80
	SNCA	6.73
	BPGM	4.75
	UBE216	4.10
	MKRN1	3.54
	CD24A	2.79
Smooth muscle cell	SGCG	6.93
	KCNMB1	6.71
	SUSD5	6.64
	MYH11	6.57
	PTPRZ1	6.34
	NPNT	6.32
	SOST	6.31
	CNN1	6.29
	SYNPO2	6.09
	OPTC	6.00
Endothelial cell	SOX17	8.67
	SOX18	8.64
	MMRN1	8.58
	PODXL	8.30
	BTNL9	8.21
	CLDN5	8.16
	MMRN2	8.13
	CCL21A	8.12
	MYCT1	8.04
	TMEM88	7.95
Macrophage	MZB1	10.52
	JCHAIN	9.73
	BCL2A1B	9.13
	LY86	8.75
	CD84	8.72
	MS4A6B	8.62
	FCGR2B	8.43
	MS4A4B	8.36
	FCMR	8.36
	CD74	8.34

Figure 5.

The t-SNE plot of expression of CANX (a), CD44 (b), DAXX (c), and STAT1 (d) within all clusters in AAA. t-SNE, t-distributed stochastic neighbor embedding; AAA, abdominal aortic aneurysm.

Discussion

The most significant features of AAA include inflammation and ECM remodeling. The inflammatory cells produce proteinases and MMPs, which promote ECM remodeling by degrading the collagen fibers and elastin. The matrix can then be cleaved into short fragments resembling bioactive chemokines, which recruit immune cells to active the immune response.[20] Thus, crosstalk between leukocytes and mesenchymal cells is mediated, establishing an interactive circle of inflammation and ECM remodeling. The current study identified some new candidate genes and signaling pathways closely associated with AAA using an integrated informatics analysis.

The GO and KEGG analyses revealed that the DEGs identified in the current study were enriched in pathways related to inflammation and ECM remodeling. A previous microarray-based expression profiling of AAA and non-aneurysmal tissues revealed that the complement cascade pathway was significantly altered in AAA.[10] A recent study demonstrated that DEGs were enriched in proteolysis, inflammation, and apoptotic processes.[21] Proteoglycan syndecan-mediated signaling, glypican pathway, and Alpha9 beta1 integrin signaling may be involved in adhesion and signaling conduction between cells and ECM. This aberrant expression of collagens and fibronectins induces myofibroblast transition, activates adhesion molecules, and enhances macrophage infiltration.[20] The TRAIL pathway is linked to cell apoptosis and matrix degradation; TRAIL and its receptors show increased expression in AAA, accompanied by vascular calcification.[22] Osteoprotegerin (OPG), the TRAIL inhibitor, protects against AAA formation by decreased expression of MMP9 in vascular smooth cells.[23] Increased vascular endothelial growth factor (VEGF) and interferon (IFN)-gamma signaling may be attributed to enhanced angiogenesis and transmural macrophage migration.[24] The ErbB pathway, also known as epidermal growth factor receptor (EGFR) signaling, plays a role in cell proliferation, differentiation, and survival. Activation of EGFR affects the physiology and pathophysiology of the cardiovascular system. AngII can induce AAA formation by activating ErbB signaling in SMCs mediated by ADAM17, and erlotinib, an EGFR inhibitor, can protect mice from AAA formation induced by AngII.[25,26] Inhibition of EGFR activity is emerging as a potential therapeutic strategy to treat AAA.

In accordance with previous reports, some proinflammatory cytokines, chemokines, proteolytic proteins, and aneurysm-related GO entries were identified in the pathogenesis and development of AAA.[27,28] By using the emerging single-cell sequencing technique, we further explored the cell-specific expression of these hub genes. The PPI and single-cell sequencing analysis identified four hub genes involved in the progression of AAA. CD44, a polymorphic hyaluronate receptor, may participate in chronic inflammation. The level of CD44 was high and positively correlated with macrophage content in pathological tissues. Proinflammatory factors can induce shedding of CD44 from macrophages, and soluble CD44, in turn, stimulates endothelial cells to secrete interleukin (IL)-1β, enhancing local inflammatory responses.[29] Another report confirmed that CD44 was highly expressed in human AAA tissues, based on bioinformatics and quantitative PCR assays.[30] In addition, osteopontin was shown to have a pro-autophagy effect on vascular SMCs mediated by the integrin/CD44 pathway.[31] In accordance with our t-SNE plot, CD44 was expressed in SMCs and macrophages. These results suggested that CD44 may act on macrophages to enhance inflammation responses and induce loss of SMCs by proapoptotic effects. Meanwhile, the role of CD44-positive T cells in AAA is worth investigating further. CANX encodes a member of the calnexin family. As endoplasmic reticulum-associated proteins, calnexins can, upon binding with calcium, interact with newly synthesized glycoproteins, facilitating protein folding. CANX may participate in regulating endoplasmic reticulum stress and unfolding protein response for cell survival as a universal chaperone.[32] CANX also interacts with STAT,[33] another hub gene in our study. A report showed that CANX interacted with NOX4, a NADPH oxidase, to form a macromolecule required for maturation and function of NOX4 in the endoplasmic reticulum.[34] However, there are no reports about its role in AAA. Based on our results, we speculated that CANX regulates oxidative stress in vascular cells in AAA progression. The signal transducers and activators of transcription (STAT) are a family of TFs associated with numerous cytokines, growth factors, and interferon receptors, including seven members (STAT-1, -2, -3, -4, -5α, -5β, -6). Activation of STAT1 is involved in regulation of apoptosis and matrix remodeling. Increased expression of STAT1 mRNA and protein occurred in AAA tissues, with its highest activity in adventitial inflammatory cells. The expression of STAT1 was also high in peripheral blood cells of AAA. These local and systemic changes indicate that STAT1 is closely associated with AAA.[35] Hypoxia and the subsequent inflammation have been shown to induce overexpression of MCP-1 in aortic SMCs, which activate macrophages to secrete IL-6, which in turn promotes aortic SMC apoptosis through STAT1.[36] In addition, STAT1 interacts with the toll-like receptor 2[37] and nuclear factor-kappaB[38] signaling pathway to regulate inflammation and vascular remodeling of AAA; in contrast, loss of STAT1 is linked to an increased chance of AAA rupture.[17] These reports, combined with our results, suggest that STAT1 has a complex role in AAA that requires further investigation. Another hub gene, DAXX, serves as a histone chaperone and may participate in the epigenetic modification of vascular cells of AAA.[39,40] Even though our single cell analysis revealed that DAXX was expressed in all vascular cells, this result must be validated in additional in vitro and in vivo experiments.

Conclusions

By combining a microarray profile and a single-cell RNA-seq dataset of AAA tissues, we identified four hub genes, CANX, CD44, STAT1, and DAXX, and some signaling pathways closely related to AAA. These key genes provide insight into the mechanism of AAA initiation and progression and thus are potential therapeutic targets for AAA.