Identification of crucial genes mediating abdominal aortic aneurysm pathogenesis based on gene expression profiling of perivascular adipose tissue by WGCNA.

BACKGROUND: With a mortality rate of 65-85%, a ruptured abdominal aortic aneurysm (AAA) can have catastrophic consequences for patients. However, few effective pharmaceutical treatments are available to treat this condition. Therefore, elucidating the pathogenesis of AAA and finding the potential molecular targets for medical therapies are vital lines of research.
METHODS: An mRNA microarray dataset of perivascular adipose tissue (PVAT) in AAA patients was downloaded and differentially expressed gene (DEG) screening was performed. Weighted gene co-expression networks for dilated and non-dilated PVAT samples were constructed via weighted correlation network analysis (WGCNA) and used to detect gene modules. Functional annotation analysis was performed for the DEGs and gene modules. We identified the hub genes of the modules and created a DEG co-expression network. We then mined crucial genes based on this network using Molecular Complex Detection (MCODE) in Cytoscape. Crucial genes with top-6 degree in the crucial gene cluster were visualized, and their potential clinical significance was determined.
RESULTS: Of the 173 DEGs screened, 99 were upregulated and 74 were downregulated. Co-expression networks were built and we detected 6 and 5 modules for dilated and non-dilated PVAT samples, respectively. The turquoise and black modules for dilated PVAT samples were related to inflammation and immune response. MAP4K1 and PROK2 were the hub genes of these 2 modules, respectively. Then a DEG co-expression network with 112 nodes and 953 edges was created. PLAU was the crucial gene with the highest connectivity and showed potential clinical significance.
CONCLUSIONS: Using WGCNA, gene modules were detected and hub genes and crucial genes were identified. These crucial genes might be potential targets for pharmaceutic therapies and have potential clinical significance. Future in vitro and in vivo experiments are required to more comprehensively explore the biological mechanisms by which these genes affect AAA pathogenesis. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Abdominal aortic aneurysm (AAA) is the true aneurysm of the abdominal aorta and the diameter of AAA is over 1.5 times than normal (1). Risk factors for AAA include smoking, male sex, family history, old age, and other factors (2). An AAA rupture has catastrophic consequences for patients, and the mortality rate is 65–85%. Currently, the primary treatment for AAA includes open surgical repair and endovascular repair (3); however, few effective pharmaceutical treatments are available. Therefore, elucidating the pathogenesis of AAA and finding the potential molecular targets for medical therapies are quite important. Previous studies have indicated that inflammation, immune response, extracellular matrix (ECM) degradation, vascular smooth muscle phenotypic switching, etc., promote the development and progression of AAA (4-6).

Perivascular adipose tissue (PVAT) is the adipose tissue surrounding most blood vessels (7). It represents approximately 3% of total adipose tissue and was initially thought only to support the vasculature, although it has now been shown to play important roles in vascular pathologies as well, such as in atherosclerosis, hypertension, and AAA (8-10). Dysfunctional PVAT is associated with inflammation through the production of inflammatory cytokines and chemokines (11). Piacentini et al. were the first to use a high throughput approach via an mRNA microarray to examine PVAT in AAA. Through functional enrichment analysis, they found that an altered immunological response in PVAT is involved in the pathogenesis of AAA (12). However, weighted correlation network analysis (WGCNA) has rarely been used with high throughput PVAT data to investigate crucial genes that mediate AAA pathogenesis.

WGCNA is an algorithm first created by Zhang and Horvath (13). Using gene expression data, WGCNA is used to build a weighted gene co-expression network, detect gene modules, and identify genes with critical biological functions (14). Non-coding RNA data, proteomic data, metabolomic data, and even microbiome data can be used in WGCNA analyses (15-17). An R package has also been developed by Langfelder and Horvath (18) and is available for researchers on the R official website (https://cran.r-project.org/) to perform WGCNA analysis.

In our study, mRNA microarray data for the PVAT of AAA patients was obtained from the Gene Expression Omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo/). This includes data from dilated abdominal aorta PVAT (dilated PVAT) samples and non-dilated abdominal aorta PVAT (non-dilated PVAT) samples. We constructed weighted networks and detected gene modules for both dilated PVAT and non-dilated PVAT samples. Functional annotation analysis was conducted using genes in modules. We then mapped differentially expressed genes (DEGs) into the co-expression network of dilated PVAT samples and created a DEG co-expression network to identify crucial genes. The crucial genes with top-6 degree in the crucial gene cluster were visualized, and their potential clinical significance was determined. We present the following article in accordance with the MDAR reporting checklist (available at: http://dx.doi.org/10.21037/atm-20-3758).

Methods

Medical ethics

The raw datasets were available from the GEO database (http://www.ncbi.nlm.nih.gov/geo/; GSE119717). In our study, neither human trials nor animal experiments were executed.

Data preprocessing

We selected the GEO dataset GSE119717 for our study. The series matrix files and data table of the GPL10558 platform were obtained from the GEO website (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119717). The probe ID of the probe matrix from the series matrix files was annotated with gene symbols using the platform data table, then the gene expression matrix was generated. Samples from dilated PVAT and non-dilated PVAT were included for further analysis.

DEG screening

The complete differentially expressed (DE) transcripts list was obtained from the online-only data supplement 3 of the AAA PVAT study conducted by Piacentini et al. This list is based on paired-sample tests. Confounder effects and patient heterogeneity were eliminated (12). The threshold to select DE transcripts and DEGs in our study was a |log2(Fold-change)| of >0.5 (dilated PVAT/non-dilated PVAT) and an adjusted P value of <0.05.

Weighted gene co-expression network construction

The top 25% of genes with the highest expression variance were selected to construct the co-expression network. Pearson’s correlation was used to calculate the correlation coefficients between genes and produce a correlation matrix. The cutoff R2 was set as 0.9 for the soft threshold β in both dilated PVAT and non-dilated PVAT samples. The adjacency matrix, which stored co-expression network information, was then built using the following formula: amn = |cmn|β (amn: element in the adjacency matrix; cmn: correlation coefficient between genes; β: soft threshold). A topological overlap measure (TOM) matrix was then generated based on the adjacency matrix. Using the TOM matrix, average linkage hierarchical clustering and a dynamic tree cut were performed to detect gene modules. The minimal module size was set as 100 and similar modules were merged with a threshold of 0.25.

Functional annotation analysis

Upregulated genes and downregulated genes were used to conduct a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and a Gene Ontology Biological Processes (GO-BP) analysis using The Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 (https://david.ncifcrf.gov/) (19). These analyses were also performed for each gene module. The threshold for the functional annotation analysis was set as P<0.05 with a count number of >2. The patterns of co-expression between dilated PVAT samples and non-dilated PVAT samples were then compared using KEGG pathway and GO-BP analysis.

Hub gene and crucial gene mining

The intramodular connectivity within gene modules was determined by the WGCNA algorithm. The gene with the highest connectivity in each module was considered the hub gene. A DEG co-expression network was then created by mapping DEGs into the co-expression network of dilated PVAT samples using Cytoscape 3.7.2. The Molecular Complex Detection (MCODE) algorithm was used within Cytoscape to detect crucial gene clusters based on the DEG co-expression network. Crucial genes with top-6 degree in the crucial gene cluster were visualized, and their potential clinical significance was determined by receiver operating characteristic (ROC) analysis.

Statistical analysis

R v3.6.2 was used to conduct data preprocessing, DEG screening and WGCNA analysis. DAVID v6.8 was used to perform functional annotation analysis, and MCODE in Cytoscape v3.7.0 was selected to mine crucial genes. We have described the details of these bioinformatic analyses in corresponding subsections. IBM SPSS 25.0 was used to conduct ROC analysis, which showed the potential diagnostic value of crucial genes. A P value of <0.05 was considered significant.

Results

Study workflow

The flowchart of our study is shown in . Data were downloaded from the GEO website and preprocessed. DEGs were then screened and weighted co-expression networks were constructed. KEGG pathway and GO-BP analyses were conducted for the DEGs and genes in each module. The co-expression patterns between dilated PVAT samples and non-dilated PVAT samples were then compared based on the results of the GO-BP and KEGG analyses. A DEG co-expression network was then constructed by mapping DEGs into the co-expression network of dilated PVAT samples. Crucial genes were identified by MCODE using the DEG co-expression network. The potential diagnostic value of crucial genes with top-6 degree was determined by ROC analysis.

Figure 1

Workflow for the whole study. DEG, differentially expressed gene.

In our study, 194 DE transcripts were screened with an adjusted P value of <0.05 and a |log2 (Fold-change)| of >0.5 () as the threshold. These DE transcripts corresponded to 173 DEGs, of which 99 genes were upregulated and 74 were downregulated (). The 10 DEGs with the highest | log2(Fold-change) | are shown in .

Figure 2

Patterns of DE transcripts and DEGs. (A) The volcano plot showed the number of DE transcripts. The red dots represent the up-regulated transcripts while the green dots represent the down-regulated transcripts. (B) The heatmap showed the expression patterns of DEGs between dilated PVAT samples and non-dilated PVAT samples. DE, differentially expressed; DEG, differentially expressed gene; PVAT, perivascular adipose tissue.

Table 1

DEGs with top-10 |log2(Fold-change)| (dilated PVAT/non-dilated PVAT)

Gene symbol	Official full name	Log2(Fold-change)	Adjusted P value
Up-regulated
FOS	Fos proto-oncogene, AP-1 transcription factor subunit	3.98	1.20E-14
FOSB	FosB proto-oncogene, AP-1 transcription factor subunit	3.68	3.04E-11
EGR1	Early growth response 1	2.47	1.94E-09
CYR61	Cellular communication network factor 1	2.11	1.01E-08
DUSP1	Dual specificity phosphatase 1	2.01	1.94E-09
ZFP36	ZFP36 ring finger protein	1.82	1.01E-08
JUNB	JunB proto-oncogene, AP-1 transcription factor subunit	1.41	1.01E-08
JUN	Jun proto-oncogene, AP-1 transcription factor subunit	1.35	7.60E-07
NR4A2	Nuclear receptor subfamily 4 group A member 2	1.21	2.18E-04
S100A9	S100 calcium binding protein A9	1.19	8.30E-04
Down-regulated
SLN	Sarcolipin	−1.76	4.73E-07
DES	Desmin	−1.54	2.67E-05
MYH11	Myosin heavy chain 11	−1.38	2.97E-05
CNN1	Calponin 1	−1.23	8.71E-05
MYLK	Myosin light chain kinase	−1.17	2.63E-04
KCNMB1	Potassium calcium-activated channel subfamily M regulatory beta subunit 1	−1.11	2.78E-05
PDLIM3	PDZ and LIM domain 3	−1.06	8.87E-05
SYNM	Synemin	−1.03	7.81E-05
NTS	Neurotensin	−0.97	1.54E-03
PDLIM3	PDZ and LIM domain 3	−0.96	5.36E-05

PVAT, perivascular adipose tissue.

Sample clustering detected no outliers in either dilated or non-dilated PVAT samples, and all were included in the analysis ( and Figure S1A). With a cutoff R2 value of 0.9, a soft threshold β value of 8 was selected for dilated PVAT samples, while a value of 10 was selected for non-dilated PVAT samples ( and Figure S1B,S1C). The histogram and the linear plot showed that we constructed 2 networks with scale-free topology ( and Figures S1D,S1E). After a dynamic tree cut and merging, we discovered 6 and 5 gene modules for dilated PVAT and non-dilated PVAT sample, respectively.

Figure 3

Weighted co-expression network construction for dilated PVAT. (A) No outliers were detected by sample clustering. (B,C)The soft-threshold β was chosen to be 8. (D,E) The network we constructed met the requirements of scale-free topology. (F) After dynamic tree cut and merging, 6 gene modules were detected. PVAT, perivascular adipose tissue.

The upregulated and downregulated DEGs were used in functional annotation analyses; the results are shown in Table S1 and Table S2. Upregulated DEGs were related to inflammation and the immune response; downregulated DEGs were related to muscle contraction. The 5 KEGG pathways and GO-BP terms with the most significant p values were visualized ().

Figure 4

GO-BP and KEGG pathways with top-5 significant p value for up-regulated DEGs and Down-regulated DEGs. (A,B) Up-regulated DEGs are mainly associated with inflammation and immune response. (C,D) Down-regulated DEGs are mainly associated with muscle contraction. The left axis represents count number and the right axis represent −logP. GO-BP, Gene Ontology Biological Processes; DEG, differentially expressed gene.

The results of the KEGG pathway and GO-BP analyses of gene modules are presented in Table S3 and Table S4. The black and turquoise gene modules for dilated PVAT samples were associated with inflammation and the immune response (). However, these inflammation and immune response-associated pathways were scattered across the modules of non-dilated samples.

Figure 5

GO-BP terms and KEGG pathways with top-5 significant p value for black and turquoise module of dilated PVAT. The left axis represents count number and the right axis represent −logP. GO-BP, Gene Ontology Biological Processes; KEGG, Kyoto Encyclopedia of Genes and Genomes; PVAT, perivascular adipose tissue.

Hub genes for each module are listed in Table S5. MAP4K1 and PROK2 were the respective hub genes for the turquoise and black PVAT sample gene modules.

After removing isolated nodes and node pairs, a DEG co-expression network with 112 nodes and 953 edges was obtained (). Crucial genes were mined using MCODE, and crucial genes with top-6 degree were visualized (). Plasminogen activator urokinase (PLAU) was the crucial gene with the highest connectivity. All these crucial genes showed potential clinical significance ( and ).

Figure 6

DEG co-expression network and crucial gene cluster. (A) The DEGs co-expression network with 112 nodes and 953 edges were generated. The red box represents the up-regulated genes and the green box represents the down-regulated genes. (B) Crucial genes with top-6 degree were visualized. DEG, differentially expressed gene.

Figure 7

The ROC curve of the crucial genes with top-6 degree. ROC, receiver operating characteristic.

Table 2

AUC of the crucial genes with top-6 degree

Gene symbol	Official full name	AUC	P value
PLAU	Plasminogen activator, urokinase	0.819	<0.001
CD14	CD14 molecule	0.741	0.001
SPON1	Spondin 1	0.673	0.021
TCEAL1	transcription elongation factor A like 1	0.670	0.024
CD248	CD248 molecule	0.720	0.003
CCL2	C-C motif chemokine ligand 2	0.789	<0.001

Discussion

Obesity has become a major public health concern worldwide and significantly increases the risk of AAA (20). PVAT expands and becomes dysfunctional during obesity, which contributes to AAA progression (10,21). By using gene co-expression networks, our study aimed to provide further insight into the associations between obesity and AAA pathogenesis. We selected 173 DEGs, of which 99 were upregulated and 74 were downregulated. The upregulated DEGs were related to inflammation and the immune response, while downregulated DEGs were related to muscle contraction. By constructing weighted gene co-expression networks for dilated and non-dilated PVAT samples, 6 gene modules were discovered for dilated samples and 5 were discovered for non-dilated samples. The black and turquoise modules for dilated PVAT samples were also associated with inflammation and the immune response, and the hub genes of these 2 modules were PROK2 and MAP4K1, respectively. A DEG co-expression network was then created, and the crucial genes were identified. Crucial genes with top-6 degree in the crucial gene cluster were visualized. PLAU was the gene with the highest connectivity. Finally, by conducting a ROC analysis, we found that all 6 crucial genes had potential diagnostic significance.

Piacentini et al. were the first to use an mRNA microarray to study PVAT in human AAA. Using a |log2(Fold-change)| of >0.38 and an adjusted P value of <0.05, they screened 193 upregulated DE transcripts and 142 downregulated transcripts (12). We chose to use a |log2(Fold-change)| of >0.5 and an adjusted P value of <0.05 as a more stringent threshold, and we screened fewer DE transcripts and DEGs than Piacentini et al.

The results of our functional annotation analysis showed that the co-expression patterns between dilated and non-dilated PVAT were different. Inflammation and the immune response were clustered in the black and turquoise modules of dilated PVAT samples. Using microarray data from AAA tissue, researchers have found that immune and inflammation-related processes, such as B cell and T cell signaling pathways, cytokine–cytokine receptor interactions, natural killer (NK) cell–mediated cytotoxicity, etc. are involved in AAA pathogenesis (22-24). Immune cells such as macrophages, B and T cells, NK cells, etc. are present in AAA tissue (25-27). Studies investigating PVAT tissue have made similar discoveries. PVAT is rich in interleukin-6 (IL-6), which might upregulate CD14 expression in the adventitia throughout the formation of an AAA (28,29). Li et al. found that endothelin-1 increases reactive oxygen species production, as well as CD4+ T cell and monocyte/macrophage infiltration in Apoe(−/−) mouse models (30). Folkesson et al. also discovered that increased T cell, neutrophil, and mast cell counts in the PVAT of AAA patients (31). In addition, PAVT contributes to AAA formation by promoting angiotensin II type 1 receptor and platelet-derived growth factor D–mediated vascular inflammation and the binding of IL-18 to corresponding receptors (32-34). These studies, in addition to the research by Piacentini et al., suggest that PVAT mediates the pathogenesis of AAA via inflammation and the immune response; our results were consistent with these findings. Therefore, targeting inflammation and modulating the immune response at the PVAT level, might inhibit AAA progression. AMP-activated protein kinase (AMPK) activators, such as diosgenin, salicylate, metformin, etc., can inhibit inflammation induced by palmitic acid and downregulate the expression of pro-inflammatory cytokines CCL2, IL-6, and TNF-α in PVAT (35,36). Additionally, Skiba et al. found that AVE0991, an agonist of the Mas receptor, can reduce the recruitment of monocytes/macrophages to PVAT, and Mikolajczyk et al. showed that 1,2,3,4,6-penta-O-galloyl-β-d-glucose (PCG) can decrease leukocyte and T cell infiltration in PVAT (37,38). Although few studies of PVAT in AAA have been carried out, these findings still suggest that anti-inflammatory and immunomodulating strategies might be effective treatments for AAA. We also found that downregulated DEGs were primarily associated with muscle contraction. Although the detailed function of muscle contraction in AAA pathogenesis has not yet been verified, Chew et al. found that certain proteases, such as matrix metalloproteinase and elastase, have inhibitory roles in vascular smooth muscle contraction by inhibiting the passage of Ca2+, thereby promoting early aortic dilation (39-41). Our previous study that integrated mRNA microarray data from AAA tissues obtained similar findings (42).

PROK2 and MAP4K1 were the hub genes of the black and turquoise gene modules for dilated PVAT samples, respectively. PROK2 is a chemokine that interacts with 2 types of G protein-coupled receptors (PKR1 and PKR2) and is associated with inflammation and the immune response (43). PROK2 is upregulated in granulocytes and macrophages found in inflamed tissue (44). Previous studies suggest that PROK2 plays important roles in mediating inflammation-associated disorders, such as inflammatory pain, psoriasis, and orchitis (44-46). Choke et al. also found that PROK2 was upregulated with 7-fold at rupture sites of AAA (47). MAP4K1 is usually expressed by T and B cells and participates in the regulation of inflammation and immune-related processes, such as cytokine production, immune cell adhesion, and T and B cell activation et al. Our previous study also found that MAP4K1 is the hub gene of the module related to inflammation and the immune response in AAA samples (42,48-50). Although the detailed mechanisms of these 2 PVAT genes in AAA remain unknown, these genes might be potential targets for the mediation of AAA pathogenesis because they are important mediators of inflammation and the immune response.

PLAU, also named u-PA, is the gene with the highest connectivity among the crucial genes. It is a serine protease that cleaves plasminogen into plasmin. Once generated, plasmin activates metalloproteinases (MMPs), which promote the formation and progression of AAA by degrading ECM (1,2,51). Therefore, u-PA is closely associated with the pathogenesis of AAA. Using Apoe(−/−) Plau(−/−) mouse models, Carmeliet et al. found that aneurysm formation and ECM degradation were inhibited when PLAU genes were deficient (52). These studies suggest that u-PA expression in the PVAT of AAA patients might mediate AAA pathogenesis by indirectly degrading ECM via plasmin and MMPs. Furthermore, Qian et al. had locally overexpressed human plasminogen activator inhibitor-1 (PAI-1), a natural inhibitor of u-PA, in the adventitia of AAA mouse models and found that AAA formation and early progression were prevented. Reduced activity of MMPs and inflammatory cells were also observed. The blockage of u-PA/u-PA receptor (u-PAR) mediated inflammation by PAI-1 might account for the therapeutic effects observed when PAI-1 is locally overexpressed (53). The findings of Qian et al. implicated that local inhibition of u-PA/u-PAR might be a new medical therapy to prevent AAA formation and stabilize early stage AAA.

Other genes with high degree in the crucial gene cluster, such as CD14 and CCL2, also play important roles in AAA pathogenesis. CD14 is preferentially expressed by monocytes and macrophages and, as previously discussed, was found to be upregulated in the PVAT of AAA patients. CD14(−/−) mice have been shown to exhibit both reduced macrophage infiltration and AAA formation (28). CCL2, also called monocyte chemoattractant protein 1 (MCP1), is a C-C motif chemokine ligand and participates in inflammation and immune-related processes. It performs chemokine functions for monocytes and basophils and is upregulated in AAA (54-56). These crucial genes are also potential molecular targets to mediate AAA pathogenesis and have potential diagnostic value.

Our study was the first to construct weighted gene co-expression networks for dilated and non-dilated PVAT samples from AAA patients. We also combined DEGs and the weighted co-expression network to identify crucial genes with potential biological functions in the pathogenesis of AAA. However, we did not analyze the correlation between gene modules and the clinical traits of patients due to the clinical information provided on the GEO website being incomplete.

The co-expression network patterns found between the dilated PVAT samples and non-dilated PVAT samples were different. These differences imply that inflammatory and immunological responses in PVAT are involved in the pathogenesis of AAA Then crucial genes were identified, which might potentially play critical roles in AAA formation.

Conclusions

Using WGCNA, 6 and 5 gene modules were detected for dilated PVAT and non-dilated PVAT samples, respectively. Module hub genes were detected, and then crucial genes were identified using our DEG co-expression network. These crucial genes might be potential molecular targets for pharmaceutical therapies and have potential diagnostic value for AAA. Future in vitro and in vivo experiments are required to clarify the detailed molecular mechanisms by which these genes contribute to the pathogenesis of AAA.

The article’s supplementary files as