Bioinformatics-Based Identification of Potential Hypoxia-Related Genes Associated With Peyronie’s Disease

Abstract

Hypoxia is one of the most important predisposing conditions for Peyronie’s disease (PD) and the pathogenetic mechanism is yet to be completely elucidated. This study applied bioinformatic approaches to select candidate hypoxia-related genes involved in the pathogenesis of PD. The Gene Expression Omnibus (GEO) data set GSE146500 was introduced to compare the transcriptional profiling between normal and PD samples. The differential expression of hypoxia-related gene was determined with R software. On the selected candidate genes, further functional analyses were applied, including protein–protein interactions (PPIs), gene correlation, gene ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. A total of 66 candidate genes (24 candidates overexpressed in PD and 42 showing reduced expression in PD) were distinguished according to the differential expression between human fibroblast cells from normal and PD patients. The interactions among these candidate genes were recognized according to PPI analysis. The functional enrichment analyses revealed the potential modulatory functions of the candidate genes in some major biological processes, especially in glycolysis/gluconeogenesis and carbon metabolism. The findings would facilitate further study on the pathogenesis of PD, which might consequently promote the improvement of clinical strategies against PD.

Introduction

Peyronie’s disease (PD) is a condition of erectile dysfunction (ED) in which the formation of fibrotic plaque can be resulted from the abnormal healing of the tunica albuginea (TA), which affects 5% to 9% of men (Gholami et al., 2003). Typically, the PD-inducing fibrosis of connective tissue is characterized by intracellular collagen accumulation. The upregulation of collagen production relies on the activation by profibrotic factors, including transforming growth factor (TGF) β1, plasminogen activator inhibitor–1 (PAI-1), and excessive reactive oxygen species (ROS; El-Sakka, 2011). Consequent to the progressive fibrotic process, apart from the reduced elasticity of TA directly causing ED, the major symptoms of PD also include penile deformity and intense pain, which can impede penetrative intercourse and therefore induce relationship difficulties and mental diseases. Compared with the physical discomfort, the complications of PD can ruin the life quality of males to a remarkable extent (Gaffney & Kashanian, 2020). To relieve the threats from PD to public health, the researchers have managed to partially recognize the predisposing conditions, such as inflammatory infiltration, oxidative stress, and hypoxia (Zhang et al., 2021). To date, the conventional treatments for PD have limited curative effects and therefore are considered unsatisfactory. The hope of improving the clinical strategies against PD relies on the progress in exploring its pathogenetic mechanisms, especially the molecular signaling processes and the underlying transcriptional changes.

Hypoxia inducible factor–1 (HIF-1), an oxygen-sensitive transcriptional activator, is a crucial mediator in the process of proper tissue repair (Young & Moller, 2010). Consequent to even minor tissue injury, the oxygen supply via vascular perfusion for the local cells can be instantly impeded. The cellular responses to the hypoxic state involve multiple genes, thereby activating the tissue repair processes to restore the biological functions of the damaged tissue (Chiche et al., 2009). To facilitate the regeneration of damaged tissue, appropriate stimulating conditions are required, such as the transient hypoxia caused by minor tissue injury and the subsequent HIF-1-mediated activation of the downstream pathway. Tissue regeneration can at least partially retrieve the original structural organization of the injured site and maintain the normal functions as much as possible. However, once the severity or the duration of the injury exceeds the healing power of tissue regeneration, the resultant chronic hypoxia may activate an alternative pathway that induces fibrosis, thereby resulting in scar formation and irreversible loss of normal biological functions (Tanaka, 2016).

As the clinical significances of the pathogenetic factors in PD were realized among the researchers, some key regulatory factors for the relevant signaling pathways have been investigated. Cirino et al. demonstrated that hypoxia, HIF-1, and the target genes of HIF-1 could contribute to PD development in Tsk mice (Lucattelli et al., 2008). However, the correlation of hypoxia-related genes to the development of PD had not been discussed in depth. The selection of hypoxia-related genes significant in PD development would help to establish an advance model for PD pathogenesis and suggest some potential targets for therapeutic measures.

In this study, based on Gene Expression Omnibus (GEO) data set GSE146500, the hypoxia-related genes were selected according to the altered expression levels in PD samples. Then, further functional analyses were applied, including protein–protein interactions (PPIs), gene correlation, gene ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The functional enrichment analyses revealed the potential modulatory functions of the candidate genes in some major biological processes (BPs), especially in glycolysis/gluconeogenesis and carbon metabolism. The findings would facilitate further study on the pathogenesis of PD, which might consequently promote the improvement of clinical strategies against PD.

Materials and Methods

Microarray Data and the Selection of Hypoxia-Related Genes

Through the GEO website (http://www.ncbi.nlm.nih.gov/geo/), the data set GSE146500 was obtained, which contained the comprehensive transcriptional profiles of normal and PD patient samples. The expression profiling was conducted on the GPL16791 platform, the high-throughput sequencing of which utilized Illumina HiSeq 2500 (Homo sapiens) Technology type. The analytical objects included four samples of human normal fibroblasts (undergoing penoplasty for congenital curvature) and four samples of human PD fibroblasts (the fibroblasts were taken from the plaques). Gene Set Enrichment Analysis (GSEA; http://www.gsea-msigdb.org/gsea/index.jsp) provided the list of genes (n = 200) that had confirmed association with hypoxia. As the database is available to the public by means of open access, authorization from the local ethics committee was not necessary.

Differential Expression Analysis of Hypoxia-Related Genes

Based on the normalized expression profiling provided by the data set GSE146500, the differential expression between normal and PD samples could be determined. The probe annotation was applied according to the corresponding files included in GSE1465000. The principal components analysis (PCA) was conducted to exclude batch effect. The “limma” package of R software was used to identify candidate hypoxia-related genes with differential expression levels in PD. The changes in expression levels with adjusted absolute fold-change > 1 and p < .05 were considered as significant differences. The visualizing interpretation of the results utilized the “heatmap” and “ggplot2” packages of R software.

PPI and Gene Correlation Analysis of the Candidate Genes

PPI analysis relied on the STRING database (https://string-db.org/) and Cytoscape software (version 3.8.1). The gene correlations were checked by the means of Spearman correlation in the “corrplot” package of R software.

GO and KEGG Pathway Enrichment Analysis of the Candidate Genes

GO and KEGG pathway enrichment analysis utilized the “GO plot” package of R software. The setting for GO aspects included molecular function (MF), BP, and cellular component (CC).

Statistical Analysis

R software (version 3.6.2) was utilized to acquire the statistics. Student’s t test was assigned to evaluate the significances for differences. The p value of <.05 was considered as statistically significant.

Results

Retrospective Identification of Hypoxia-Related Genes Differentially Expressed in PD

According to the PCA results, the repeatability of intragroup data in GSE146500 passed the validity check (Figure 1A). By analyzing the GSE146500 data set with R software and comparing the expression profiling of human fibroblasts from normal and PD penile TA, the genes were selected according to differential expression and presented in volcano plot (Figure 1B). Therefore, the differences in the expression patterns of 200 hypoxia-related genes were evaluated by comparing the samples from PD patients and healthy individuals. According to our selection criteria (absolute fold-change > 1, adjusted p < .05), 66 hypoxia-related genes with differential expression in PD were distinguished, including 42 genes showing reduced expression and 24 showing overexpression in PD (Figure 1C). Next, the 66 candidate genes were plotted in a heatmap (Figure 1D) generated with R software. The box plots also allowed the visualized comparison of the differential expression levels in PD among the 66 candidate genes (Figure 2A and B). The top five hypoxia-related genes overexpressed in PD were IL6, F3, CSRP2, SLC2A3, and VLDLR, and the top five hypoxia-related genes suppressed in PD were PGF, EFNA1, SELENBP1, PPARGC1A, and GPC4.

Figure 1.

Differentially Expressed Hypoxia-Related Genes in PD and Healthy Samples. (A) Principal Components Analysis for GSE146500. (B) Volcano Plot of the Differentially Expressed Genes. The Red Dots Represent the Significantly Upregulated Genes and the Blue Dots Indicate the Significantly Downregulated Genes. (C) 66 Hypoxia-Related Genes Were Identified Using the Criteria of Absolute Fold-Change Value > 1 and Adjusted p Value < .05, Including 42 Downregulated Genes and 24 Upregulated Genes. (D) Heatmap of the 66 Differentially Expressed Hypoxia-Related Genes in PD and Healthy Samples.

Figure 2.

The Boxplot of 66 Differentially Expressed Hypoxia-Related Genes in PD and Healthy Samples. (A) The Boxplot of Top 24 Differentially Expressed Hypoxia-Related Genes in PD and Healthy Samples. (B) The Boxplot of Last 42 Differentially Expressed Hypoxia-Related Genes in PD and Healthy Samples.

PPI and Correlation Analysis on the Candidate Genes

In this study, the construction of PPI network for the candidate genes also relies on the bioinformatic approach. By analyzing the data from the STRING with the Cytoscape software, the detection of PPI suggested the interaction among the candidate genes (Figure 3A). Figure 3B showed the numbers of detected interactions for the proteins encoded by the candidate genes. Gene correlation analysis directly assessed the coregulated expression of the candidate genes (Figure 4).

Figure 3.

PPI Analysis the 66 Differentially Expressed Hypoxia-Related Genes. (A) The PPI Among 66 Differentially Expressed Hypoxia-Related Genes. (B) The Interaction Number of Each Differentially Expressed Hypoxia-Related Gene.

Figure 4.

Spearman Correlation Analysis of the 66 Differentially Expressed Hypoxia-Related Genes.

Functional Enrichment Analyses on the Candidate Genes

GO and KEGG enrichment analysis could provide more details about the functions of the 66 candidate genes. In the aspect of BPs, the most intensively enriched GO terms indicated the potential functions of the transcriptional products in monosaccharide (e.g., hexose) biosynthesis, gluconeogenesis, and ADP metabolism. Also, the enrichment of GO terms indicates the protein products of the candidate genes could function as the CCs of the collagen-containing extracellular matrix, vacuolar lumen, Golgi lumen, and lysosomal lumen. In addition, the expression of the candidate could have influences on carbohydrate binding and monosaccharide binding (Figure 5A–C). According to the KEGG pathway analysis, the candidate genes mainly displayed their regulatory roles in glycolysis/gluconeogenesis, carbon metabolism, pentose phosphate pathway, and galactose metabolism (Figure 6).

Figure 5.

GO Enrichment Analysis of 66 Differentially Expressed Hypoxia-Related Genes. (A), (B), and (C) Bubble Plot of Enriched GO Terms.

Figure 6.

KEGG Analysis of 66 Differentially Expressed Hypoxia-Related Genes.

Discussion

In most cases, fibrotic diseases and various predisposing factors for tissue injuries (such as mechanical stress, infectious pathogens, and autoimmune disorders) are closely related to the sustained damage of cells and tissues. Characterized by the chronic fibrosis of TA, PD can be the direct consequence of repeated trauma associated with sexual activities or other vigorous exercises, thereby leading to penile malformation, corporovenous occlusive ED, psychological problem, and disrupted intimate relationship (Garaffa et al., 2013). To preliminarily establish the theoretic model for PD pathogenesis, two dominant hypotheses have been proposed. As the first hypothesis, the inner layer of TA and the sinusoidal tissue can be separated by minor trauma, followed by the formation of micro hematoma in the connective tissue sleeve between TA and corpus cavernosum. Accordingly, the lesion should involve the dorsomedial side of TA and the distinctive extension into the corpus cavernosum (Milenkovic et al., 2019). The second hypothesis suggests that the TA is separated between its inner and outer layers, especially in the dorsomedial side formed by the septum between the cavernous bodies, due to repeated trauma (Graziottin, 2015).

Reported by preceding experimental studies, as the result of primary injury to the penile tissue, reactive oxygen species (ROS) and PAI-1 significantly accumulated in the fibrotic plaque and thereby led to the increase of oxidative stress (Davila et al., 2003). In response to the enhanced fibrogenic factors, nitric oxide (NO) production, mediated by inducible nitric oxide synthase (iNOS), showed significant increase both in human and in animal plaques, which can inhibit ROS activity, collagen accumulation, plaque formation, and its product cGMP (Vernet et al., 2002). Previous articles investigated the effects of iNOS silencing in mice, thereby confirming its anti-fibrotic effect and highlighting its significance in protecting smooth muscle cells (SMCs) of corpus cavernosum (Ferrini et al., 2010). Depending on animal model, iNOS and hypoxia-inducible factor–1 (HIF-1) were upregulated in PD-like lesions (Krakhotkin et al., 2020). Therefore, inhibiting the expression or the product activity of hypoxia-inducible profibrotic gene may be an exciting supplement to PD treatment in the future.

The preceding urological studies had partially revealed the significances of hypoxia-related genes, especially in malignancies. However, the correlation of hypoxia-related genes to the development of PD had not been discussed. In this article, 66 candidate hypoxia-related genes differentially expressed in PD were distinguished using bioinformatic approaches. Some of our candidate genes had been explored to some extent by other research groups. For example, Atar et al. (2017) reported that the serum concentration of interleukin-6 (IL-6) showed remarkable increase in PD patients compared with the control group. The IL-6 levels detected in ED patients were obviously elevated. Zimmermann et al. (2008), analyzing IL-6 serum levels in 91 patients with stable phase of PD, reported that the IL-6 levels were significantly lower compared with those in healthy controls. The main difference between the two articles is that, because the former included only patients in the acute phase of PD, IL-6 levels were significantly higher possibly due to the ongoing local inflammation. PD is generally divided into two different phases: acute and chronic. Plaque formation generally occurs during the acute phase, whereas during chronic phase pain usually tends to complete resolution and penile deformity stabilizes. PD’s pathophysiology is still subject of great discussion (Di Maida et al., 2021). In addition, according to valid evidence, after bilateral cavernous nerve injury (BCNI), increased expression and activity of penile lysyl oxidase (LOX) could be detected, which promoted penile fibrosis, reduced SMC content, and induced ED (Wan et al., 2018). Besides, Mateus et al. (2018) demonstrated the effect of an adenosine receptor A2B (ADORA2B) agonist on TGF-β1-induced myofibroblast transformation.

According to the GO enrichment analysis, the potential biological functions of the candidate genes were mainly displayed in the following aspects: ADP metabolic process, monosaccharide biosynthetic process, hexose biosynthetic process, and gluconeogenesis (BP); collagen-containing extracellular matrix, vacuolar lumen, Golgi lumen, and lysosomal lumen (CC); and carbohydrate binding and monosaccharide binding (MF). Suggested by the results of KEGG enrichment analysis, the differentially expressed hypoxia-related genes are mainly involved in the process of glycolysis/gluconeogenesis, carbon metabolism, pentose phosphate pathway, and galactose metabolism. Based on two published articles, hypoxia could affect the progress of PD. One of the studies proved that hypoxia, HIF-1, and HIF-1 target genes could notably contribute to PD development in Tsk mice (Lucattelli et al., 2008). Moreland and his colleagues (1995) reported that TGF-β1 promoted collagen accumulation in human corpus cavernous SMCs in response to hypoxia. In addition, hypoxia could enhance TGF-β1 expression and suppress prostaglandin E (PGE). Accordingly, the increased PGE1 level could inhibit TGF-β1-induced collagen synthesis (Moreland et al., 1995). Future experiments are needed to verify the potential biological functions of our candidate hypoxia-related genes in PD pathogenesis.

Among the known biomolecules modulating the signaling pathways relevant to hypoxia responses, the transcription factor HIF-1 is noticeable. It regulates the transcriptional activity of multiple essential genes in tissue repair and the target genes may have pro-inflammatory, angiogenic, and potentially profibrotic properties (Wang et al., 1995). Consequent to the enhanced expression of these genes, the cells could adapt to the hypoxic environment via enhanced cell proliferation and migration, accelerated glucose metabolism, the promotion of angiogenesis, and increased cellular survival. Our KEGG enrichment analysis demonstrated that HIF-1 signaling pathway was also one of the most significant processes in which the differentially expressed hypoxia-related genes were mainly involved.

There are still some limitations to our study. First, as the bioinformatic results were obtained from four samples of human fibroblasts from normal and PD penile TA, respectively, the number of clinical samples is limited. Besides, the included data set failed to describe the phases of PD. Whether men undergone any prior injection therapy for PD prior to the collection of fibroblasts was also not mentioned. As no experiment was conducted, this is a major limitation of this study. Second, because most of our candidate genes had not been reported to be PD-related, it was definitely required to conduct in-depth investigation on the PD-related functions of these genes. Setting up mouse or other animal models, in which these genes can be conditionally knocked out, may facilitate the efforts to elucidate their regulatory roles in PD development. Therefore, further research needs to be explored in the future.

Conclusion

Taken together, based on the bioinformatic approaches, the selection of candidate hypoxia-related genes further clarified the pathogenic mechanism underlying PD. The hypoxia-related genes significant in PD development would help to establish an advance model for PD pathogenesis and suggest some potential targets for therapeutic measures. However, the thorough evaluation of clinical potentials for the candidate genes required in-depth investigation on their regulatory functions and the corresponding signaling pathways.