Identification and Validation of an Aging-Associated circRNA-miRNA-mRNA Network in Neovascular Age-Related Macular Degeneration Open Access

Neovascular age-related macular degeneration (NVAMD) is one of the leading causes of irreversible vision impairment in the elderly population [1]. It is characterized by choroidal neovascularization (CNV) and secondary intraretinal or subretinal leakage and hemorrhage [2]. The occurrence of NVAMD increases substantially with age. As predicted, the prevalence of AMD increases from 0.04 to 5.48% in patient groups aged 50 and 85, respectively [3]. Anti-vascular endothelial growth factor (VEGF) treatment is currently the first-line therapy; however, not all patients respond well [4]. Therefore, a comprehensive investigation of AMD pathologies is warranted to develop novel or complementary therapies.

As a multifactorial disease, many pathological mechanisms, such as inflammation and oxidative stress, have been demonstrated to cause retinal pigment epithelium (RPE) cell dysfunction and death in previous studies [5, 6]. Aging is a predisposing factor correlating with these pathophysiologies. Senescent cells are detected in the retinas and RPEs of older people and primates [7, 8]. In addition, the cytokine profile is specific for aging in human primary RPE cell cultures and RPE cells from donor AMD eyes [9]. Moreover, eliminating senescent RPE cells alleviates the CNV volume in laser-induced CNV mouse models [10]. Therefore, RPE cell senescence is a promising target in novel NVAMD strategies.

Noncoding RNAs, including microRNAs (miRNAs), long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs), are vital regulators of aging [11]. For example, ectopic expression of miR-182 triggers epithelial cell senescence through the p53/p21 pathway [12]. CircRNAs are noncoding RNAs with a specific loop structure that act by sponging miRNAs. For example, circRNA CCNB1 sponges miRNA449a and further prevents CCNE2 expression, inhibiting the senescence process [13]. These studies suggest that aging-based noncoding RNA networks play a potential role in aging.

This study analyzed the NVAMD microarray dataset from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) in this analysis were intersected with a list of aging-associated genes to identify aging-associated DEGs. Functional analyses were performed to elucidate the functions of these genes. Then, public datasets and online tools were evaluated to explore potential miRNAs and circRNAs. Overall, a regulatory circRNA-miRNA-mRNA network associated with aging was established to study its potential influence on AMD pathology.

Male C57BL/6J male mice (6–8 weeks old) were purchased from the Shanghai SLAC Laboratory Animal Company. Laser photocoagulation was performed as described previously [14]. Briefly, we anesthetized the mice with intraperitoneal injection of 1.5% sodium pentobarbital (100 μL/20 g) and dilated the pupils with tropicamide. A cover glass was applied as a contact lens to flatten the cornea, and the laser beam was focused on the RPE layer. Each eye was created with four laser spots between the major vessels in a relatively symmetric manner using a 532-nm laser with a power of 120 mW and duration of 100 ms (Visulas 532S; Carl Zeiss Meditec, Dublin, Ireland) around the optic nerve via a slip lamp delivery system. After laser induction, a drop of artificial tear solution was applied to keep the eye hydrated and prevent potential cataract development. Only laser spots showing white bubbles without hemorrhage were included in the study. Three or 7 days after laser treatment, the mice were sacrificed, and the RPE/choroid complexes were dissected around the corneal limbus for RNA samples. The groups were as follows: control group (2 months, n = 6), CNV-3d (2 months, n = 3), CNV-7d (2 months, n = 6), and CNV-7d (6 months, n = 3).

The total RNA of RPE/choroid complexes was isolated according to the RNAsimple Total Kit protocol (Tiangen Biotech, Beijing, China) following the manufacturer’s instructions and further quantified using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). cDNA was synthesized according to the RT Master Mix protocol (Takara Bio Inc., Dalian, China).

Raw data from the public datasets utilized in our study were obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/). The GSE146887 [15], based on the GPL18573 platform, is an mRNA profile consisting of RPE/choroid complexes from four healthy controls and four NVAMD patients. GSE131646 [16] is an miRNA profile with three normal controls and three RPE/choroid complexes from laser-induced CNV mice. This dataset was established on the GPL19057 platform (Table 1).

Raw data was log₂ transformed before further analysis. The “limma” package in R software (version 4.2.0, New Zealand) was used to identify DEGs in sequencing data [17]. During the analysis in this study, the thresholds of differentially expressed mRNAs/miRNAs were set as p < 0.05 and log₂|fold change| > 1. DEGs were further visualized with both heatmaps and volcano plots with “heatmap” and “ggplot” packages in R software.

The aging-related genes were mined from the Human Aging Genomic Resources (https://genomics.senescence.info/). Overlapping was performed between the DEGs in microarray analysis and the aging-related gene list with a Venn plot (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Functional analyses of these aging-related DEGs were performed with the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool (https://david.ncifcrf.gov/) and shown as bar graphs. The Gene Ontology (GO) analysis consisted of three parts: molecular function, cellular component, and biological process. Only the terms with p < 0.05 were included as significantly enriched.

The Search Tool for the Retrieval of Interacting Genes/Proteins website (https://cn.string-db.org/) was used to analyze and visualize the protein-protein interaction (PPI) association of these aging-related DEGs. This relationship was further exported to Cytoscape software (version 3.8.2) and analyzed with the Cytohubba plug-in unit, by which all the genes were ranked according to the maximal clique centrality algorithm.

The miRNet database (https://www.mirnet.ca/) was utilized to predict the miRNAs with potential binding positions with the identified DEGs. Considering that an miRNA can target multiple genes, we further overlapped these potential miRNAs with those identified in GSE131646 to narrow down the range. A Venn plot was used to visualize the 12 common miRNAs. Then, a miRNA-mRNA regulatory network was established to describe the interactions between them via Cytoscape software.

We used the starBase online tool (https://starbase.sysu.edu.cn/) to predict circRNAs with potential relationships with these miRNAs [18]. To narrow the range, only circRNAs overlapping with at least five miRNAs were included. Similarly, a circRNA-miRNA regulatory network was established to describe the interactions.

By inputting and merging the miRNA-mRNA and circRNA-miRNA networks in the Cytoscape software, we constructed and visualized a circRNA-miRNA-mRNA network.

The RNAs from control group (2 months, n = 6), CNV-3d (2 months, n = 3), CNV-7d (2 months, n = 6), and CNV-7d (6 months, n = 3) were reverse-transcribed to cDNAs by HiScript III RT SuperMix (Vazyme, China) with random hexamers as the primer. The first-strand cDNAs were quantified with a SYBR green master kit (Vazyme, China) on a real-time PCR machine ViiA^TM 7 (Applied Biosystems) using the standard curve and relative quantitation method with GAPDH as a control. Data are presented as means ± SD.

R software (version 4.2.0, New Zealand) was used for statistical analysis. The two-tailed Student’s t-test was used to compare the statistical significance between two groups, and a p value <0.05 was considered statistically significant.

With a cutoff of p < 0.05 and log₂|fold change| > 1, 1,175 genes were identified as differentially expressed (402 downregulated and 773 upregulated) in GSE146887. We visualized the DEGs using a heatmap and volcano plot (Fig. 1a, b). Human Aging Genomic Resources (HAGR) is a public database with a detailed profile of known aging-related genes. By overlapping the HAGR gene list with DEGs in sequencing, 23 aging-related DEGs were identified with a Venn diagram (Fig. 1c). We exhibited the expression of these 23 DEGs in the sequencing results via a heatmap and volcano plot (Fig. 2a, b). Figure 2c shows the expressions of these 23 genes in detail. ATR, ELN, LRP2, and PPM1D were downregulated, whereas the other 19 genes were downregulated (Table 2).

A PPI network was constructed via the Search Tool for the Retrieval of Interacting Genes/Proteins database (Fig. 3a). A total of 67 edges existed between these 23 aging-related DEGs, with an average node degree of 5.83 and a PPI enrichment p value of 1.11e−16. These were imported into the Cytoscape software and analyzed with the Cytohubba plug-in unit (Fig. 3b). Subsequently, the genes were ranked according to the gene ontology algorithm. AKT1 ranked the highest, followed by JUN and APOE (Fig. 3c).

Then, we utilized the DAVID tool to carry out GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for 23 aging-related DEGs (Fig. 4, online suppl. Table 1; for all online suppl. material, see https://doi.org/10.1159/000531287). Chagas disease was the most significantly enriched pathway in KEGG analysis, with a 25.3-fold enrichment. Other enriched diseases included colorectal cancer, renal cell carcinoma, and chronic myeloid leukemia. Significant pathways included the TNF, IL-17, and p53 signaling pathways (Fig. 4a). When considering reactome pathways, the molecules associated with elastic fibers were the most enriched items. The following items were the Toll-like receptor 5 (TLR5) cascade, Toll-like receptor 10 (TLR10) cascade, MyD88 cascade initiated on the plasma membrane, and TRAF6-mediated induction of NFkB and MAP kinases upon TLR7/8 or 9 activation. Senescence-associated secretory phenotype was also clustered in reactome analysis (Fig. 4b). In the GO analysis, the most enriched biological function was the positive regulation of cellular protein metabolic processes. The terms significantly involved included the response to muscle stretch, positive regulation of protein metabolic processes, positive regulation of fibroblast proliferation, and cellular response to reactive oxygen species. Aging was also enriched in biological functions, with a fold enrichment of 38.8 (Fig. 4c).

The miRNet online tool was utilized to predict the potential miRNAs targeting these aging-related genes, and as many as 867 miRNAs with 3,112 edges were obtained. To narrow that down, we further mined differentially expressed miRNAs in the public GSE131646 dataset, based on laser-induced NVAMD models. With a threshold of p < 0.05 and log₂|fold change| > 1, 52 miRNAs (four upregulated and 48 downregulated) were identified in this dataset. Among them, 12 miRNAs overlapped with those predicted in miRNet. A Venn plot was created to present these commonly differentially expressed miRNAs (Fig. 5a, b). A volcano plot (Fig. 5c) and heatmap (Fig. 5d) were further used to visualize their expressions. As the sponging relationship between miRNAs and mRNAs requires an opposite expression in diseases, we excluded the miRNA-mRNA pairs with the same trend in NVAMD, and nine miRNAs were finally included (Table 3). Subsequently, we constructed an aging-related miRNA-mRNA network using Cytoscape (Fig. 6). Among them, miR-129-2-3p had the highest score of 11, followed by miR-183-5p and miR-182-5p.

We further explored the potential circRNAs for these nine miRNAs via the starBase tool, and 1,751 circRNAs were predicted. To construct the most promising network, only 13 circRNAs with binding sites for at least five miRNAs were included for analysis, namely, circ-CLIP2, circ-TCONS_l2_00025633, circ-STAG3L3, circ-GTF2I, circ-POM121C, circ-BAZ1B, circ-PTK2, circ-SACS, circ-HUWE1, circ-NCKAP1, circ-HECTD1, circ-DYNC1H1, and circ-TOP2A. Next, the miRNA-mRNA and circRNA-miRNA networks were merged to construct a circRNA-miRNA-mRNA network, containing 281 pairs of circRNA-miRNA and 30 pairs of miRNA-mRNA (Fig. 7).

To validate the sequencing results of the identified genes, we conducted real-time quantitative PCR (RT-qPCR) using laser-induced NVAMD mouse models. Nine of the 15 genes were significantly differentially expressed in laser-induced CNV mouse models compared to control samples (Fig. 8a). Among them, eight genes presented consistent trends with sequencing, namely, FOS, JUN, GRN, NFKBIA, CEBPB, IGFBP2, PML, and TGFB1. When laser treatment was applied to both young and old mice, five genes including AKT1, NFKBIA, VCP, DGAT1, and EFEMP1 were even more upregulated in older mice than in younger ones (Fig. 8b).

To our knowledge, this is the first study to construct an aging-related circRNA-miRNA-mRNA network to investigate AMD etiology. Public microarray datasets in the GEO database were used to determine DEGs in NVAMD. After intersecting with the aging-related gene list from the HAGR database, 23 genes were identified. Nine miRNAs potentially targeting these genes were further identified with both the miRNet and public datasets. Subsequently, circRNAs with binding positions for all at least 5 miRNAs were predicted using the starBase tool, and 13 circRNAs were mined. Thereby, an aging-related circRNA-miRNA-mRNA network was established to explore the promising regulatory targets and pathways in NVAMD.

As one of the leading causes of irreversible visual impairment, NVAMD has resulted in a severe financial burden and time constraints on the geriatric population [19]. CNV is a critical pathological event of NVAMD. Angiogenic growth factors, such as VEGF, upregulate and further promote pathological angiogenesis. Therefore, anti-VEGF treatment was introduced and became the first-line therapy. However, as many as 45% of NVAMD patients do not respond to anti-VEGF therapy, according to a previous report [20]. Thus, the development of alternative therapies targeting other disease-associated mechanisms is warranted.

As its name implies, the prevalence of AMD increases with age, and cell senescence generally occurs with aging [21]. Multiple etiologies, such as inflammation, oxidative stress, and metabolic dysregulation, promote AMD progression. These factors also have a cause-and-effect relationship with aging, resulting in macular degeneration. For example, oxidative stress extensively attacks DNA and many other cellular organelles. Hydrogen peroxide increases ROS generation, DNA damage, and further cell senescence in RPE cells [22]. Furthermore, senescent cells secrete many chemokines, cytokines, and proteases, such as IL-1beta, IL-8, and matrix metalloproteinases (MMPs), further inducing senescence-associated secretory phenotype. The secreted inflammatory factors result in a low-grade and chronic inflammatory environment, which is a great elicitation for AMD pathology. Proteases, such as MMP-9, can destroy the RPE-formed outer blood-retinal barrier [23].

In order to elucidate the role of aging in NVAMD comprehensively, we screened public sequencing data from NVAMD patients and identified hundreds of DEGs. HAGR is a public database with a detailed profile of known aging-related genes. After the intersection, 23 genes were commonly screened; their expressions were both significantly altered in NVAMD and associated with aging. Additionally, RT-qPCR was conducted to validate their expression in mouse models. CDKN1A, C1QA, GRN, TGFB1, and CEBPB presented significant changes among them. CDKN1A (also known as p21) is an important cyclin-dependent kinase inhibitor, which is known to generate two transcripts. CDKN1A transcript variant 2 is identified as a more sensitive marker of aging and cellular senescence than total p21 in mice [24]. In p21-deficient mice, inflammatory responses and oxidative stress mediated by cigarette smoke are attenuated [25]. As a major constituent of complement C1q, C1QA gene disruption could relieve aging-associated impairments [26]. C1Q stimulates tube formation in vitro as well as vessel sprouting ex vivo. In addition, C1QA−/− mice present defective angiogenesis in wound healing in vivo [27]. Expression of C1QA detected by RNAscope is conspicuous in atrophy areas of AMD patient samples [28]. In malignant melanoma, C1QA, CSNK1E, and SOD-2 form an aging risk model to predict prognosis [29]. GRN gene variants are widely studied in aging phenotypes, including hippocampal sclerosis and the human cerebral cortex [30, 31]. GRN depletion exacerbates the inflammatory response following acute injury. Aged GRN−/− mice present robust neuroinflammation [32]. The function of TGFB in NVAMD remains controversial; it promotes the expression of VEGF and induces the fibrosis process secondary to NVAMD. However, some studies found TGFB to be anti-angiogenic [33, 34]. The expression of TGFB in aqueous samples of NVAMD showed a decreasing trend compared to that in healthy controls [35]. Suppression of CEBPB promotes angiogenesis and inhibits inflammation in ischemic stroke [36].

As noncoding RNAs are important regulators in multiple pathological processes, the public dataset GSE131646 and miRNet online tool were combined to search for mutual miRNAs sponging these aging-related genes, namely, miR-129-2-3p, miR-183-5p, miR-182-5p, miR-1224-5p, miR-383-5p, miR-155-3p, miR-770-5p, miR-129-1-3p, and miR-138-5p. Numerous miRNAs have been revealed to influence senescence by regulating aging-related genes. Among the miRNAs we identified, miR-183 and miR-138 were reported before [37]. The miR-183 cluster is highly expressed in the retina. By generating a mouse model in which the miR-183 cluster is inactivated, researchers find that miR-183 cluster inactivation results in retinal function defects and susceptibility to light damage [38]. miR-138-5p aggravates the decrease in aged osteoblast differentiation and aging-related bone loss by targeting downstream genes [39]. Furthermore, an miR-138-5p agonist inhibits ROS accumulation and ferroptosis in high glucose-treated RPE cells [40]. miR-155-3p is an inflammation-related miRNA, which upregulates TNF-α, IL-6, and IL-8 in IL-1-β/IFN-γ-activated astrocytes [41]. miR-383-5p aggravates lactate dehydrogenase-mediated oxidative stress and inflammation in spinal cord injury by regulating LDHA [42]. miR-1224-5p is upregulated in high glucose-treated human retinal microvascular endothelial cells. Its upregulation promotes cell proliferation, migration, and tube formation by targeting IRS-1 [43]. miR-182-5p regulates pathological angiogenesis in several cancer types as well as retinal neovascularization [44, 45]. Overexpression of miR-182-5p suppresses the expression of angiogenic factors significantly in both retinal neovascularization mouse and cell models [45]. On the contrary, exosomal miR-182-5p promotes tumor angiogenesis by accumulating VEGFR [46]. An miR-129-2-3p inhibitor promotes tube formation and VEGFA expression in breast cancer cells [47].

CircRNAs are a special type of noncoding RNAs with multiple binding motifs for miRNAs; therefore, they usually function by sponging miRNAs. We applied the starBase tool to predict the potential circRNAs with binding positions with these miRNAs mentioned above, and 13 circRNAs were screened. Circ-CLIP2 promotes glioma progression by sponging miR-195-5p and further HMGB3 [48]. Circ-GTF2I regulates miR-590-5p to induce inflammation and apoptosis factors, such as IL-6, TNF-α, Bax, and Bcl-2, in myocardial ischemia-reperfusion injury [49]. Circ-PTK2 has been studied in various diseases, including multiple myeloma, acute myeloid leukemia, and ovarian cancer. Generally, it promotes cell proliferation and migration [50, 51]. Circ-HUWE1 has been studied in aging-related Alzheimer’s disease. Its absence relieves amyloid-β-Induced cell viability degradation, apoptosis, and inflammatory responses in an AD cell model [52]. Circ-HECTD1 was studied in vascular diseases, including cerebral ischemia injury and cerebral infarction [53, 54]. Circ-TOP2A is also studied in glioma, exacerbating disease progression by sponging miR-345 [55]. As circRNAs were not widely studied until 2003, their roles are not yet well elucidated. In addition, the role of these circRNAs in aging and NVAMD remains unclear, requiring further study.

Our study has some limitations. First, only a few datasets based on NVAMD patients are available to the public currently. Therefore, we did not employ multiple datasets to validate our results. Second, we validated the expression of these genes in laser-induced CNV mouse models. This experimental model is the most widely used one for NVAMD, which recapitulates subtretinal angiogenesis. However, the mouse models available now failed to imitate the long-term disease progression of NVAMD. Therefore, further experimental investigations with more appropriate models would be meaningful to extend our study findings.

In this study, we identified and validated aging-related genes in NVAMD. Furthermore, an aging-related circRNA-miRNA-mRNA network was constructed for NVAMD etiology for the first time. This network provides novel insights into the noncoding RNAs and aging mechanisms of NVAMD, which could also be promising targets and biomarkers for further therapeutic exploration.

Animal experiments were approved by the Animal Care and Use Committee at Shanghai Jiao Tong University, and the procedures were conducted according to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research (permit number: 2019AW055).

The authors declare that they have no competing interests.

This study was supported by the National Natural Science Foundation of China (82171076) and Shanghai Hospital Development Center (SHDC2020CR5014).

Jiali Wu, Junran Sun, and Xiaodong Sun established the concept of this research. Jiali Wu drafted the manuscript. Yuxin Jiang performed the experiments. Junran Sun and Xiaodong Sun reviewed the manuscript. All authors have read and approved the manuscript.

Jiali Wu and Yuxin Jiang contributed equally to this work.

All data generated or analyzed during this study are included in this published article. Further inquiries can be directed to the corresponding author.