Identification of crucial genes in abdominal aortic aneurysm by WGCNA

Introduction

Abdominal aortic aneurysm (AAA) is a type of true aneurysm located at the abdominal aorta, and it is defined as permanent and irreversible dilation of the abdominal aorta to 50% more than the normal aortic diameter (Sakalihasan et al., 2018). The prevalence of AAA in populations above 65 years old is approximately 4–8%, and it is the main cause of death in the senior population. Furthermore, the mortality of a ruptured AAA is almost 100% if left untreated (Nordon et al., 2011; Sakalihasan et al., 2018). Open surgical repair and endovascular repair remain two main methods to treat AAA. However, there are temporarily few effective medical therapies to treat AAA (Golledge, 2019). Therefore, elucidating the molecular mechanism of AAA pathogenesis and exploring the potential molecular targets of medical therapies for AAA is of vital importance. Although the molecular mechanism of AAA pathogenesis is not yet clear, it may involve biological processes such as immune response, chronic inflammation, oxidative stress, phenotypic switching of smooth muscle, and degradation of extracellular matrix, etc. (Golledge et al., 2006; Raffort et al., 2017; Weintraub, 2009). Abundant studies have been performed on the molecular mechanism of AAA pathogenesis (Golledge, 2019; Raffort et al., 2017; Sakalihasan et al., 2018), but weighted gene co-expression network analysis (WGCNA) has not been used to construct a gene co-expression network in AAA.

Weighted gene co-expression network analysis was developed by Horvath & Zhang in 2005 (Langfelder & Horvath, 2008; Zhang & Horvath, 2005), and an R package is available on the official website of R (https://cran.r-project.org/). WGCNA can be used to construct a weighted gene co-expression network, detect gene modules, correlate gene modules with clinical traits and identify intramodular hub genes (Langfelder & Horvath, 2008; Van Dam et al., 2018; Zhang & Horvath, 2005). Furthermore, WGCNA can be used not only to construct co-expression networks based on coding RNA expression data but also to construct noncoding RNA co-expression networks, such as miRNA co-expression networks and lncRNA co-expression networks (Hu et al., 2019; Ma et al., 2019).

Previous studies using microarray data of AAA have conducted differentially expressed gene (DEG) analysis and functional enrichment analysis. Gene expression profiles were compared between AAA and normal abdominal aorta, and validated genes such as CTLA4, NKTR and CD8A were identified in DEG analysis (Biros et al., 2015; Lenk et al., 2007). Nevertheless, DEGs and functional enrichment analysis cannot reveal connections and interactions among genes that are crucial in biological processes.

In our study, three Gene Expression Omnibus (GEO) mRNA microarray datasets (Table 1) were used to conduct WGCNA. Gene co-expression networks and gene modules were constructed in both AAA expression data and normal expression data. Due to the lack of clinical trait data, we used genes in all gene modules to conduct functional enrichment analysis and then mapped the DEGs to whole gene expression networks of both the AAA and normal expression data. Furthermore, we identified 10 hub genes and showed their diagnostic values by receiver operating characteristic (ROC) analysis.

Materials and Methods

Results

Workflow

The workflow is shown in Fig. 1. DEG identification was conducted, and 143 DEGs were screened. WGCNA on the AAA and normal samples was performed. Co-expression networks for the AAA and normal conditions were constructed, and gene modules were also detected for further hub gene identification and comparison of expression patterns between the AAA and normal samples. Then, the DEGs were mapped to the AAA co-expression network and a DEG co-expression network was obtained. Crucial genes were identified based on this network, and ROC analysis was conducted using these crucial genes.

Identification of DEGs

A total of 143 DEGs were identified with the threshold at |log₂(fold-change)| > 0.8 and p < 0.05 and consisted of 125 down-regulated genes and 18 up-regulated genes. The expression patterns of the DEGs are shown in Fig. 2. Down-regulated genes and up-regulated genes with a top 10-fold change ratio are shown in Table S1.

Construction of the co-expression network of the AAA and normal conditions

The sample clustering dendrograms of the AAA and normal conditions are shown in Figs. 3A and 4A, respectively. In both conditions, all samples were included in the clusters. The soft-power threshold β was determined by the function “sft$powerEstimate”; β = 5 and β = 16 were selected for further analysis of the AAA and normal conditions, respectively (Figs. 3B–3E; Figs. S2B–S2E). Then, gene modules were detected based on the TOM matrix. In the analysis, six modules were detected under the AAA condition, while 10 modules were detected under the normal condition (Figs. 3F and 3G; Figs. S2F and S2G).

Mining crucial genes mediating AAA

Differential expressed genes were mapped to the whole co-expression network in the AAA condition using Cytoscape, and a DEG co-expression network was obtained. The threshold of weighted edges was set as 0.2 and 102 nodes, and 303 edges are in the DEG co-expression network (Fig. 5A). Then, GO-BP and KEGG analyses were conducted using genes from the DEG co-expression network, and the results are shown in Table S7. These genes were mainly associated with oxidative stress and inflammation. The most significant cluster of the DEG co-expression network detected by MCODE consisted of 10 genes and was visualized (Fig. 5B; Table 3). Yip1 domain family member 6 (YIPF6), a protein associated with clathrin-derived vesicle budding and vesicle-mediated transport, was the gene with the highest degree in the cluster.

Table 3:

Genes in the hub gene cluster of DEG co-expression network.

	Gene ID	Gene symbol	Official full name
1	286451	YIPF6	Yip1 domain family member 6
2	23637	RABGAP1	RAB GTPase activating protein 1
3	157567	ANKRD46	Ankyrin repeat domain 46
4	23171	GPD1L	Glycerol-3-phosphate dehydrogenase 1 like
5	10424	PGRMC2	Progesterone receptor membrane component 2
6	25994	HIGD1A	HIG1 hypoxia inducible domain family member 1A
7	2762	GMDS	GDP-mannose 4,6-dehydratase
8	4256	MGP	Matrix Gla protein
9	291	SLC25A4	Solute carrier family 25 member 4
10	116496	FAM129A	Niban apoptosis regulator 1

DOI: 10.7717/peerj.7873/table-3

All 10 genes were selected for ROC analysis. The area under the curve (AUC) of each gene is listed in Table 4. HIG1 hypoxia inducible domain family member 1A (HIGD1A) has the highest AUC value (AUC = 0.886, SE = 0.037, p < 0.001), while matrix Gla protein (MGP) has the lowest AUC value (AUC = 0.742, SE = 0.050, p < 0.001). The ROC curves of the genes with AUC exceeding 0.8 were plotted (Fig. 6). These genes were YIPF6, PGRMC2, ANKRD46, GPD1L, HIGD1A, GMDS, FAM129A and SLC25A4. Thus, these genes have potential diagnostic value and may become biomarkers for AAA; however, they need to be further validated in future studies.

Table 4:

The area under the curve (AUC) of 10 genes in the cluster detected by MCODE.

Gene symbol	AUC	p-value
YIPF6	0.807	<0.001
RABGAP1	0.797	<0.001
PGRMC2	0.845	<0.001
ANKRD46	0.832	<0.001
GPD1L	0.831	<0.001
HIGD1A	0.866	<0.001
GMDS	0.845	<0.001
MGP	0.742	<0.001
FAM129A	0.830	<0.001
SLC25A4	0.811	<0.001

DOI: 10.7717/peerj.7873/table-4

Discussion

In our present study, 143 DEGs were screened (Fig. 2; Table S1). WGCNA was used to construct a co-expression network and detect gene modules in both AAA and normal conditions. In our analysis, six modules in the AAA condition and 10 modules in the normal condition were detected based on the co-expression network. The gene modules were subjected to GO-BP and KEGG pathway analyses, and co-expression patterns between the AAA and normal conditions were compared (Fig. 4; Tables S2 and S3). In the AAA condition, the black module was mainly associated with the immune response and T/B cell activation, while the turquoise module was mainly associated with the inflammatory response and signaling pathways in inflammation. The red module was mainly associated with muscle contraction and focal adhesion. In the normal condition, the GO-BP terms and KEGG pathways associated with immune response, inflammatory response, muscle contraction and focal adhesion were scattered in each module and were not enriched in one or two modules. These differences suggest that the immune response, inflammation and vascular smooth muscle play important roles in the pathogenesis of AAA. For each module, the genes with the highest intramodular connectivity were selected as the hub genes (Table S6). Then, DEGs were mapped to the co-expression network of AAA, and a DEG co-expression network with 102 nodes and 303 edges was obtained (Fig. 5A). The Hub gene cluster was the most significant cluster in the DEG co-expression network selected by MCODE. Ten genes were in the cluster: YIPF6, RABGAP1, ANKRD46, GPD1L, PGRMC2, HIGD1A, GMDS, MGP, FAM129A and SLC25A4 (Fig. 5B). These genes are potentially crucial genes in the pathogenesis of AAA and showed diagnostic value (Fig. 6).

A previous study using the same GEO datasets that we used screened 1,199 DEGs (Wan et al., 2018), which is far more than the number of DEGs we screened. However, unlike the previous study, we used log₂|FoldChange| > 0.8 as another threshold, and the adjusted p-value was used. Furthermore, the “sva” R package was used to normalize the three datasets before DEG analysis in our study. The differences in choosing parameters and normalization may account for the different number of DEGs screened.

Previous genome-wide association studies (GWASs) reported some significant AAA-associated Single-nucleotide polymorphism (SNPs) and genes from genomic data. Helgadottir et al. (2008) discovered one common variant (rs10757278-G on chromosome 9p21, OR = 1.31, p = 1.2 × 10⁻¹²) that is associated with the risk of AAA. Other SNPs, such as rs7025486 on 9q33 in the DAB2IP gene (OR = 1.21 and p = 4.6 × 10⁻¹⁰) (Gretarsdottir et al., 2010), rs1466535 on 12q13.3 in the LRP1 gene (OR = 1.15 and p = 4.52 × 10⁻¹⁰) (Bown et al., 2011), rs599839 on 1p13.3 in the SORT1 gene (OR = 0.81 and p = 7.2 × 10⁻¹⁴) (Jones et al., 2013) and rs6511720 in LDLR (OR = 0.76, p = 2.08 × 10⁻¹⁰) (Bradley et al., 2013), are significantly associated with AAA. Jones et al. (2017) conducted a meta-analysis of GWAS for AAA and identified 1q32.3 (SMYD2), 13q12.11 (LINC00540), 20q13.12 (near PCIF1/MMP9/ZNF335), and 21q22.2 (ERG) as four new loci associated with AAA and found that these four loci may be specifically associated with AAA compared with other cardiovascular diseases. We screened DEGs based on transcriptomic data, and there was some overlap between previous GWASs and DEGs we screened. MMP9 is in the DEG list with a fold-change of 2.03, which is in accord with the study by Jones et al. (2017) which found that one SNP near PCIF1/MMP9/ZNF335 was significantly associated with AAA. Furthermore, one SNP in LRP1 was associated with AAA, and Boucher et al. (2007) found that the LRP1 gene plays a protective role in vascular disease via TGF-β signaling (Bown et al., 2011). Thus, the TGFBR3 gene with a fold-change of −1.01 in our DEG list may be associated with LRP1 gene function and pathogenesis of AAA, although the detailed mechanism has not been clearly elucidated. Furthermore, we obtained genes closest to SNPs in those literature, respectively (i.e., GWAS genes) (Table S8) using dbSNP (https://www.ncbi.nlm.nih.gov/snp/) and found that there were two GWAS genes (GSN and MMP9) overlapped with DEGs. Then we compared the GWAS genes with genes enriched in GO-BP and KEGG pathway terms of DEGs in the co-expression network (Table S7) and there were also no overlap. Therefore, we inferred that GWAS genes were not “close” to DEGs in the co-expression network in the functional level. We also evaluated the association between SNPs and DEGs using eQTL calculator in GTEx website (https://gtexportal.org/home/testyourown) (Brown et al., 2017). Although some normalized effect size between some SNPs and genes were significant (Table S9), the biological significance is not clear and further validation is needed. Some researchers combined genomic data and patient history of AAA. Li et al. integrated genomic data of AAA based on whole genome sequencing and electronic health record data using a machine learning approach. Their study provided a paradigm for both understanding the pathogenesis of complex diseases such as AAA and managing the personal health of patients (Li et al., 2018). Duren et al. (2017) developed a novel method (paired expression chromatin accessibility) to construct enhancer gene network and enable us to view how the regulatory elements works. Both of these studies provided future directions for studying the pathogenesis of AAA.

The comparison of gene modules based on GO-BP and KEGG between the AAA and normal conditions indicates that these two conditions have different co-expression patterns. Using RNA microarray, Spin et al. (2011) showed that in AAA mice, immune genes were widely upregulated. Furthermore, evidence for the presence of innate and adaptive immune systems in the pathogenesis of AAA has been provided by many studies (Coscas et al., 2018; Lareyre et al., 2017; Senemaud et al., 2017; Trachet et al., 2016). Our finding that the immune response was clustered in the black module is consistent with the studies mentioned above. Inflammation also plays an important role in the pathogenesis of AAA. In addition, Spin et al. (2011) found that inflammatory genes are also up-regulated in AAA mice. Cells related to inflammation such as B cells, T cells, natural killer cells and macrophages and their products have been found in AAA tissue (Biros et al., 2014; Furusho et al., 2018; Meng et al., 2016; Tsuruda et al., 2008; Zhang et al., 2011). Studies have also shown that regulatory T cells can inhibit AAA development in human and animal models such as mice by inhibiting the inflammatory response in AAA (Ait-Oufella et al., 2013; Yodoi et al., 2015; Zhou et al., 2015), which further supports the inflammation-mediated theory of AAA. Our findings on the turquoise module agree with the results reported by these studies. Our study has also shown that vascular smooth muscle contraction is associated with AAA pathogenesis. Unlike the immune response and inflammation, little is known about the relationship between vascular smooth muscle contraction and the pathogenesis of AAA. However, Chew, Conte & Khalil (2004) and Chew, Orshal & Khalil (2003, 2004) found that elastase and matrix metalloproteinase (MMP, mainly MMP2 and MMP9) both promote early dilation of AAA by inhibiting vascular smooth muscle contraction through the mechanism of Ca²⁺ entry inhibition of vascular smooth muscle. Based on the GO-BP and KEGG pathway analyses, we found that the genes in the DEG co-expression network were mainly associated with oxidative stress and inflammation. Emeto et al. (2016) summarized that inflammation, oxidative stress and reactive oxygen species play an important role in the pathogenesis of AAA and believed that oxidative stress induces inflammation and that tissue injury is caused in the pathogenesis of AAA.

The hub genes of the black, turquoise and red modules in the AAA condition were MAP4K1, STYX and NFIB, respectively. MAP4K1, also named hematopoietic progenitor kinase 1, is expressed in T cells, B cells, and macrophages, as well as hematopoietic progenitor cells. It serves as a regulator of JNK and T and B cell signaling. It can also regulate immune cell adhesion (Chuang, Wang & Tan, 2016) and is associated with T and B cell activation and cytokine production (Alzabin et al., 2009; Shui et al., 2007), which is consistent with our findings in the functional enrichment analysis of the black module in the AAA condition. STYX is a pseudophosphatase with no actual catalytic activity. It serves as an anchor for MAPKs and plays an important role in their function (Reiterer, Pawlowski & Farhan, 2017). In addition, STYX regulates apoptosis in breast cancer cells and HeLa cells through cross-talk with FBXW7 (Reiterer et al., 2017). Therefore, STYX may also play a role in regulating vascular smooth muscle apoptosis in the pathogenesis of AAA. NFIB is a nuclear transcription factor; however, its detailed functions are still not known (Ono & Okada, 2018). It can be downregulated by microRNA-29, which is an immune system threshold. In conclusion, these three hub genes from the three key WGCNA gene modules may be of vital importance in the pathogenesis of AAA.

Yip1 domain family member 6 was the gene with the highest degree in the hub gene cluster generated from the DEG co-expression network. YIPF6 is a member of the Yip1 family of proteins and is located in the Golgi apparatus. It is believed to participate in vesicular transport (Kranjc et al., 2017). Brandl et al. (2012) discovered that one type of null mutation of YIPF6 induces spontaneous intestinal inflammation in mice, but the mechanism seems to be related to granule secretion by Paneth and goblet cells. Whether YIPF6 is associated with chronic inflammation in AAA is unknown, but it may play a role in the balance between proinflammatory and anti-inflammatory cytokine production. However, the ROC curve showed that YIPF6 as well as other genes in this cluster has potential diagnostic value.

In this study, we used WGCNA to construct a co-expression network, detect gene modules and identify hub genes in AAA for the first time. In addition, we merged and normalized three GEO datasets instead of using only one to conduct our analysis. However, due to the lack of clinical trait data in the three GEO datasets, we did not correlate gene modules with clinical traits when conducting the WGCNA.

Here, our study revealed that the co-expression patterns between the AAA and normal conditions are different. These differences indicate that immune response, inflammation and vascular smooth muscle contraction may play important roles in the pathogenesis of AAA. Then, we identified hub genes of gene modules in both the AAA and normal conditions and hub gene clusters of the DEG co-expression network. These crucial genes may have crucial biological functions in the pathogenesis of AAA.

Conclusions

Based on WGCNA, we detected six modules in the AAA condition and 10 modules in the normal condition. Hub genes of each module and hub gene clusters of the DEG co-expression network were identified. These genes may act as potential targets for medical therapy and diagnostic biomarkers. Further studies are needed to elucidate the detailed biological function of these genes in the pathogenesis of AAA.

Supplemental Information