We used 22 geriatric male SD rats (20 months) purchased from the Experimental Animal Center of Nanjing Medical University, and our research was approved by the Animal Care and Use Committee of the First Affiliated Hospital of Nanjing Medical University. The animals were anesthetized with pentobarbital natrium (40 mg/kg) by intraperitoneal injection.

After anesthesia, dissections were made along the penis midline to expose the carotid artery. A neck artery was cannulated and a heparinized 25-gauge needle was inserted into the corpora cavernosum of the penis to measure both mean arterial pressure (MAP) and ICP. Electrical stimulation was set at 5 V, 15 Hz, pulse width 0.2ms, and the stimulation lasted for 1 min. The BL-420S Biological Functional System (Chengdu Taimeng Technology Co, LTD, Chengdu, China) was used to measure pressure and generate electrical stimulation. Each experiment was done in triplicate. The ratio of the maximal ICP to the corresponding MAP (ICP/MAP) was calculated to estimate the erectile function. The corpus cavernosum tissue was isolated and stored -80°C for further experiments.

Experiments were carried out as previously described (21). Penile tissue was fixed with 4% paraformaldehyde and embedded in paraffin. Masson trichrome staining was used to observe, under 200x magnification, then evaluate the ratio of smooth muscle to collagen area in penile tissue.

The rats were divided into geriatric A-ED (erectile dysfunction) and NC (negative control) groups based on the erectile evaluation results. ICP/MAP > 35% was defined as the geriatric NC group, and ICP/MAP < 35% was defined as the geriatric ED group. After stripping the urethra and other connective tissue from the penis, the cavernous tissue was retained. Cavernous tissue was minced into small pieces with a low-temperature grinder (Servicebio, KZ-III-F) under grinder parameters of 60Hz & -10°C. Then, total RNA was extracted from pieces of cavernous tissues or cells using a TRIzol (Invitrogen, Carlsbad, CA, USA) based protocol. RNA concentration and purity were determined via a NanoDrop 2000 spectrophotometer (Madison, WI, USA). Ribosomal RNA (rRNA) was separated from total RNA and RNA integrity was analyzed via an Agilent RNA 6000 Nano Kiton an Agilent 2100 Bioanalyzer. Residual total RNA was broken into short fragments for 2 min, then reverse transcribed into cDNA with random hexamers as primers via a SuperScript III First-Strand Synthesis protocol (Invitrogen). Synthetized products were purified using a QIAquick PCR Purification Kit (QIAGEN). DNA fragment ends were repaired and ligated with poly (A) tails to connect them to Illumina sequencing adapters. Then, the second strand was digested and the appropriate fragment was identified for PCR amplification. The cDNA sequencing library was made and subsequently sequenced using an Illumina HiSeq TM 4000 platform (Gene Denovo Biotechnology Co., Guangzhou, China). Only reads that met the following criteria were used as raw data (1) reads did not contain adapters; (2) reads did not contain more than 10% of unknown nucleotides (N); (3) reads did not contain more than 50% of low quality (Q-value ≤ 20) bases. These reads were aligned to Ribosome RNA (rRNA) and mapped to a rat reference genome using HISAT2 (22).

We extracted total RNA from previously isolated corpus cavernosum tissues using a Trizol (Invitrogen, USA) based protocol. Then, cDNA was synthesized using HiScript^® III All-in-one RT SuperMix Perfect for qPCR (Vazyme, China). Next, we utilized the StepOne Plus Real-Time PCR system (Applied Biosystems, USA) to perform qPCR. Fold changes in mRNA expression were calculated using the 2−ΔΔCt method and normalized against β-actin with ABI Step One software version 2.1. PCR primer sequences were synthesized by TSINGKE Biological Technology (Nanjing, China) and are listed in Supplementary Table S1 .

Transcriptome data were normalized and analyzed in R using the edgR package. Then three validated lncRNAs were selected (ENSRNOT00000090050, ENSRNOT00000076482, and ENSRNOT00000029245). For prediction analysis, we defined differently expressed miRNAs and mRNAs as logFC >2 and p <0.05. MiRNAs target prediction tools TargetScan, miRcode, and miRDB were used to predict miRNAs targets. All relevant annotations of miRNAs were derived from miRbase.

The network was constructed based on ceRNA theory as previously described (21). (1) Spearman rank correlation coefficient (SCC) was used to evaluate expression correlation between miRNAs-mRNAs or miRNA-lncRNAs. RNA interaction pairs that had SCC<0 were defined as either negatively co-expressed lncRNA-miRNA or mRNA-miRNA pairs. All RNAs were differentially expressed. The mRNA and lncRNA were predicted as miRNA target genes. (2) Pearson correlation coefficient (PCC) was used to assess expression correlations between lncRNA and mRNA. ceRNA interaction pairs that had PCC > 0.5 were defined as co-expressed lncRNA-mRNA pairs. RNA interaction pairs of mRNA and lncRNA were termed as negatively or postively co-expressed with a common miRNA. Pairs of mRNA and lncRNA were targeted negatively or positively co-expressed with a common miRNA. (3) The significance of common miRNA sponges differences between the two genes was tested by a hypergeometric cumulative distribution function test. The filtering condition was p< 0.05. Finally, all eligible RNA interaction pairs were collected and combined into a network, which was visualized in Cytoscape 3.8.2 (http://www.cytoscape.org/).

Gene ontology (GO, http://www.geneontology.org/) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) pathway analysis of differentially expressed mRNAs were performed using the R package clusterProfiler. An alpha level of p<0.05 was used to test statistical significance.

The rat vascular endothelium cell line (RAOEC) was purchased from Shanghai Jingkang Biological Co., LTD (Shanghai, China) and cultured at 37 °CC in 5% CO₂. The endothelial cell lines were characterized by double-immunofluorescent staining for biomarkers CD31 and vWF (Proteintech, China). DMEM medium, supplemented with 10% FBS, and 100 units/mL penicillin, were used in cell culture (Gibco).

ROAEC cells were treated with hydrogen peroxide (H₂O₂) to induce senescence (23, 24). Cells were seeded onto 6-well plates (10⁶/well) in a full growth medium. We choose the optimized concentration (100 µM) to build a cellular aging model according to the references (25). After 24 hours, the cells were treated with H₂O₂ (100μM) and DMSO (0.01%) respectively for 24 hours. Senescence-associated β-galactosidase (sa-β-Gal) staining was performed using a β-galactosidase enzyme assay kit (Beyotime, China) between the H₂O₂-treated group and the DMSO group.

All data analyses were done using SPSS software (SPSS V 24.0). T-tests were used to determine the statistical significance of differences between two groups of parametric data whereas Mann-Whitney tests were used for nonparametric data. Statistical significances of GO and KEGG analyses were determined using Fisher’s exact tests. An alpha level of p<0.05 was used for all statistical tests.

Among 22 rats, we selected ICP/MAP less than 0.35 as the low-pressure group and ICP/MAP greater than 0.35 as the high-pressure group, and when combined with the results of the Masson staining, 22 rats were divided into A-ED group and NC group (n=7, n=15). There was no significant difference in MAP between the A-ED and NC groups (p=0.3777). However, the ratio of ICP/MAP was lower in the A-ED than the NC group ( Figures 1A, B ). Masson trichromatic staining results showed that the smooth muscle/collagen ratio was significantly lower in the A-ED than the NC group ( Figures 1C, D ). Thus, this A-ED rat model was valid for studies on age-related erectile dysfunction.

A total of 5,307 lncRNAs, 1,322 miRNAs, and 22,050 mRNAs were identified by sequencing. The expression levels of lncRNAs in six samples were assessed by FPKM, and did not show a deviated distribution of mRNAs nor lncRNAs ( Figures 2A, B ). Scatterplots ( Figures 2C, D ) demonstrated high reproducibility of samples from each group.

Differentially expressed mRNAs and lncRNAs were analyzed using volcano plots ( Figures 3A, B ). Upregulated and downregulated mRNAs and lncRNAs, when compared to the geriatric NC group, were shown using orange and blue dots respectively. Ultimately, a total of 69, 364 and 29 differential expressed lncRNAs (49 up-regulated, 20 down-regulated), mRNAs (143 up-regulated, 221 down-regulated), and miRNAs (18 down-regulated, 11 up-regulated) respectively, were identified ( Figures 3C–E ).

Using an inclusion criteria for lncRNAs of 1<count number<120, fold change > 2, p< 0.05, seven lncRNAs were identified and used to validate RNA-seq results in seven pairs of geriatric ED and NC samples (MSTRG.3646.1, ENSRNOT00000085383, ENSRNOT00000093493, ENSRNOT00000081965, ENSRNOT00000029245, ENSRNOT00000090050, ENSRNOT00000076482, ENSRNOT00000085383). According to qPCR results, ENSRNOT00000076482(p=0.0175), ENSRNOT00000090050(p=0.007) and ENSRNOT00000029245(p=0.0006) were highly expressed in the ED group, corroborating sequencing results ( Figure 4 ). However, expression levels of the other four lncRNAs did not statistical differ between groups.

To further explore the specific function of these lncRNAs, related mRNAs were used to carry out GO and KEGG enrichment analyses and the top 20 results of enrichment analysis (including KEGG and GO) according to p values, are shown in Figure 5 . In doing KEGG analyses, we identified pathways for PI3K-AKT & Jak-STAT signaling, terpenoid backbone biosynthesis, synthesis and degradation of ketone bodies, EGFR tyrosine kinase inhibitor resistance, fluid shear stress and atherosclerosis.

From GO enrichment analyses, cellular component (CC) annotations showed that the following components played key roles: coated pit (GO:0005905), striated muscle myosin thick filament (GO:0005863), dehydrodolichyl diphosphate synthase complex (GO:1904423), muscle myosin complex (GO:0005859), and phosphatidylinositol 3-kinase complex, class IB (GO:0005944). Corresponding annotations of molecular functions included translation initiation factor binding (GO:0031369), hydroxymethylglutaryl-CoA synthase activity (GO:0004421), and interleukin-6 receptor activity (GO:0004915).

Moreover, from the annotations of biological processes, the following were enriched ( Figure 5 ): response to purine-containing compound (GO:0014074), regulation of long term synaptic depression (GO:1900452), muscle cell cellular homeostasis (GO:0046716), regulation of platelet-derived growth factor receptor signaling pathway (GO:0010640), and organophosphate biosynthetic process (GO:0090407).

From the results of RT-qPCR and RNA-seq, three verified lncRNAs, 29 differentially expressed miRNAs, and 364 differentially expressed mRNAs were selected, and the regulatory network was constructed according to its hypothesis. Figure 6 illustrates the lncRNA-miRNA-mRNA regulatory network of geriatric ED rats. In our network, the up-regulated lncRNA ENSRNOT00000029245 upregulated the expression of 23 mRNAs (Tspan18, Mtss1, Szrd1, Pik3cg, NEWGENE_619861, Fgf1, Btrc, B2m, Ggta1, kdm4a, Fxyd7, Naf1, Vps9d1, Ppm1l, Bcl2l1, RGD1560281, Scamp1, Sele, Nrcam, Hsp90b1, Cmpk1, Fam20b, Hiatl3) by miR484x axis. Similarly, the lncRNA ENSRNOT00000029245 regulated the 14 mRNAs (Sorcs3, LOC100910882, Morf4l2, Ptprd, Hmgcs1, Lox, Gpm6a, Trim32, Tnfrsf11b, Phf14, Ncapg2, Cyyr1, Reep3, Nus1) by rno-miR-653-5p axis. Interestingly, according to our network, the mRNA Il6r was upregulated by both the miR484x axis and the rno-miR-653-5p axis.

The RAOEC cells were identified by staining with anti-CD31 antibody (red), anti-vWF antibody (green), and cell nuclei were identified using DAPI (blue) ( Figure 7A ). The sa-β-Gal staining demonstrated that the RAOEC cells were induced senescence at 100μM H₂O₂ ( Figure 7B ). Interestingly, the expression of the three lncRNAs (ENSRNOT00000029245, ENSRNOT00000076482, and ENSRNOT00000090050) between cell lines showed that only ENSRNOT00000029245 had a statistical significance ( Figure 7C ).