Screening and identification of critical biomarkers in erectile dysfunction: evidence from bioinformatic analysis

Purpose Erectile dysfunction (ED) is one of the most common male-disease globally. Despite efforts to explain its pathogenesis, the molecular mechanisms of ED are still not well understood. Methods The microarray dataset GSE10804 was downloaded from the Gene Expression Omnibus (GEO) to find candidate genes in ED progression. After differentially expressed genes (DEGs) were identified, functional enrichment analysis was performed. In addition, a protein-protein interaction network (PPI) was established and module analysis was performed through the STRING and Cytoscape. Results and Conclusions A total of 618 DEGs were identified in all, containing 430 downregulated genes and 188 upregulated genes. The enriched functions and pathways of the DEGs include transcription from RNA polymerase II promoter, cell adhesion, calcium ion binding, receptor binding, Akt signaling pathway, receptor interaction, protein digestion, and absorption. We picked out twenty-five hub genes, with biological process (BP) analyses revealing that the genes were principally associated with cellular responses to amino acid stimuli, extracellular matrix structural constituent, collagen trimer, protein digestion and absorption, ECM-receptor interaction and PI3K-Akt signaling pathway. To sum up, DEGs and hub genes distinguished in this study not only help us understand the molecular mechanisms behind the carcinogenesis and progression of ED, but also play a part in the diagnosis and treatment of ED by providing candidate targets.


INTRODUCTION
Erectile dysfunction is a medical condition in which a man is unable to maintain an erection for vaginal intercourse (1993). Actually, a description of poor penile erection was firstly mentioned as early as 5,000 years ago in an ancient Egyptian literature (Shamloul & Ghanem, 2013). Erectile dysfunction is known to increase with age. In a cross-sectional study based on community, the prevalence of this condition, in men in the 40-49 age group, was about 5% for severe dysfunction (Feldman et al., 1994). For men in the 70-79 age group, the prevalence rates were 15% and 34%, for complete or severe dysfunction, respectively. By 2025, it is speculated that the global prevalence of ED will reach up to 322 million cases (Bacon et al., 2003). Men aged older than 40 years are more susceptible to erectile dysfunction. Several other medical diseases are linked to erectile dysfunction such as diabetes mellitus (Brown et al., 2005), tobacco use, central neuropathologic conditions, cardiovascular disease (Billups, 2005), lower urinary tract symptoms of benign prostatic hyperplasia (McVary, 2005), and metabolic syndrome (Kupelian et al., 2006).
The vessels of the vascular system (and lymphatics) contain a single-cell-layered lining called the endothelium. Based on its uniqueness and functions, the endothelium is thought to be an organ (Aird, 2004). The parts of the endothelium are diverse regarding biochemical properties, functional characteristics, and cell surface markers (Ho et al., 2003;Kallmann et al., 2002). The understanding of the heterogeneity of this factor in the endothelial cells (ECs) mainly originates from several examinations of various body tissues and organ specimen by microarray analysis (Chi et al., 2003;Ho et al., 2003;Kallmann et al., 2002;Podgrabinska et al., 2002).
Endothelial dysfunction is considered as the genesis of ED, especially that of the systemic vasculature (Gerber et al., 2015). Plenty of studies have been carried out to investigate the function of the endothelium both during physiological erection and also in the ED state (Hurt et al., 2002). The endothelium not only functions as a passive barrier but also as a critical regulator of muscle tone and vascular circulation, which are coordinated by signals originating from mechanical stimuli, humoral and neural sources.
In order to facilitate research into the progression of ED in terms of the functional pathways and differentially expressed genes (DEGs), many technologies such as bioinformatics analysis and microarray technology have been widely used. However, because it is exceedingly difficult to obtain penile tissue specimens from patients with ED, there are few studies of human tissue. Thus, here, the DEGs between non-ED tissues and ED tissues were searched and obtained by screening the Gene Expression Omnibus (GEO) for mRNA microarray datasets. Subsequently, protein-protein interaction (PPI) networks, pathway enrichment analyses of the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO), were performed to explore the mechanisms that underlie the pathogenesis of ED. In conclusion, we identified 618 DEGs and 25 hub genes as possible biological markers for ED. GEO (http://www.ncbi.nlm.nih.gov/geo) (Barrett et al., 2013;Edgar, Domrachev & Lash, 2002), a public functional genomics data repository, provides chips, microarrays, and high throughput gene expression data. Three gene expression datasets GSE10804 (Wessells et al., 2009) were downloaded from GEO (Affymetrix GPL571platform, Affymetrix Human Genome U133A 2.0 Array). The probes were transformed into the corresponding gene symbols on the basis of annotation information on the platform. The GSE10804 datasets had 3 groups: 5 samples are cultured human corpus cavernosum endothelial cells (HCCEC) from a donor with erectile dysfunction; 3 samples are cultured human umbilical vein endothelial cells (HUVEC) from a donor without erectile dysfunction, and 4 samples are cultured human coronary artery endothelial cells (HCAEC) from a donor without erectile dysfunction. We have standardized that the 5 ED samples compared to the 7 non-ED samples as a whole.

Identification of DEGs
The DEGs in the ED and non-ED specimen were picked out from GEO2R (http: //www.ncbi.nlm.nih.gov/geo/geo2r), which is a platform for examing DEGs across experimental conditions by comparing multiple datasets in a GEO series. The genes with multiple probes were averaged, and the probes that lacked gene symbols were removed. P-value <0.01 and logFC (fold change) >1 represented statistical significance.

Analysis of GO enrichment and KEGG networks for DEGs
The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov) (version 6.8) (Huang da, Sherman & Lempicki, 2009) offers an analytical platform with a comprehensive source of annotated information of proteins and genes which can be extracted and analyzed. On the other hand, KEGG (Kyoto Encyclopedia of Genes and Genomes) is a capable platform that is a knowledge base for systematic analysis of gene functions, linking genomic information with higher-order functional information (Kanehisa & Goto, 2000;Kanehisa et al., 2019;Kanehisa, 2019). GO, a significant bioinformatics tool enables us to annotate genes according to biological processes (Gene Ontology, 2015). In order to study the functions of DEGs, DAVID online database was used to analyze the biological. P < 0.01 was taken to represent statistical significance. ImageGP (http://www.ehbio.com/ImageGP/) is used to draw Enrichment Plot for GO and KEGG.

Construction of PPI networks and module analyses
In this study, the STRING database (STRING; http://string-db.org) (version 10.0) (Szklarczyk et al., 2017) was taken to predict the PPI network and to determine proteinprotein interactions and to investigate the molecular basis of diseases. The interactions with an average score of greater than 0.4 were chosen a the cut-off for statistical significance. To visualize molecular interaction networks, the Cytoscape (version 3.4.0), a bioinformatics software, is often utilized (Saito et al., 2012). This software contains an APP known as The plug-in Molecular Complex Detection (MCODE) (version 1.4.2), which is used to cluster-specific networks according to the on the topology to reveal the regions that are densely connected (Saito et al., 2012). Cytoscape was used to draw the PPI networks, and the MCODE helped to identify the most remarkable module in the PPI networks. The selecting criteria was: k-score = 2, Max depth = 100, node score cut-off = 0.2, degree cut-off = 2, and MCODE scores >5. Then, the DAVID tool was applied to analyze the GO and KEGG genes in the module.

Selection and analysis of Hub genes
Hub genes with a degree greater than or equal to 10 were selected. Analysis and visualization of the hub genes were conducted by the Biological Networks Gene Oncology tool (BiNGO) (version 3.0.3) plugin of Cytoscape (Saito et al., 2012).

Identification of DEGs in the ED
DEGs (618 in GSE10804) were identified ( Fig. 1A) after standardizing the microarray data. The dataset consisted of 430 downregulated genes and 188 upregulated genes between ED endothelial cells and non-ED endothelial cells.

KEGG and GO enrichment analyses of DEGs
DAVID tool is used for the analysis of DEGs classification, functional and pathway enrichment. Results of GO analyses revealed that changes in BP of DEGs had great enrichment in increasing gene expression from RNA polymerase II promoter, cell adhesion, extracellular matrix organization, cell-cell signaling and positive control of cell migration (Fig. 1B). Molecular function (MF) mainly had changes in calcium-binding, receptor binding, sequence-specific DNA binding, actin binding and receptor activity (Fig. 1C). Cell component (CC) of DEGs primarily had changes in the cell surface, proteinaceous extracellular matrix, focal adhesion and plasma membrane (Fig. 1D). KEGG pathway analyses showed that DEGs mainly had enrichment in Focal adhesion, the Akt signaling pathway, receptor interaction, Protein digestion and absorption, and Pathways in cancer (Fig. 1E). Functional and pathway enrichment analyses were conducted on downregulated and upregulated DEGs respectively for further analyses of the biological classification. In the upregulated DEGs, the Results of GO analyses demonstrated that the gamma-aminobutyric acid signaling pathway and osteoclast differentiation had changes in BP (Table 1). KEGG analysis showed the main enrichment in the Renin secretion of downregulated DEGs (Table 1). MF and CC showed no results in the upregulated DEGs. In the downregulated DEGs, the Results of GO analyses showed that the remarkable enhancement of changes in BP was in the collagen catabolic process, extracellular matrix organization, and cell adhesion (Table 2). Changes in MF were greatly improved in calcium ion binding, platelet-derived growth factor binding, and heparin-binding (Table 2). Changes in CC primarily emerged in the extracellular matrix, extracellular region and extracellular matrix (Table 2). KEGG pathway analysis revealed that the downregulated DEGs were mainly enhanced in the PI3K-Akt signaling pathway, ECM-receptor interaction, and Focal adhesion ( Table 2).

Construction of PPI networks and module analyses
The constructed PPI networks of DEGs are shown in (Fig. 2A), and the most significant module was obtained using Cytoscape (Fig. 2B). The DAVID tools were applied in the functional analyses of genes associated with this module. We found that the genes primarily had enrichment in cell responses to amino acid stimuli, extracellular matrix structural constituent and collagen trimer (Table 3).

DISCUSSION
Men aged older than 40 years are more susceptible to erectile dysfunction. Several other medical diseases are linked to erectile dysfunction such as diabetes mellitus (Brown et al., 2005), tobacco use, central neuropathologic conditions (e.g., hemorrhagic or ischemic    et al., 2006). More than 70% of coronary artery disease patients present with moderate to severe erectile dysfunction which might occurred before cardiovascular events. Moreover, endothelial dysfunction is also linked with atherosclerosis (Billups, 2005). The likelihood of developing erectile dysfunction tends to increase with prolonged diabetes probably due to a high level of glycated hemoglobin. Some drugs are also suspected to cause erectile dysfunction based on a clinical diagnosis of 25% of the patients diagnosed with erectile dysfunction (Thomas et al., 2003). However, little is understood about the molecular mechanism of ED at present. The NO-cGMP signaling pathway has been shown to be critical in the pathogenesis of ED (Garcia-Cardoso et al., 2010). These results have also been applied in clinical practice, as selective inhibitors of 5-type phosphodiesterase specific to cGMP, such as sildenafil and tadalafil have achieved significant clinical effects in the treatment of ED. RhoA/Rho-kinase contributes to the development and progression of ED. It has been demonstrated that the vasodilation of NO may be due to the reduction of calcium levels by NO-cGMP-PKG on one hand, and the relaxation effect by blocking the RhoA/Rho-kinase calcium-sensitive pathway on the other hand (Jin et al., 2006). Recently, it has also been found that the MAPK signaling pathway was demonstrated to affect the NO-cGMP signaling pathway by regulating related kinases, which indirectly leads to the occurrence of ED. The primary treatments for erectile dysfunction are oral PDE5-Is. Other interventions are psychotherapy, penile devices, testosterone therapy, injection therapies, and psychotherapy, etc (1993). But there is no curable medications for ED that can be completely cured, which may be one of the reasons for the poor patient prognosis. Therefore, potential markers for diagnosis and treatment with high efficiency are urgently demanded. Microarray technology allows us to explore the genetic alterations in the ED and there is evidence surporting it to be a useful method for identifying new biomarkers in other diseases. An article conducted a large-scale genome-wide association study of erectile dysfunction in the multiethnic Kaiser Permanente Northern California Genetic Epidemiology Research in Adult Health and Aging cohort and UK Biobank. They find the single-minded family basic helix-loop-helix transcription factor 1 (SIM1) gene is a mechanism that is specific to sexual function (Jorgenson et al., 2018).
Here, we analyzed the mRNA microarray datasets to identify DEGs in ED cells and non-ED cells. 618 DEGs were obtained in all datasets, including 430 downregulated genes and 188 upregulated genes. The reciprocities among the DEGs were searched out by KEGG enrichment analyses. The DEGs largely had enrichment in positively regulating transcription from RNA polymerase II promoter, cell adhesion, calcium ion binding, receptor binding, Akt signaling pathway, receptor interaction, Protein digestion, and absorption. Lianjun Pan and Jiehua Ma speculated that many lncRNAs are involved in the regulation of muscle contraction and relaxation through ion channel activity. They suggested that a muscle contraction disorder induced by abnormal ion channel activity, collagen deposition and fibrosis might play a critical role in the pathogenesis of ED (Pan et al., 2015). Changes to the endothelial cell-cell junction integrity like in the case of increases in vascular permeability and decreased cell-cell junction proteins were observed in the cavernous tissue of the diabetic mouse. The presence of hollow junctions in diabetes state may explain the high prevalence of erectile dysfunction that occurs before other cardiovascular conditions (Ryu et al., 2013). Large conductance, voltage, and Ca2+ activated K+ channels are widely present in many cells of the body, which play a part in regulating the neuronal excitability and smooth muscle tone. When their function is compromised, they may bring about the development of erectile dysfunction, hypertension and overactive bladder (Carmona et al., 2016). Akt plays an essential role in regulating cell homeostasis through its effects on several downstream effectors (Carmona et al., 2016).
Prior studies have shown that phosphorylation of eNOS by Akt causes prolonged NO production and maximal construction and influences penile erection (Zhang et al., 2011). Similarly, it was found that PI3K/Akt-eNOS activation mediated the effects of the A2B adenosine receptor on penile erection (Ryu et al., 2013). We selected 25 DEGs as hub genes with difference degrees ≥10. Among these hub genes, COL3A1 showed the highest node degrees with 13. And other hub genes also confirmed the top node degrees with 12 (COL1A1, COL1A2, COL21A1, COL2A1, COL4A1, COL4A2, COL4A4, COL5A3, COL6A1, COL6A2, COL6A3, CXCL12, LEPREL1) and 11 (C3, CCR5, CXCL11, CXCL5, CXCL6, DRD2, FPR2, GABBR2, GNAI1, GNG7, LPAR1). A study identification of oxidative stress-induced gene expression profiles in cavernosal endothelial cells, cxcl12 showed a corresponding variation in the CECs and ED rat model compared with the results of the gene microarray analysis (Hu et al., 2015). While ORP in NHPs produced persistent erectile and urinary tract dysfunction. Periurethral injection of CXCL-12 was feasible and improved both urinary incontinence and erectile dysfunction and suggests (Zambon et al., 2018). A study about five immune agents (IgG, IgM, IgA, C4, and C3) was collected on the basis of the Fangchenggang Area Male Health and Examination Survey (FAMHES)and methods of traditional cross-sectional analysis were used. However, there was no clear relationship between ED risk in the baseline analyses and the tested indexes (C3: P = 0.737) (Chen et al., 2014). Another study on sexual dysfunction in male schizophrenia showed that correlation between two polymorphisms in the genes for the DRD2 and scores on questionnaires for erectile and sexual dysfunction were studied. The functional −141Ins/Del promoter region polymorphism of DRD2 was remarkably correlated with sexual dysfunction (Zhang et al., 2011). The remaining genes did not retrieve research literature related to erectile dysfunction.

CONCLUSIONS
In conclusion, 618 DEGs and 25 hub genes with the possibility to be considered as diagnostic biomarkers for ED were identified However, additional investigations should be carried out to explore the biological function of these genes in the ED.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by grant (No. 81571433) from the The National Natural Science Foundation of China and (No. 2017A030313453) from the the Science and Technology Projects of Guangdong Province. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.