Screening and identification of NOTCH1, CDKN2A, and NOS3 as differentially expressed autophagy-related genes in erectile dysfunction

Background Loss of function of key autophagy genes are associated with a variety of diseases. However specific role of autophagy-related genes in erectile dysfunction ED remains unclear. This study explores the autophagy-related differentially expressed genes (ARGs) profiles and related molecular mechanisms in Corpus Cavernosum endothelial dysfunction, which is a leading cause of ED. Methods The Gene Expression Omnibus (GEO) database was used to identify the key genes and pathways. Differentially expressed genes (DEGs) were mined using the limma package in R language. Next, ARGs were obtained by matching DEGs and autophagy-related genes from GeneCard using Venn diagrams. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of ARGs were described using clusterProfiler and org.Hs.eg.db in R. Moreover, hub ARGs were screened out through protein-protein interaction (PPI), gene-microRNAs, and gene-transcription factors (TFs) networks then visualized using Cytoscape. Of note, the rat model of diabetic ED was established to validate some hub ARGs with qRT-PCR and Western blots. Results Twenty ARGs were identified from four ED samples and eight non-ED samples. GO analysis revealed that molecular functions (MF) of upregulated ARGs were mainly enriched in nuclear receptor activity. Also, MF of downregulated ARGs were mainly enriched in oxidoreductase activity, acting on NAD(P)H and heme proteins as acceptors. Moreover, six hub ARGs were identified by setting high degrees in the network. Additionally, hsa-mir-24-3p and hsa-mir-335-5p might play a central role in several ARGs regulation, and the transcription factors-hub genes network was centered with 13 ARGs. The experimental results further showed that the expression of Notch1, NOS3, and CDKN2A in the diabetic ED group was downregulated compared to the control. Conclusions Our study deepens the autophagy-related mechanistic understanding of endothelial dysfunction of ED. NOTCH1, CDKN2A, and NOS3 are involved in the regulation of endothelial dysfunction and may be potential therapeutic targets for ED by modulating autophagy.


INTRODUCTION
Erectile dysfunction (ED), also known as impotence, is the inability to maintain enough erections to achieve the satisfactory sexual activity (Muneer et al., 2014). ED is a common disease with high incidence worldwide, and there are few therapeutic methods. It is associated with aging, diabetes, and postoperative complications. The etiology of ED is mainly vascular, caused by nitric oxide (NO) metabolism disorder in the cells. NO enhancing drugs like phosphodiesterase 5 inhibitors (PDE-5i), which alleviate NO's decay, are widely utilized but still fail sometimes (Altabas & Altabas, 2015). The dysfunction of endothelial cells (ECs) plays a vital role in the pathophysiological process of ED. The complex pathogenesis of diabetic ED, including nerves, blood vessels, endothelium, smooth muscle, other tissues, and organs, have been long studied. As an independent factor, hyperglycemia leads to vascular endothelial cell injury and dysfunction. Vascular endothelial integrity and barrier function are destroyed through oxidative stress injury, leading to vascular disease.
Autophagy is an evolutionarily conserved eukaryotic self-catabolic mechanism. It balances energy sources, clears misaggregated or folded proteins, and clears damaged organelles (such as mitochondria, endoplasmic reticulum, and peroxisome). It is involved in many biological functions including development, in response to nutrient stress, cellular differentiation, and resistance to pathogens. Therefore, autophagy serves as a survival mechanism, although its deregulation is a related to non-apoptotic cell death. There is growing evidence that mutations or suppression of key autophagy genes are associated with cancer, neuropathy, cardiovascular disease, autoimmune diseases, and other diseases (Glick, Barth & Macleod, 2010;Levine & Kroemer, 2019). Studies have confirmed that endothelial autophagy plays a critical role in the occurrence and development of cavernous endothelial cell lesions in ED (Lin et al., 2018;Zhang et al., 2019a). However, the role and specific regulatory mechanism of the autophagy gene in ECs remains unclear.
In this study we expected to analyze the human ED dataset in the GEO database using bioinformatics to explore the autophagy-related genes and related molecular mechanisms connected with endothelial dysfunction. Then, NOTCH1, PPARG, NOS3, KEAP1, CDKN2A, and IRS1 were predicted to affect endothelial dysfunction of ED as autophagy-related hub genes. Additionally, NOTCH1, CDKN2A, and NO3 were validated by establishing rat model of diabetic ED (DMED), and may be the underlying target molecules to treat ED by modulating autophagy.

Identification of differentially expressed genes (DEGs)
The dataset was classified into two groups: the ED group (four samples of human corpus cavernosum ECs from donors with ED) and the non-ED group (including one sample of human corpus cavernosum ECs from donors without ED, and seven samples of human arteriovenous ECs from donors without ED). The limma package in R was used with the Wilcoxon test to screen out the significant DEGs. According to the annotation information of GPL571 platform, the probes were transformed to the corresponding gene symbol. Probe sets without corresponding gene symbols were excluded and genes with more than one probe set were averaged. The cut-off values were set according to the parameters, an absolute logFC >0.5, and false discovery rate (FDR) < 0.05.

Enrichment analysis and gene-concept network
We explored the correlation between all ARGs, up-regulated ARGs, down-regulated ARGs, and AREs using R. We focused on the functional enrichment analysis of Gene Ontology (GO), including the biological process, cellular component, and molecular function. We also analyzed the Kyoto Encyclopedia of Genes and Genomes (KEGG) of 20 ARGs.
We introduced the ''cnetplot'' to depict the linkages of genes and biological concepts of GO to form the Gene-Concept network (Yu et al., 2012).

Interaction network of hub genes
The ARGs interaction network was created in STRING (http://string-db.org/cgi/input.pl) using the PPI network with the cutoff criteria as a combined score >0.4. Cytoscape (3.7.2) software was used to analyze and visualize the PPI network for hub genes with a cutoff value of >4.

Animal and experiment designs
Eight-week old male Sprague-Dawley (SD) rats, weighing approximatley 200 g, with normal erectile function were provided by the Experimental Animal Center of Southern Medical University. The Institutional Research Ethics Committee of Nanfang Hospital of Southern Medical University provided full approval for this research (NFYY-2020-64). The 3R principle was used in our animal experiments: reduction of animal numbers, replacement of animals using other methods, and refinement of animal welfare. The experimental rats were randomly divided into the DMED group and the control group (NC). The rats with DMED were modeled and identified as previously mentioned (Zhang et al., 2019b). Briefly, healthy adult male SD rats were intraperitoneally injected with a 1% streptozotocin solution (65 mg kg −1 ). Diabetes was determined by measuring random tail vein blood glucose levels 72 h after injection. Rats with random blood glucose concentrations >16.7 mmol L −1 were diagnosed as diabetic. Six DMED rats were included in this study; six control rats were not treated. We used a blood glucose meter (Anwen, Changsha, China) to measure random blood glucose in the tail vein blood every two weeks. Weight was also evaluated every two weeks.
All rats were reared in the Specific Pathogen Free (SPF) temperature-controlled animal house of Nanfang hospital's animal center with a 12-hour light-dark cycle. A specialized breeder changed the padding once a day and fed the rats with warm water that had been boiled and ordinary food (produced by Yancheng Biotechnology Co., Ltd., Guangzhou, China).

Erectile function evaluation
The mean ICP and ICP/MAP ratio were applied to evaluate erectile function, as mentioned above (Zhang et al., 2017). All rats were anesthetized with an intraperitoneal injection of sodium pentobarbital (30 mg/ kg). The cavernous nerves were completely exposed and the corpus cavernosum was carefully isolated. A 25-G needle containing a 100 U/ml heparin solution was slowly inserted into proximal corpus cavernosum to ensure that the rat was securely connected with the sensor and amplifier (MP150 Biopac System; Biopac Systems Inc., CA, USA). Both were connected using AcqKnowledge R (V4.4) software. We recorded the intracavernous pressure (ICP) when the cavernous nerve was stimulated with bipolar stainless steel electrodes. The mean atrial blood pressure (MAP) was measured by intubation in the biological signal system above after exposing the left carotid artery. A successful DMED rat model should meet an ICP of less than 60 mmHg and a ICP/MAP ratio less than 0.5. Data were recorded and the penis shaft was collected for future study. All rats were euthanized using a 10-fold anesthetic dose (300 mg/kg sodium pentobarbital).

Histology
Freshly dissected tissue was fixated. H&E and Masson's trichrome staining was performed to determine tissue structure changes based on the manufacturer's instructions. We also determined the ratio between smooth muscle and collagen in the corpus cavernosum. Images were captured using Olympus BX63 microscopy (Olympus, Tokyo, Japan).

Immunofluorescence analysis and laser confocal microscope
Tissue slides were prepared in a smiliar manner to those used in H&E staining. After antigen repair, goat serum was inoculated at room temperature for 30 min, the CD31 antibody (ab182981; Abcam) was incubated overnight, and a secondary antibody solution of Goat Anti-Rabbit IgG (HRP) (ab205718; Abcam) was added. Antigen repair was performed after the first dyeing and cleaning. Donkey serum was inoculated at room temperature for 30 min, the LC3 antibody (1:200; Proteintech, USA) was incubated overnight, and Donkey Anti-Rabbit 594 secondary antibody (Abcam, ab150075) was added through a drip. We used 4 , 6-Diamidino-2-phenylindole (DAPI, #C0060; Solarbio) to stain the cell nuclei and images were captured using laser confocal microscopy (Nikon, Tokyo, Japan).

Statistical analysis
The gray value and fluorescence intensity of image were further analyzed using Image J software. The outcomes were shown as the mean ± standard error. One-way ANOVA was applied to perform statistical comparisons among the groups. R software (version 4.0.0) and GraphPad Prism 8.3.0 were used to perform all statistical analyses. A P-value of <0.05 was considered statistically significant.

Identification of differentially expressed autophagy-related genes (ARGs)
We retrieved a total of four ED samples and eight non-ED samples with gene expression profiles from the GEO dataset. The RNA degradation plot was introduced to measure the quality of the RNA in the raw-data sample. The fluorescence intensity at the 5 end of the chip was much lower than that at the 3 end, and the slopes of the curves ranked in a similar order (Fig. 1A). We normalized the raw data for further analysis (Fig. 1B). We identified 1,236 significant DEGs between the ED and non-ED samples. Among these DEGs, 423 genes were up-regulated in ED tissue compared with non-ED tissue; the other 815 genes were down-regulated via Volcano Plot and heatmap (Fig. 1C). We matched DEGs with 380 AREs from the GeneCard (File S1) using a Venn diagram (Fig. 1D), resulting in 20 ARGs. Thirteen ARGs were down-regulated in ED tissues; the other seven genes were up-regulated compared with non-ED tissues (Table 2, Fig. 1E).

Functional enrichment analysis and interactions of ARGs
All twenty ARGs were used for functional enrichment analysis to explore the biological functions and pathways of ARGs in autophagy. GO results showed that the biological process (BP) of ARGs were mainly enriched in autophagy, a process utilizing autophagic mechanism, the regulation of autophagy, positive regulation of proteolysis, and the positive regulation of the cellular catabolic process. The cellular component of ARGs was enriched in the nuclear outer membrane. ARGs' molecular functions were enriched with oxidoreductase activity, acting on NADH, heme protein as acceptor, enzyme inhibitor activity, Hsp70 protein binding, ligand-activated transcription factor activity, and nuclear receptor activity (Figs. 2A, 2C). The results of KEGG enrichment of ARGs were significantly enriched in the mTOR signaling pathway, PI3K-Akt signaling pathway, longevity regulating pathway, AMPK signaling pathway, and platelet activation, excluding autophagy-related pathways (Fig. 2B).  The enriched GO pathways of molecular functions of seven up-regulated ARGs were primarily nuclear receptor activity, ligand-activated transcription factor activity, and steroid hormone receptor activity (Fig. 2D). GO enrichment of 13 down-regulated ARGs showed that molecular functions were enriched in oxidoreductase activity, acting on NAD(P)H, heme protein as acceptor, enzyme inhibitor activity, and Hsp70 protein binding (Fig. 2E). The GO enrichment of 380 AREs showed that molecular functions were enriched in ubiquitin-protein ligase binding, ubiquitin-like protein ligase binding, protein serine/threonine kinase activity, phosphatase binding, protein phosphatase binding, and heat shock protein binding (Fig. 2F). Interestingly, heat shock protein (HSP) binding appeared in three groups.

Blood glucose, weight, ICP/MAP, and morphological changes of penis corpus cavernosum of DMED model
The blood glucose of DMED rats was significantly higher than that of the NC rats, while the weight of DMED rats was significantly lower than that of the NC rats at the 8th week after modeling (Fig. 4A). We successfully established the DMED rat model with ICP less than 60 mmHg and ICP/MAP ratio less than 0.5 (Fig. 4B). Results from H&E and Masson's trichrome staining revealed that the corpus cavernosum was disordered and the ratio of fibroblast to collagen decreased when compared with the NC group. The endothelium and smooth muscle layer in the DMED group were thinner than the NC group (Figs. 4C-4D).   Node color reflects the degree of connectivity (dark-green color represents a higher degree, and lightgreen color represents a lower degree)-networks of (C) target gene-miRNA and (D) target gene-TF. The red circle nodes are the genes, green hexagon nodes are the miRNAs, and green diamond nodes are the TFs.

Enhanced autophagy of cavernous tissue and Endothelium in DMED model
The protein level of LC3-II and LC3-I, and the ratio of LC3-II/LC3-I were upregulated, and the grayscale value of Beclin-1 was downregulated in DMED group compared with the NC group (Fig. 5A). The fluorescence intensity of LC3 in the DMED group was significantly higher than that in the NC group. The increase of LC3 represented an increase in autophagosomes. The discontinuity of CD31 labeled ECs, the decrease of fluorescence intensity, and the disappearance of fluorescence on smooth muscle surface suggested that the ECs in DMED tissue were injured (Fig. 5B).

The expression level of hub ARGs
The results of qRT-PCR showed that the expressions of Notch1, PPARG, NOS3, CDKN2A, and IRS1 were down-regulated, however, there was no statistical significance in the expression of KEAP1 compared to NC group (Figs. 6A-6F). The gene expression of NOTCH1, CDKN2A, and NOS3 was consistent with the data in Table 2. The protein level of NOTCH1, CDKN2A, and NOS3 were downregulated in DMED group compared to NC (Fig. 6G).

DISCUSSION
Endothelial dysfunction is one of the leading causes of ED. ECs play a critical role in regulating inflammation, platelet aggregation, vascular smooth muscle hyperplasia, and thrombosis. Thus, targeting endothelial dysfunction may be a more effective treatment for ED. Studies have suggested that endothelial dysfunction may be alleviated by inhibiting autophagy (Altabas & Altabas, 2015;Verma & Anderson, 2002). Autophagy is a multi-step collaborative process regulated by autophagy-related genes (Sarkar et al., 2014). However, the pathway of autophagy-related genes and ED remains unclear. In this study, we utilized bioinformatics to mine twenty autophagy-related genes and their related molecular functions in an ED database. Moreover, six hub genes were further screened out using the PPI co-expression network. We further verified these genes by constructing rat model of diabetic ED relying on the homology of human and mammalian genes. The LC3-II/LC3-I ratio and LC3-II/ β-actin ratio in the DMED group was relatively higher than the control group. LC3-II is bound to the cell membrane in mature autophagosomes, is released upon fusion with the lysosome, and is considered to be a marker for autophagy detection (Ni et al., 2011). This suggests that autophagy is enhanced in DMED tissue. In addition, the enhanced autophagy of endothelial cells was confirmed by CD31 and LC3 co-staining. Excessive ROS accumulation promotes compensatory autophagy during the pathogenesis of diabetic ED, which protects the body from oxidative stress by clearing damaged intracellular substances (Scherz-Shouval et al., 2007). The Beclin1/ β-actin ratio in the DMED group was lower than in the control. This protein plays a central role in the formation and maturation of autophagosomes. The conditions that disrupt or promote the Beclin1-Bcl2 complex play a crucial part in determining whether cells undergo autophagy or apoptosis. A recent study showed that when Beclin1 was induced it facilitated autophagy in various cell types, including ECs (Leonard et al., 2019). We found that the gene expression of NOTCH1, CDKN2A, and NOS3 was consistent with the high throughput sequencing results from GEO. These may become underlying target molecules in the regulation of EC autophagy in ED.
Notch is a classical signaling pathway, including four Notch receptors (Notch 1, 2, 3 and 4) as well as five notch ligands (delta-like receptors 1, 3, 4 and jagged 1 and 2) in mammals (Gordon, Arnett & Blacklow, 2008). The pathway is activated when a Notch receptor interacts with a ligand. Studies suggested that the inhibition of autophagy leads to the activation of Notch pathway. Autophagy suppresses the Notch1 pathway by up-regulating the degradation of Notch1 (Li et al., 2016;Wu et al., 2016;Zeng et al., 2016). The reduction of NOTCH1 in ECs is a predisposing factor for vascular inflammation and atherosclerosis (Briot et al., 2015). Endothelial dysfunction may be associated with a decrease of Notch1. Moreover, the Notch pathway was associated with senescence which the mainstream factor of ED. A recent study showed that the Notch-1 inhibitor, Tangeretin, alleviated ED by increasing the maximum ICP/MAP and up-regulating the expression of eNOS in hypertensive rats (Chiangsaen et al., 2020). Thus, NOTCH1 may be a candidate biomarker for ED.
CDKN2A, also called p16 INK4a (p16), has been shown to increase with age in several rodent and human tissues (Baker et al., 2011). As a member of the INK4 family of cyclin-dependent kinase inhibitors, p16 plays a critical role in cell-cycle regulation. p16 expression inhibits cellular proliferation by downregulating cyclin-dependent kinases 4 and 6 (CDK4/6) (Serrano, 1997). As a tumor suppressor, p16 INK4a (p16) is a well-studied maturation marker, inducing cellular senescence (permanent cell-cycle arrest) in response to stress. p16 accumulates in tissues with age, causing cellular senescence. The elimination of p16 expression in senescent cells is associated with a prolonged life span and a reduction of tumorigenesis (Baker et al., 2011). Studies have suggested that autophagy inhibitors bafilomycin A1 and chloroquine induced p16 accumulation in stagnant vesicles containing the autophagy marker LC3. The accumulation of p16 in these vesicles was consistent with the increase of LC3-II. The knockout of autophagosome chaperone p62 attenuated p16 aggregates in lysosomes, indicating that p16 targets these vesicles through p62. As the regulator of autophagy, p16 performs a key role in the etiology of cancer and dementia (Coryell et al., 2020). Moreover, advanced age is a critical risk factor for most chronic diseases and functional defects in humans. A significant increase has been shown in ED incidence in men over the age of 40 (Shamloul & Ghanem, 2013). Senescence is closely associated with ED. We found that the expression of CDKN2A was decreased. This result suggests that not all cells with high expression of p16 Ink4a are senescent cells and not every senescent cell has high expression of p16 Ink4 (Hall et al., 2017). There should be further studies conducted on CDKN2A as an autophagy-related molecule.
Endothelial nitric oxide synthase (eNOS) is also known as NOS3, one of three subtypes of nitric oxide synthase (NOS). It is involved in the synthesis of nitric oxide (NO) with L-arginine and molecular oxygen as substrates. It is also involved in the regulation of physiological and pathological functions. The pathophysiology of endothelial dysfunction involves multiple mechanisms, such as the dysregulation of NO by vascular/endothelial eNOS in ED.
eNOS performs an essential role in erectile response. eNOS activity and endothelial NO bioavailability in the penis are jointly regulated by multiple post-translational molecular mechanisms, including eNOS phosphorylation, eNOS interaction with upstream proteins and contractile pathways, and the generation of reactive oxygen species (ROS). The activation of the PI3K/Akt/eNOS pathway was considered to be one of classic pathways to regulate endothelial dysfunction and participate in the cellular and molecular biology of ECs (Chen et al., 2021;Chen et al., 2019;Duan et al., 2019;Zhang & Zhang, 2020). The regulation of the AMPK pathway and the mTOR pathway corresponded to endothelial autophagy (Xiong et al., 2014). Our results confirmed that autophagosome increased in the cavernous ECs of DMED rats, suggesting an increase in autophagy. The participation of eNOS in the autophagy pathway in endothelial dysfunction requires further experimental study.
Previous literature reporting on mammalian cells has indicated that protein degradation pathways are commonly classified in the ubiquitin-proteasome system. There are three main types of autophagy: macro-autophagy, micro-autophagy, and chaperone-mediated autophagy. Hsp70 acts a central hub in protein degradation through the ubiquitinproteasome system and different autophagy pathways (Fernández-Fernández et al., 2017). Our GO analysis suggested downregulated ARGs and 380 AREs were both associated with HSP protein binding. The Hsp70 family performs central roles in every aspect of proteostasis, from protein folding to disaggregation and degradation (Mayer & Bukau, 2005). Targeting HSP70 binding in ED may be a new direction for future research. Furthermore, the upregulated ARGs (AR and PPARG) were mainly enriched in nuclear receptor activity and downregulated ARGs (TSC1 and SNCA) were mainly enriched in oxidoreductase activity. However, their relationship with ED requires further study.
The top three targeted ARGs in the miRNA-gene network were AR, NGFR, and FLCN. It has been suggested that nerve growth factor (NGF) may induce nerve regeneration by activating the autophagy of Schwann cells. It may also promote the maintenance, survival, and function of vascular endothelial cells (Tanaka et al., 2004). Low doses of rapamycin, an autophagic inducer, accelerated autophagy by activating the NGFR promoter. Therefore, NGFR may be involved in the regulation of endothelial autophagy. Hsa-mir-24-3p and Hsa-mir-335-5p, which control more ARGs, perform a key role in angiogenesis (Esquinas et al., 2017;Li et al., 2020). The relationship between AR, FLCN, and endothelial dysfunction requires further study. NOS3, VAMP8, and MLST8 were ranked as the top three targeted ARGs in the TF-gene network. A previous study showed that platelets stimulated BMDC homing and promoted angiogenesis. Thus they are a key mediator between hypoxic tissue and bone marrow, in which VAMP8 plays an important role (Feng et al., 2011). In endothelial cells, MLST8 is involved in modulating the activation of Akt and reducing the phosphorylation of Akt targets such as the eNOS and FOXO subfamily (Partovian et al., 2008). The role of NOS3 in endothelial dysfunction has been discussed.
Recently there has been great interest in the relationship between autophagy and ED. Glucagon-like peptide-1 (GLP-1) receptor agonists, defocused low-energy shock wave therapy, and stem cell therapy have been suggested to alleviate ED by promoting endothelium autophay in DMED rats (Yuan et al., 2020;Zhang et al., 2019a;Zhu et al., 2018). Autophagy was shown in the ED model of BCNI and aging. Thus ED may be allieviated by improving autophagy (Tang et al., 2018;Ye et al., 2021). However, there is no research on the expression changes and molecular functions of ED-related autophagy genes. Our study may provide new directions for future research and novel molecular targets for the treatment of ED.
Our research was limited by the small number of microarray samples which may lead to statistical bias. Secondly, we screened genes based on microarray data sourced from human tissue in GEO. However, this was verified on SD rats which may cause species bias. We have also not specifically studied whether these genes can regulate autophagy.

CONCLUSIONS
Our study mined endothelial-related ARGs to identify the possible pathogenesis of ED related to autophagy. Consistent with the fact that high glucose accelerates the senescence of endothelial cells by inhibiting autophagy, two of three identified genes, NOTCH1 and CDKN2A, are involved in senescence. NOTCH1, CDKN2A, and NO3 may be the underlying therapeutic targets to treat ED. Moreover, the study of upstream transcription factors and miRNAs of target genes may improve the mechanistic understanding of endothelial dysfunction in ED.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the National Natural Science Foundation of China (82060809). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.