Development and Validation of Hub Genes for Adrenal Aldosterone-Producing Adenoma by Integrated Bioinformatics Analysis

Background

Primary hyperaldosteronism (PA) is characterized by spontaneous secretion of excessive aldosterone and inhibition of plasma renin activity.¹ The pathogenesis of adrenal aldosterone-producing adenoma (APA) involves the abnormal proliferation of adrenal cortex cells and the excessive secretion of aldosterone, accounting for nearly 30% of PA. Excessive secretion of aldosterone can lead to retention of sodium, suppression of plasma renin levels, increased excretion of potassium, secondary hypertension (HTN) and end-organ damage.^2,3 APA is known as a heterogeneous disease. The occurrence, development and metastasis of tumors are usually attributed to the combined effects of genetic aberrations, cell environment changes and internal and external environmental influences. Therefore, clarifying the cause, molecular mechanisms and key pathways of carcinogenesis of APA is of great significance for improving the diagnosis and treatment of APA.

At present, there are few studies focusing on comprehensive analysis of changes in mRNA expression involving the development of APA and PA. Previous studies have discovered several new genes in APA, but there are some limitations that cannot be ignored. Recently, weighted correlation network analysis (WGCNA) is considered to be an effective method for processing gene expression data and discovering new key genes linked to the clinicopathological characteristics of interest.^4,5 According to the pairwise correlation, highly co-expressed genes can constitute multiple gene modules. Then, the genes enriched in the identified modules are screened for further analysis. Therefore, WGCNA has been widely used to determine indicators in various cancers, including glioblastoma,⁴ bladder cancer,⁶ breast cancer,⁷ and adrenocortical carcinoma.⁸ However, there have been no previous studies using WGCNA to study the underlying mechanism of APA and PA.

In the present study, we downloaded the gene expression profiles of GSE60042 and GSE8514 from the Gene Expression Collection (GEO) dataset.^9,10 All APA samples were taken from patients with adrenal adenoma. We compared the gene expression profiles of APA samples and normal adrenal gland samples to identify differentially expressed genes (DEGs). WGCNA was performed to screen key modules and genes. Further bioinformatics analysis was used to clarify potential annotations and pathways of APA and PA. Protein–protein interaction (PPI) network module analysis was also constructed. Finally, the gene expression profiles of GSE10927 and GSE33371 were used for further validation.^11,12

Materials and Methods

Microarray Data

A total of four datasets of gene expression profiles related to aldosterone-producing adenoma (including GSE60042, GSE8514, GSE10927 and GSE33371) were downloaded from Gene Expression Omnibus (GEO) datasets (https://www.ncbi.nlm.nih.gov/geo/). There were 7 samples of APA and 7 samples of adjacent adrenal tissue in GSE60042, 10 samples of APA and 5 samples of normal adrenal tissue in GSE8514. R packages “limma” was used to identify DEGs between APA and normal adrenal tissue in GSE60042 and GSE8514. We set false discovery rate (FDR) < 0.05 and |log fold change (FC)| ≥ 1 as the cut-off criterion.

Weighted Gene Co-Expression Network Analysis (WGCNA) of DEGs

“WGCNA” package in R software was used to perform WGCNA. Firstly, data pre-processing and quality assessment were performed, and the Pearson’s correlation matrices were constructed for all paired genes. Furthermore, a weighted adjacency matrix was established utilizing a power function a_mn =|c_mn|^β (c_mn = Pearson’s correlation between gene m and gene n; a_mn = adjacency between gene m and gene n). The parameter β emphasized a strong correlation between genes and penalized a weak correlation. We selected an applicable β value to increase the similarity matrix and achieve a scale-free co-expression network. Next, we converted the adjacency matrix into a topological overlap matrix (TOM). Gene significance (GS) was defined as the log10 transformation of the P-value (GS = lgP) in the linear regression between gene expression and the clinical features. Module significance (MS) was defined as the average GS within the module. Module eigengene (ME) was considered to be the major component in the principal component analysis for each gene module. We explored the relationship between clinical features and ME to filtrate modules. Finally, the genes in module highly related to aldosterone-producing adenoma were selected for subsequent analysis.

Protein–Protein Interaction (PPI) Network Construction

Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, https://string-db.org/cgi/input.pl) database was used to evaluate the protein–protein interaction (PPI) information. The genes were mapped in STRING database, and interactive relationships were then produced and downloaded. PPI networks were then visualized by using the Cytoscape software. The plug-in “CentiScape” was used to calculate centrality parameters for screen the most important nodes. Degree score > 5 was set as the criteria. The plug-in Molecular Complex Detection (MCODE) was used to screen the top key modules in PPI network.

Functional Enrichment Analyses

Functional enrichment analyses, including Gene Ontology (GO) analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were performed using R package “clusterProfiler”, “org.Hs.eg.db”, “enrichplot” and “ggplot2”. The results of GO analysis included biological processes (BP), molecular functions (MF), and cellular component (CC).

Validation of Identified Intersection Genes

In order to further verify the results, two independent datasets, including GSE10927 and GSE33371, were used as an external validation. There were 6 samples of APA and 10 samples of adjacent adrenal tissue in GSE10927, 8 samples of APA and 10 samples of normal adrenal tissue in GSE33371. Wilcoxon test and R packages “limma” was used to compare the expression levels of these seven hub genes between APA tissues and normal adrenal tissues. P < 0.05 was considered statistically significant.

Results

Identification of DEGs in GSE60042

By comparing APA samples and adjacent adrenal tissue in GSE60042, a total of 892 DEGs were identified for the subsequent analysis, consisting of 569 significantly down-regulated DEGs and 323 significantly upregulated DEGs. The volcano plot and heatmap of these DEGs are demonstrated in Figure 1. All the 892 DEGs were extracted for WGCNA analysis.

Figure 1 (A) Heatmap and (B) Volcano plot of differentially expressed genes from GSE60042. Red represents upregulated differentially expressed genes; green represents downregulated differentially expressed genes.

Weighted Gene Co-Expression Network Analysis DEGs from GSE60042

WGCNA was applied to identify stable co-expression modules related to APA. Cluster analysis and dynamic tree cut algorithm were performed to generate six gene modules, including blue, brown, green, grey, turquoise and yellow. Each color represented one gene module. Each gene module consisted of a number of similar gene expression patterns. Finally, we found that the blue module has the most positive correlation with APA, while the turquoise module has the most negative correlation with APA. A total of 558 genes (including 255 genes in blue module and 303 genes in turquoise module) were considered to be most relevant to APA in Figure 2.

Figure 2 The results of WGCNA. (A) The calculation diagram of the weight parameter (power) of the adjacency matrix. The X axis represents the weight parameter (power). The Y axis represents the square of the correlation index from log(K) and log(P(K)). (B) WGCNA reveals clustering and modular screening based on gene expression patterns. The top is the gene tree diagram, and the bottom is the gene modules in different colors. (C) Dendrogram of characteristic genes of consensus modules obtained by WGCNA. The red line is the merging threshold. Feature genomes below the threshold represent the modules whose expression profiles should be merged due to their similarity. Clarified the overlap between the modules. Each row and each column corresponds to a module. The colors in the table indicate the gene counts at the intersection of the corresponding modules. (D) The correlation between MEblue membership and gene significance. (E) Correlation between MEturquise membership and gene significance.

Next, we mapped these 558 genes identified from WGCNA into STRING database to construct PPI network. The PPI network consisted of 396 nodes and 923 edges. Using plug-in “CentiScape” with degree of connectivity >5 as selection criteria, we eventually identify a total of 134 genes of great importance in APA development.

Functional Enrichment Analyses

GO term enrichment analysis results varied from GO classification and expression change of DEGs. As to biological process, the DEGs significantly enriched in cellular divalent inorganic cation homeostasis, calcium ion homeostasis, heart contraction and heart process. For cellular component, the DEGs significantly enriched in collagen−containing extracellular matrix, transport vesicle and transport vesicle membrane. About molecular function, the DEGs significantly enriched in metal ion transmembrane transporter activity, monovalent inorganic cation transmembrane transporter activity, channel activity, passive transmembrane transporter activity, substrate−specific channel activity and cation channel activity. More detailed GO enrichment analysis results are shown in Figure 3A. KEGG pathway analysis found five significantly enriched pathways. The DEGs significantly enriched in neuroactive ligand−receptor interaction, calcium signaling pathway, adrenergic signaling in cardiomyocytes, phagosome, hematopoietic cell lineage. More detailed KEGG pathway analysis results are shown in Figure 3B. These significantly enriched pathways and terms could help us a lot to further understand the role which DEGs played in APA occurrence and progress.

Figure 3 (A) GO and (B) KEGG enrichment analysis of 558 genes (including 255 genes in blue module and 303 genes in turquoise module) identified from WGCNA.

Identification of DEGs in GSE8514

By comparing APA samples and adjacent adrenal tissue in GSE8514, we screened a total of 1167 DEGs including 559 down-regulated DEGs and 608 upregulated DEGs. The volcano plot and heatmap of these DEGs in GSE8514 are demonstrated in Figure 4. All the 1167 DEGs were extracted for further PPI network construction.

Figure 4 (A) Heatmap and (B) Volcano plot of differentially expressed genes from GSE8514. Red represents upregulated differentially expressed genes; green represents downregulated differentially expressed genes.

Protein–Protein Interaction Network Analysis of DEGs from GSE8514

The relationship between protein–protein interactions (PPI) was downloaded from the STRING database, consisting of 720 nodes and 2781 edges. Using plug-in MCODE, we identify the two most significant module from the PPI networks. The first module has 17 hub genes, including GPER1, TAS2R5, ACKR3, S1PR4, CCL4, NPY, HCAR1, GRM3, HTR5A, PENK, SST, CORT, CHRM4, CXCL12, CCL21, TAS2R10, C3 in Figure 5A. The second module has 37 hub genes, including SOCS3, CSF2, CCND1, JUN, PTPRC, IL1R1, IL10, PTGS2, FOS, NOD2, CD19, CD40LG, PPARG, CCL2, LEP, ADIPOQ, IL6, WFS1, IGFBP1, SCG3, FGA, SCG2, FGF23, FGG, CHRDL1, CHGB, CCKBR, AVPR1A, PTGFR, TAC1, PLCB4, CCK, NTS, NTSR2, TACR2, LPAR4, NMBR in Figure 5B.

Figure 5 Two most significant module in the protein–protein interaction (PPI) network analysis (A) First module in the PPI network. (B) Second module in the PPI network. (C) Venn diagram used to identify cross genes.

Identification of Intersection Genes

Finally, we identify 134 genes from GSE60042 and 54 genes from GSE8514 as the candidate genes for APA. The Venn diagram showed that there were 11 intersection genes, including CCL21, CXCL12, C3, GRM3, AVPR1A, IGFBP1, WFS1, PTGFR, NTSR2, CCK, JUN in Figure 5C. These 11 genes were considered as the hub genes involved in the development and progression of primary aldosteronism (PA) and aldosterone-producing adenoma.

Validation of Hub Genes in Two Independent Cohort

Moreover, we utilized two independent datasets (GSE10927 and GSE33371) to validate the differential expression of these intersection genes between APA tissues and normal adrenal tissues. Eventually, a total of seven genes were successfully verified, including C3, GRM3, AVPR1A, WFS1, PTGFR, NTSR2, JUN. The expression levels of these seven hub genes between APA tissues and normal adrenal tissues from GSE10927 and GSE33371 were demonstrated in Figures 6 and 7, respectively.

Figure 6 Validation of the expression levels of seven hub genes between APA tissues and normal adrenal tissues in GSE10927. (A) AVPR1A. (B) C3. (C) GRM3. (D) JUN. (E) NTSR2. (F) PTGFR. (G) WFS1.

Figure 7 Validation of the expression levels of seven hub genes between APA tissues and normal adrenal tissues in GSE33371. (A) AVPR1A. (B) C3. (C) GRM3. (D) JUN. (E) NTSR2. (F) PTGFR. (G) WFS1.

Discussion

PA is the most common cause of secondary hypertension and accounts for approximately 20% of resistant hypertension.¹³ It has been known to us that bilateral adrenal hyperplasia (BAH) and APA are considered as the most common causes of PA. According to the latest guidelines, unilateral adrenalectomy (UA) not only is the preferred treatment for APA but also has been shown to be effective in preventing adverse consequences of PA. Nevertheless, the assessment of postoperative benefit of UA is still difficult, and the investigation of prognostic factors is processing continuously.¹⁴

In this study, comprehensive bioinformatics approaches were used to identify genetic variance between APA and normal adrenal gland. The gene co-expression network was generated; PPI network analyses were performed, and key hub genes were highlighted. Interestingly, we successfully filtrated several genes, which might be associated with the development and progression of APA. GO term enrichment analysis revealed that, during the development of APA, cellular divalent inorganic cation homeostasis, calcium ion homeostasis, collagen−containing extracellular matrix, transport vesicle, metal ion transmembrane transporter activity and monovalent inorganic cation transmembrane transporter activity played essential roles. KEGG pathway analysis revealed that neuroactive ligand−receptor interaction, calcium signaling pathway, adrenergic signaling in cardiomyocytes, phagosome, and hematopoietic cell lineage might be crucial in APA. The alterations in calcium ion homeostasis, an important component of tumor microenvironment, is considered as hallmarks of carcinoma.¹⁵ These hallmarks include evasion of apoptosis, invasion and metastasis, unlimited replication, insensitivity to antigrowth signals and sustained angiogenesis.^16–18

From GSE60042 and GSE8514, we identified 11 candidate genes. Further validation in two independent datasets revealed that there are seven hub genes involving in the development and progression of PA and APA, including C3, GRM3, AVPR1A, WFS1, PTGFR, NTSR2, and JUN. Felizola et al¹⁹ reported that the level of GRM3 mRNA was higher in APA than that in cortisol-producing adenoma (CPA) or normal adrenal cortex (NAC) and considered that GRM3 regulates steroidogenesis in adrenocortical tissues. Kiyotaka et al²⁰ revealed that, in microarray, GRM3 was highly expressed in APA; however, in qPCR analysis, there are no different expressions in APA. The difference between the two methods requires further research to confirm. Zhao et al^21,22 revealed that AVPR1A played important roles in CRPC progression and metastasis. Akiyama et al²³ suggested that PTGFR is a novel marker for tumor endothelial cells (TECs) in renal cell carcinoma. Stockwin et al²⁴ demonstrated that NTSR2 was involved in the growth and invasion in alveolar soft-part sarcoma (ASPS). Ye et al²⁵ demonstrated that C3 is correlated to gastric cancer prognosis. JUN is a well-known oncogene and has been reported to be involved in various tumors, such as breast cancer,^26,27 renal cell carcinoma,²⁸ and ovarian cancer.²⁹ The relationship between WFS1 expression and cancer development has not been fully illuminated yet. Recent studies have revealed that WFS1 might have associations with diabetes.^30,31 Metabolic reprogramming has gradually reported to be a potential and vital mechanism for cancer^32,33 and diabetes could increase cancer progression.³⁴ But little is aware of how high glucose may intrigue the genetic alterations resulting in tumor phenotypes. Hence, WFS1 might be a hub gene for APA metabolic reprogramming. In conclusion of all the results above, they indicated that above mentioned 7 hub genes might be potential biomarkers in the development of APA and PA. However, we have to admit the limitations of this study. Firstly, our sample size may be too small to observe some connections. More study with large volume of candidates is necessary. Secondly, experimental evidences are still missing. It can be more persuasive if we further validate our findings through IHC.

Strengths and Limitations of This Study

We comprehensively analyze hub genes through bioinformatics method.

These hub genes, which may involve in the development and progression of APA and PA, possibly act as potential biomarkers and novel therapeutic targets for APA and PA.

Further studies and experiments are still required to provide solid foundation of the function of the identified genes and pathways in vitro and in vivo.

Conclusions

To summarize, we comprehensively analyze hub genes through bioinformatics method. These hub genes, which may involve in the development and progression of APA and PA, possibly act as potential biomarkers and novel therapeutic targets for APA and PA. Nevertheless, to be cautious, further studies and experiments are still required to provide solid foundation of the function of the identified genes and pathways in vitro and in vivo.