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

ED datasets

Microarray data were obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) (Edgar, Domrachev & Lash, 2002). Specifically, we downloaded the gene expression profile of GSE10804 (Affymetrix GPL571platform, Affymetrix Human Genome U133A 2.0 Array) (Table 1). This dataset contained four samples of human corpus cavernosum ECs from donors with ED, one sample of human corpus cavernosum ECs from donors without ED, and seven samples of human arteriovenous ECs from donors without ED (three human umbilical vein ECs and four human coronary artery ECs).

RNA quality was measured using the “RNA degradation map” provided by the Affy package in R language (4.0.0) (Gautier et al., 2004; Wilson & Miller, 2005). The dataset was standardized usung the limma package in R.

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.

Identification of autophagy-related DEGs (ARGs)

We obtained the autophagy-related genes (AREs) from GeneCard (https://www.genecards.org/). We entered the keywords “autophagy”, used the relevance score >5.0 as the cutoff value, and downloaded them in the txt file format. Then, the DEGs and autophagy-related genes were entered into a Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) to obtain the ARGs (Wang, Thilmony & Gu, 2014).

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.

Construction of target gene-miRNA network and target gene- transcription factor (TF) network

The expression of target genes in the post-transcriptional stage was regulated by miRNA or TF under specific disease conditions. Baldwin Jr, 2001; Liu et al., 2020; Soifer, Rossi & Saetrom, 2007) NetworkAnalyst (https://www.networkanalyst.ca/) was used to analyze gene-TF interaction (ENCODE database) and gene-miRNA interaction (miRTarBase v8.0 database) (Xia, Gill & Hancock, 2015) and was visualized using Cytoscape.

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⁻¹). 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⁻¹ 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^® (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).

Real-time quantitative PCR (qRT-PCR)

Cavernous tissue RNA was extracted and used for qRT-PCR analysis. The primer sequences we used were: β-Actin, sense: 5′-GATCA-AGATCATTGCTCCTCCTG-3′, anti-sense: 5′-AGGGTGTAAAACGCAG-CTCA-3′; NOTCH1, sense: 5′-CAGTACAACCCGCTAAGGC, anti-sense: 5′-GGACAAGG-TATTGGTGGAGA-3′; NOS3, sense: 5′-CCGATTACACGAC- ATTGAGA-3 ′, anti-sense: 5′-TGGTCCAGTTGGG-AGCAT-3′ (anti-sense); KEAP1, sense: 5′-GGTCGC-CCTGTGCCTCTAT-3′, anti-sense: 5′-CACGCTGCTGTGGTGGAT-3′; CDKN2A, sense: 5′-GAGGGCTTCCTAGACACTCTG- 3′, anti-sense: 5′-CGCAAATACCGCA-CGAC- 3′; PPARG, sense: 5′-CCTCCCTGATGAATAAAGA-3′, anti-sense: 5′-AAC-TCAAACTTAGGCTCCA-3′; IRS1, sense: 5′-GAGTGGTGGAGTT- GAGTTGG-3′, anti-sense: 5′-GTCCGCATGTCAGCATA-3′.

Western blot (WB)

WB was carried out for protein expression analysis as previously described (He et al., 2012). We detected LC3 and Beclin-1 using a rabbit anti-LC3 antibody (14600-1-AP) (1:1,000), anti-Beclin-1 antibody (11306-1-AP; Proteintech) (1:1,000), anti-NOTCH1 antibody (10062-2-AP; Proteintech) (1:1000), and anti-eNOS (Abcam, ab199956) (1:1000). Secondary antibodies were purchased from LI-COR Biosciences (catalog #D00825-14) and were diluted at 1:15,000. β-actin (AC026; ABclonal) was used as a loading control.

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.

Construction of the PPI network, gene-microRNA network, and gene-TF network

We explored the twenty ARGs interactions using the STRING website (Fig. 3A). The gene network demonstrated that IRS1, CDKN2A, NOTCH1, PPARG, KEAP1, and NOS3 were the hub genes using Cytospace software (Fig. 3B). The top three ARGs for miRNAs were AR, NGFR, and FLCN. Hsa-mir-24-3p and Hsa-mir-335-5p may control the largest number of ARGs (Fig. 3C). NOS3, VAMP8, and MLST8 were the top three target ARGs for TFs, modulated by 78, 77, and 70 TFs (Fig. 3D).

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).

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).