Identification of differentially expressed and methylated genes and construction of a co-expression network in age-related macular degeneration

Introduction

With the aging population, age-related macular degeneration (AMD) is the main cause of blindness worldwide for people over 50 years old (1). The most frequent AMD phenotype exhibits a feature of an increased number and diameter of extracellular retina deposits called drusen, pigmentary irregularity, progressive atrophy of the retinal pigmented epithelium (RPE) and retina, as well as a graded loss in visual acuity (2). It is estimated that by 2040, the number of AMD patients worldwide will increase to 288 million (3). However, the current treatment of AMD is limited. A thorough understanding of its molecular mechanisms may help promote the progression of effective therapies against AMD.

Various risk factors such as genetic factors as well as epigenetic modification are the cause of AMD (4). With the advance of microarray or RNA sequencing technology, gene expression profiling has been used to identify genes and pathways related to the pathogenesis of AMD, which helps to uncover its underlying mechanism. More than 50 genetic susceptibility loci have been identified, especially the CFH and ARMS2 genes (5). In addition to genetic risk factors, environmental factors such as smoking contribute to the progression of AMD. It has been thought that environmental factors could accelerate the accumulation of epigenetic changes over time (6). Nevertheless, the exact pathogenesis of AMD caused by DNA methylation is still unclear.

Epigenetic modification is a key factor of AMD (7), which is composed of DNA methylation, histone modification, and chromatin remodeling. The DNA methylation may mediate gene transcription and cell differentiation in association with gene silencing while demethylation induces gene transcription in AMD (8). Increasing evidence suggests that DNA methylation is involved in complex biological processes such as aging, inflammation, and apoptosis in AMD (9). Extensive research has been conducted on epigenetics for AMD. For example, a previous study constructed genome-wide DNA methylation profiles of retinal tissues and blood samples from AMD patients (10). Another study identified 2 novel DMGs, SKI and GTF2H4, for AMD using genome-wide methylation analysis (11). Besides, Xu et al. identified five hub genes might serve as aberrant methylation-based candidate biomarkers for AMD, including NOP56, EZR, IGF2, SLC2A1, CDKN1C (12). Weighted gene co-expression network analysis (WGCNA) has been widely used to explore modules and hub genes related to clinical features for various diseases (13,14). Nevertheless, little is known about the WGCNA co-expression network characteristics in AMD.

The Gene Expression Omnibus (GEO) database is a gene expression database created and maintained by the National Center for Biotechnology Information (NCBI). It was founded in 2000 and contains high-throughput gene expression data submitted by research institutions around the world. Based on the gene expression profiles, a systematic network can help identify driver genes or disease-related pathways. Herein, this study aimed to identify a subset of differentially methylated genes (DMGs) and differentially expressed genes (DEGs) that could regulate the biological processes and pathways involved in AMD and to construct an AMD-related co-expression network based on WGCNA algorithm. The expression of these hub genes was further verified by RT-qPCR. Our research could provide new targets for AMD in the future therapeutic intervention.

We present the following article in accordance with the STREGA reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-7043/rc).

Methods

Data acquisition

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Microarray expression (GSE29801 dataset) and DNA methylation (GSE102952 dataset) profiles were retrieved from the GEO (http://www.ncbi.nlm.nih.gov/geo/) database. The GSE29801 dataset was based on the GPL4133 platform, composed of 27 normal extramacular retina samples, 31 AMD extramacular retina samples, 28 normal macular retina samples, and 32 AMD macular retina samples (2). The GSE102952 dataset included 9 control retina tissues and 9 AMD retina tissues on the GPL13534 platform (12,15).

Data preprocessing

The microarray expression profile data were normalized by quartile using the normalize Between Arrays function of the limma package in R, which were then transformed into log2 to obtain the gene expression matrix of the samples (16,17). The Infinium Methylation 450K was used to profile. The standard index of DNA methylation at any specific CpG site is β = M/(M + U + α), where M and U are methylated and unmethylated signal intensities and α is an arbitrary offset (usually 100) intended to stabilize β values where fluorescent intensities are low. Absolute methylation values (β-values) were extracted with the GenomeStudio™ software and then subjected to detailed analysis (18). The β value of Illumina 450K DNA methylation was normalized using the beta-function of the wateRmelon package in R (19).

Principal component analysis (PCA)

PCA of samples was presented for removal of outlier samples, and to identify samples with high similarity. Samples with one-way analysis of variance (ANOVA) P≤0.05 were screened out. Furthermore, the correlation coefficients of gene expression levels between samples were calculated. The correlation coefficients can reflect the degree of similarity between samples.

Differential expression and methylation analysis

Differential expression and methylation analysis was presented using the limma package in R (https://www.rdocumentation.org/packages/limma/versions/3.28.14). The cutoff value was set as |fold change (FC)| >1.2 and P<0.05. Genes with FC >1.2 and P<0.05 were highly expressed in AMD while genes with FC <−1.2 and P<0.05 were lowly expressed in AMD. Meanwhile, genes with FC >1.2 and P<0.05 were hypermethylated in AMD, while genes with FC <−1.2 and P<0.05 were hypomethylated in AMD. Then, hypermethylated and down-regulated and hypomethylated and up-regulated genes were separately identified by intersection of genes with high methylation and low expression as well as genes with low methylation and high expression.

Protein-protein interaction (PPI) network

The interactions between proteins were predicted using the Search Tool for the Retrieval of Interacting Genes (STRING; version 11.0) online database (20), with the threshold of required confidence (combined score) >0.7. A PPI network was visualized using the Cytoscape software (version 3.5.1; https://cytoscape.org/) (21). Hub genes with the highest degree were identified.

Function annotation analysis

Function annotation analysis of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes or genes in the AMD-related modules was carried out using the Gene Ontology (GO) database (http://www.geneontology.org) (22) and Kyoto Encyclopedia of Genes and Genomes (KEGG) PATHWAY database (http://www.genome.ad.jp/kegg/) (23). The GO terms were composed of biological process (BP), molecular function (MF), and cellular component (CC). Fisher’s exact test was used to determine which specific functional items were most closely related to the list of the above genes. The P value corresponding to each term was calculated. The smaller the P-value, the greater the connection between the terms and the above genes, indicating that most of the above genes had the function corresponding to the terms.

External validation

Gene expression profiles of AMD and normal samples were downloaded from the E-MTAB-7183 dataset. The expression of hypomethylated and up-regulated hub genes and hypermethylated and down-regulated hub genes was externally validated in AMD and normal samples. The differences in their expression between AMD and normal samples were compared by Student’s t-test.

WGCNA

The WGCNA package in R was utilized to construct a co-expression network (24). Following construction of a Pearson correlation matrix, the weighted adjacency matrix was presented by combining the Pearson correlation between genes and soft thresholding parameters (β). The β value can amplify the correlation between genes by enhancing high correlation and weakening low correlation. In this study, β value was set to 4 to ensure a scale-free co-expression network. After the adjacency relation was converted into a topological overlap matrix (TOM), hierarchical clustering was carried out to identify modules within similar genes (25). Then, genes with similar expression patterns were clustered into the same module. Furthermore, the average linkage hierarchical clustering was presented in line with the TOM dissimilarity measure (13). The minimum number of genes was set to 30. Next, the dendrogram was constructed. Following evaluation of the dissimilarity of module eigengenes (MEs), a cut line was determined for the module tree, and similar modules were merged.

Construction of co-expression modules

Co-expression modules were constructed by calculating the correlation between MEs and clinical traits. The P value of each gene was converted into log10, which was defined as gene significance (GS). The median GS of all genes in each module was defined as module significance (MS). The module with high MS value was considered correlated with clinical traits.

Cell culture and treatment

Human ARPE‐19 cells (American Type Culture Collection; ATCC, Manassas, VA, USA) were grown in Dulbecco’s modified Eagle medium (DMEM) plus 10% fetal bovine serum (FBS), glutamine, and 50 U penicillin/50 mg streptomycin. To construct an AMD cell model, ARPE‐19 cells were maintained in a 95% N₂ and 5% CO₂ microenvironment. Normal ARPE‐19 cells were maintained in a 95% O₂ and 5% CO₂ microenvironment.

Enzyme-linked immunosorbent assay

The cell culture medium was collected. At 4 °C, the culture medium was centrifuged at 5,000 g for 15 min. The supernatant was then harvested. The levels of HIF-1α (ab227898; Abcam, Cambridge, MA, USA) and VEGF (ab222510; Abcam, USA) in the cell supernatant were detected according to the kit instructions.

Real-time quantitative polymerase chain reaction (qPCR)

Total RNA was extracted from ARPE‐19 cells by TRIzol reagent. The concentration, purity, and integrity of the RNA were detected by a micro-spectrophotometer. The complementary DNA (cDNA) was synthesized according to the instructions of the reverse transcription kit. Using cDNA as a template, by fluorescence qPCR detection kit, amplification was carried out on the PCR instrument. The reaction conditions were as follows: pre-denaturation at 94 °C for 3 min; denaturation at 94 °C for 10 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s, for a total of 38 cycles; after the end of the cycle, extension was performed at 72 °C for 5 min. The primer sequences were as follows: FSTL1: 5'-GAGGGCAAGAGTACACCA-3' (F), 5'-TACGGCATAGACGACAGC-3' (R); CHRDL1: 5'-CCGAGGACGAGGAAGAAGA-3' (F), 5'-AGTTGTCCCATTGTACTCG-3' (R); IFITM3: 5'-ATAGCATTCGCCTACTCCG-3' (F), 5'-ACAGACAAAGCCACTGACG-3' (R); GAPDH: 5'-CAAGGTCATCCATGACAACTTTG-3' (F), 5'-GTCCACCACCCTGTTGCTGTAG-3' (R). The relative expression of FSTL1, CHRDL1, and IFITM3 was quantified with 2^−ΔΔCt method, with glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as an internal control.

Western blot

Total protein was extracted from ARPE‐19 cells using lysis (Beyotime, Shanghai, China). Above protein was separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS‐PAGE; Invitrogen, Carlsbad, CA, USA) and transferred onto polyvinylidene fluoride (PVDF) membrane (Invitrogen, USA). Afterwards, the membrane was incubated with primary antibodies against FSTL1 (1/1,000; ab223287; Abcam, USA), CHRDL1 (1/1,000; ab103369; Abcam, USA), IFITM3 (1/1,000; ab109429; Abcam, USA), and GAPDH (1/1,000; ab8245; Abcam, USA), followed by being incubation with horseradish peroxidase (HRP)‐conjugated goat anti‐rabbit IgG secondary antibodies (1/10,000; ab7090; Abcam, USA). The protein bands were visualized using enhanced chemiluminescent (ECL) kit.

Statistical analysis

Statistical analysis was presented by R language and GraphPad Prism software (GraphPad Software, La Jolla, CA, USA). Data were displayed as the mean ± standard deviation. The difference between 2 groups was evaluated with Student’s t-test, while multiple comparisons were carried out by ANOVA. A P value <0.05 indicated statistical significance.

Results

DEGs of extramacular and macular AMD

Microarray expression profiles concerning 27 normal extramacular retina samples, 31 AMD extramacular retina samples, 28 normal macular retina samples, and 32 AMD macular retina samples from the GSE29801 dataset were used for this study. After PCA, there were no outlier samples (Figure 1A). There was a distinctly high correlation between different samples (Figure 1B). The expression profiles were normalized, as shown in Figure 1C,1D. Under the threshold of |FC| >1.2 and P<0.05, 896 genes were differentially expressed in AMD extramacular retina samples compared to normal extramacular retina samples (Figure 2A). Furthermore, 1,161 DEGs were screened between AMD macular retina tissues and normal macular retina tissues (Figure 2B). A volcano diagram showed that there were 490 up-regulated genes and 406 down-regulated genes between AMD extramacular retina samples and normal extramacular retina samples (Figure 2C). The top 10 DEG with the highest |FC| were visualized, including LOC100128239, SERPINC1, LOC729684, HSPH1, HSPA1A, CP, FGFBP2, RGS1, NPVF, and RBM3. Moreover, 766 genes were up-regulated and 395 genes were down-regulated in AMD macular retina tissues compared to normal macular retina tissues (Figure 2D). The top 10 DEG with the highest |FC| included MTHFD2, CHAC1, ATF3, LOC64521, GDEP, CP, C3, RGS1, SERPINA3, and NPVF.

Figure 1 Data filtering and normalization for GSE29801 dataset. (A) PCA of normal extramacular retina samples, AMD extramacular retina samples, normal macular retina samples and AMD macular retina samples. (B) Heatmap showing the correlation between different samples. The color depth is proportional to the correlation coefficient. (C,D) Pre-normalization and post-normalization of different samples. PCA, principal component analysis; AMD, age-related macular degeneration.

Figure 2 Screening DEGs for extramacular and macular AMD. Heat map showing all DEGs between AMD extramacular retina samples and normal extramacular retina samples (A), and all DEGs between macular retina tissues and normal macular retina tissues (B). Volcano plots visualizing 490 up-regulated genes and down-regulated 406 genes for AMD extramacular retina tissues (C), and 766 up-regulated genes and 395 down-regulated genes for AMD macular retina tissues (D). Red: up-regulation and blue down-regulation. DEGs, differentially expressed genes; AMD, age-related macular degeneration.

DMGs for extramacular and macular AMD

Methylation profile data including 9 control retina tissues and 9 AMD retina tissues were retrieved from the GSE102952 dataset. No outlier samples were detected according to PCA results (Figure 3A). A high correlation between different samples was found, as shown in the heat map (Figure 3B). After screening with the threshold of |FC| >1.2 and P<0.05, 1,345 hypermethylated and 1,086 hypomethylated genes were determined for AMD retina tissues (Figure 3C). Heat maps depicted the differences in methylation of all DMGs between AMD retina tissues and control retina tissues (Figure 3D).

Figure 3 DMGs for extramacular and macular AMD. (A) PCA of 9 control retina samples and 9 AMD retina samples. (B) Heat map showed the correlation between different samples. The color depth is proportional to the correlation coefficient. (C) Volcano plots depicted 1,345 hypermethylated and 1,086 hypomethylated genes between AMD retina samples and control retina samples. (D) Hierarchical clustering analysis visualized the methylation levels of all DMGs between the 2 groups. Red: hypomethylation and blue: hypomethylation. PCA, principal component analysis; AMD, age-related macular degeneration; DMGs, differentially methylated genes.

Identification of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes for AMD extramacular retina

After interaction, 16 hypermethylated and down-regulated genes (Figure 4A) and 15 hypomethylated and up-regulated genes (Figure 4B) were identified for AMD extramacular retina samples. A PPI network was constructed based on these overlapping genes. There were 4 nodes in the network, which could be considered as hub genes (LMNB2, EMD, HLA-A, and HLA-B) for AMD extramacular retina (Figure 4C). Functional enrichment analysis was then presented. As shown in Figure 5A, these genes were distinctly related with immune-related biological processes such as cytokine-mediated signaling pathway, regulation of T cell mediated cytotoxicity, and antigen processing and presentation of endogenous antigen. These genes were involved in cell components like MHC class I protein complex, transport vesicle, and phagocytic vesicle (Figure 5B). Furthermore, they had the molecular functions of peptide antigen binding, secondary active transmembrane transporter activity, cadherin binding, and anion transmembrane transporter activity (Figure 5C). The KEGG enrichment analysis results showed that these genes were involved in several crucial pathways such as antigen processing and presentation, endocytosis, and natural killer cell-mediated cytotoxicity (Figure 5D).

Figure 4 Identification of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes for AMD extramacular retina. (A,B) Venn diagrams showing hypermethylated and down-regulated genes and hypomethylated & up-regulated genes for AMD extramacular retina. (C) A PPI network based on overlapping genes. AMD, age-related macular degeneration; PPI, protein-protein interaction.

Figure 5 Functional enrichment analysis of hypermethylated and down-regulated and hypomethylated and up-regulated genes for AMD extramacular retina. (A) GO-BP; (B) GO-CC; (C) GO-MF; (D) KEGG. AMD, age-related macular degeneration; GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Identification of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes for AMD macular retina

A total of 13 hypermethylated and down-regulated genes and 31 hypomethylated and up-regulated genes were identified for AMD macular retina tissues (Figure 6A,6B). These genes were used for construction of the PPI network. As shown in Figure 6C, there were 11 hub genes (HLA-A, HLA-B, HLA-DRB1, IFITM3, SAT1, MAOB, CHRDL1, FSTL1, HSPA1A, AR, and YAP1) in the network. Consistent with hypermethylated and down-regulated and hypomethylated and up-regulated genes, GO-BP results showed that these genes were primarily enriched in immune-related biological processes (Figure 7A). As for GO-MF, they were prominently correlated with plasma membrane, recycling endosome membrane, luminal side of membrane, and so on (Figure 7B). Moreover, these genes possessed the functions of peptide antigen binding, transition metal ion transmembrane transporter activity, calcium ion transmembrane transporter activity, transcription corepressor activity, and phosphatidylinositol binding (Figure 7C). As displayed in Figure 7D, these genes mainly participated in antigen processing and presentation, allograft rejection, graft-versus-host disease, and the like.

Figure 6 Identification of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes for AMD macular retina. (A,B) Venn diagrams demonstrating hypermethylated and down-regulated genes and hypomethylated and up-regulated genes for AMD macular retina; (C) construction of a PPI network through overlapping genes. AMD, age-related macular degeneration; PPI, protein-protein interaction.

Figure 7 Functional enrichment analysis of hypermethylated and down-regulated and hypomethylated and up-regulated genes for AMD macular retina. (A) GO-BP; (B) GO-CC; (C) GO-MF; (D) KEGG. AMD, age-related macular degeneration; GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Validation of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes in AMD

The expression of hypermethylated and down-regulated and hypomethylated and up-regulated hub genes between AMD and normal samples was validated in the E-MTAB-7,183 dataset. Among them, ATP2C2 (log2FC =0.033, P=8.67e-03), CACNG6 (log2FC =0.119, P=1.84e-02), CNKSR2 (log2FC =0.038, P=1.22e-02), FSTL1 (log2FC =0.719, P=1.31e-02), MTM1 (log2FC =0.082, P=1.42e-02), MYO10 (log2FC =0.212, P=1.71e-03), RPH3AL (log2FC =0.045, P=6.74e-04), and SORCS2 (log2FC =0.319, P=2.97e-02) were highly expressed in AMD compared to normal samples (Figure 8A). Meanwhile, B4GALNT1 (log2FC =−0.697, P=1.94e-02), LACTB (log2FC =−0.02, P=1.63e-02), and SLC35E2 (log2FC =−0.263, P=1.36e-02) were more lowly expressed in AMD than normal samples (Figure 8B).

Figure 8 Validation of hypomethylated and up-regulated hub genes and hypermethylated and down-regulated in E-MTAB-7183 dataset. (A) Expression of hypomethylated and up-regulated genes in AMD and normal samples; (B) expression of hypermethylated and down-regulated in AMD and normal samples. AMD, age-related macular degeneration.

Construction of AMD-related co-expression modules

A total of 27 normal extramacular retina samples, 31 AMD extramacular retina samples, 28 normal macular retina samples, and 32 AMD macular retina samples from the GSE29801 dataset were used for co-expression analysis. The clinical traits of these samples are shown in Figure 9A. Soft threshold β value was set as 4 to ensure a scale-free network (R²=0.005; Figure 9B). A heat map depicted the high co-expressions between all genes (Figure 9C). A total of 16 merged modules were identified (Figure 9D). As shown in Figure 9E,9F, the turquoise module was in significant correlation with macular AMD (r=0.46 and P=2e-07), normal macular (r=0.56 and P=5e-11), and extramacular AMD (r=−0.54 and P=2e-10). Furthermore, blue module was distinctly correlated with macular AMD (r=−0.45 and P=4e-07) and extramacular AMD (r=0.47 and P=6e-08). The 2 co-expression modules were selected for further analysis.

Figure 9 Construction of AMD-related co-expression modules. (A) AMD sample dendrogram and trait heatmap; (B) determination of appropriate soft threshold β value; (C) network heatmap plot. The branches correspond to different modules. The color-coded modules are displayed below the tree diagram and in the color bar on the left. The co-expression is expressed by gradually saturated yellow and red. Modules are corresponding to blocks within highly co-expressed genes; (D) construction of clustering dendrogram; (E) module-trait relationships; (F) module eigengene dendrogram and network heatmap. Green is representative of negative correlation and red is representative of positive correlation. AMD, age-related macular degeneration.

Functional enrichment analysis of genes in the turquoise and blue modules

We performed functional enrichment analysis of genes in the turquoise and blue modules. As shown in Figure 10A, we found that genes in the turquoise module were primarily enriched several biological processes such as chromosome organization, multicellular organismal reproductive process, and spermatid development. For GO-CC, these genes were significantly correlated with nucleus, photoreceptor cell cilium, nuclear lumen, and so on (Figure 10B). Furthermore, these genes had the functions of DNA binding, DNA N-glycosylase activity, neurotransmitter transmembrane transporter activity, and so on (Figure 10C). For KEGG enrichment analysis results, several key KEGG pathways were enriched, such as phototransduction, base excision repair, and RNA polymerase (Figure 10D). For genes in the blue module, GO-BP results demonstrated that they were enriched in nervous system development, regulation of ion transmembrane transport, regulation of heart contraction, and so on (Figure 11A). Moreover, these genes could participate in various cell components such as axon, cation channel complex, and dendritic tree (Figure 11B). As shown in Figure 11C, they possessed the molecular functions of calcium ion binding, tubulin binding, small GTPase binding, and so on. The KEGG pathway enrichment analysis results suggested that these genes were involved in a few pathways like axon guidance, the CAMP signaling pathway, and metabolic pathway (Figure 11D).

Figure 10 Functional enrichment analysis of genes in the turquoise module. (A) GO-BP; (B) GO-CC; (C) GO-MF; (D) KEGG. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 11 Functional enrichment analysis of genes in the blue module. (A) GO-BP; (B) GO-CC; (C) GO-MF; (D) KEGG. GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Validation of hub genes in AMD cellular models

To construct an AMD cellular model, ARPE‐19 cells were maintained in a hypoxic state. Enzyme-linked immunosorbent assay (ELISA) was used for detecting HIF1α and VEGF expression in cell supernatant. Compared to normal ARPE‐19 cells, HIF1α and VEGF levels were significantly up-regulated in cell supernatant of hypoxia-induced ARPE‐19 cells (Figure 12A), indicating that the AMD cellular model was successfully established. The real-time qPCR (RT-qPCR) results showed that hub genes including CHRDL, FSTL1, and IFITM3 displayed significantly higher expression in hypoxia-induced ARPE‐19 cells compared to normal cells (Figure 12B). Meanwhile, greater up-regulation of CHRDL, FSTL1, and IFITM3 expression was found in hypoxia-induced ARPE‐19 cells than in normal cells (Figure 12C,12D).

Figure 12 Validation of hub genes in hypoxia-induced AMD cellular models. (A) ELISA for detecting the levels of HIF1α and VEGF in cell supernatant of normal and hypoxia-induced ARPE‐19 cells; (B) RT-qPCR for the expression of hub genes including CHRDL, FSTL1, and IFITM3 in normal and hypoxia-induced ARPE‐19 cells; (C,D) Western blot of the expression of CHRDL, FSTL1, and IFITM3 in normal and hypoxia-induced ARPE‐19 cells. *, P<0.05; ***, P<0.001; ****, P<0.0001. ELISA, enzyme-linked immunosorbent assay; RT-qPCR, real time quantitative polymerase chain reaction.

Discussion

In this study, we identified 16 hypermethylated and down-regulated genes and 15 hypomethylated and up-regulated genes for extramacular AMD and 13 hypermethylated and down-regulated genes and 31 hypomethylated and up-regulated genes for macular AMD. Our findings suggested that changes in DNA methylation may affect the differential expression of genes in AMD. In the PPI network, there were 4 hub genes (LMNB2, EMD, HLA-A, and HLA-B) for extramacular AMD and 11 hub genes (HLA-A, HLA-B, HLA-DRB1, IFITM3, SAT1, MAOB, CHRDL1, FSTL1, HSPA1A, AR, and YAP1) for macular AMD. These differentially expressed and methylated genes (DEMGs) could participate in various biological processes and pathways. Xu et al. have identified five hub genes might serve as aberrant methylation-based candidate biomarkers and construct a co-expression network for AMD (12). However, WGCNA co-expression network has not been used to explore modules and hub genes related to clinical features for AMD yet. This has been the first study to analyze the co-expression networks for AMD. What’s more, aberrantly-methylated DEGs have been identified in AMD in previous studies (12,15); however, these studies have lacked external validation. In this study, we validated the expression of DEMGs in the E-MTAB-7183 dataset. Among them, ATP2C2, CACNG6, CNKSR2, FSTL1, MTM1, MYO10, RPH3AL, and SORCS2 were highly expressed in AMD compared to normal samples. Meanwhile, B4GALNT1, LACTB, and SLC35E2 were lowly expressed in AMD than normal samples. A co-expression network was constructed and 2 AMD-related modules were screened out. Thus, our results offered novel insights into the molecular mechanisms for AMD and probed into several hub genes that could become underlying therapeutic markers for AMD in further research.

The DEMGs at the center of a PPI network could possess more important functions than other genes. Functional enrichment analysis results revealed that these DEMGs were involved in various immune-related biological processes and pathways such as cytokine-mediated signaling pathway, regulation of T cell mediated cytotoxicity, antigen processing and presentation of endogenous antigen, antigen processing and presentation, endocytosis, natural killer cell mediated cytotoxicity, and so on. Excessive activation of the immune system is important in the pathogenesis of AMD. Among 4 hub genes for extramacular AMD, LMNB2 has been verified to possess the function of promoting proliferation, migration, and invasion of liver cancer cells, regulated by miR-122 (26), and it is a key factor for RNA granules as well as axonal trafficking (27). The EMD gene has been confirmed as a risk factor for AMD (28). We found that HLA and HLA-B were hub genes both for extramacular and macular AMD. A cross-sectional observation clinical study found that the most frequent HLA haplotype iPS-RPE therapy may be an underlying alternative to AMD (29). Furthermore, HLA-B rs1055821 could correlate with AMD (30). Among 11 hub genes for macular AMD, except for HLA and HLA-B, HLA-DRB1, IFITM3, SAT1, MAOB, CHRDL1, FSTL1, HSPA1A, AR, and YAP1 are the causative factors of various diseases. For example, HLA-DRB1*03 is correlated to anti-carbamylated protein for rheumatoid arthritis (31). Allele polymorphism of HLA-DRB1 is associated with systemic sclerosis (32). The IFITM3 gene can promote bone metastasis of prostate cancer (33). Alternative splicing of SAT1 could suppress melanoma progression (34). The MAO-B gene has been detected as an underlying therapeutic target for Alzheimer’s disease (35). The CHRDL1 gene may participate in the progression of pulmonary thromboembolism by regulation of angiogenesis (36). Upregulation of FSTL1 enhances glioblastoma resistance to temozolomide (37). To validate the expression of hub genes CHRDL, FSTL1, and IFITM3 in AMD, we established a hypoxia-induced AMD cellular model. Our data confirmed that CHRDL, FSTL1, and IFITM3 were significantly up-regulated in hypoxia-induced ARPE‐19 cells compared to normal cells. It has been predicted that HSPA1A is correlated to poor prognosis of epithelial ovarian cancer (38). However, the functions of these hub genes in AMD are still unclear, and in-depth verification is warranted.

Using WGNCA, we conducted a total of 16 co-expression modules for AMD, of which 2 were in significant correlation with AMD. Sor far, this has been the first study to analyze the co-expression networks for AMD. We further explored the biological significance of these 2 gene modules. For the genes in the turquoise module, several key KEGG pathways were enriched like phototransduction, base excision repair, and RNA polymerase. Our KEGG pathway enrichment analysis results suggested that genes in the blue module were involved in a few pathways like axon guidance, CAMP signaling pathway, and metabolic pathway.

In conclusion, our study probed into several hub DEMGs and pathways in AMD, which could become underlying therapeutic markers for AMD. We also constructed 2 AMD-related gene co-expression modules, deepening the underlying concerning the molecular mechanisms of AMD. However, there are some limitations about this study. For example, it was better to choose a down-regulated expression of representative gene after hypoxia. Besides, the more specific experiments about the function of hub gene based on animal model or clinical samples should be performed in the further study.

Conclusions

In this study, we identified a series of DEMGs for extramacular and macular AMD, which could be involved in a variety of biological processes and pathways. Among them, several hub genes were identified, which could become therapeutic markers of AMD. Furthermore, a co-expression network was conducted for AMD. Thus, our research probed into several key DEMGs and pathways for AMD, for which more in-depth verification through clinical and experimental research is warranted.

Acknowledgments

Funding: This work was funded by Natural Science Foundation of Guangxi (2020GXNSFAA259054); The First Batch of High-Level Talent Scientific Research Projects of the Affiliated Hospital of Youjiang Medical University for Nationalities in 2019 (R20196340, Y20196304); and The Fund of Baise Scientific Research and Technology Development Plan (Baise 20184708).

Footnote

Reporting Checklist: The authors have completed the STREGA reporting checklist. Available at https://atm.amegroups.com/article/view/10.21037/atm-21-7043/rc

Data Sharing Statement: Available at https://atm.amegroups.com/article/view/10.21037/atm-21-7043/dss

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://atm.amegroups.com/article/view/10.21037/atm-21-7043/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Lim LS, Mitchell P, Seddon JM, et al. Age-related macular degeneration. Lancet 2012;379:1728-38. [Crossref] [PubMed]
2. Newman AM, Gallo NB, Hancox LS, et al. Systems-level analysis of age-related macular degeneration reveals global biomarkers and phenotype-specific functional networks. Genome Med 2012;4:16. [Crossref] [PubMed]
3. Handa JT, Bowes Rickman C, Dick AD, et al. A systems biology approach towards understanding and treating non-neovascular age-related macular degeneration. Nat Commun 2019;10:3347. [Crossref] [PubMed]
4. Gemenetzi M, Lotery AJ. Epigenetics in age-related macular degeneration: new discoveries and future perspectives. Cell Mol Life Sci 2020;77:807-18. [Crossref] [PubMed]
5. Mitchell P, Liew G, Gopinath B, et al. Age-related macular degeneration. Lancet 2018;392:1147-59. [Crossref] [PubMed]
6. Feinberg AP. Phenotypic plasticity and the epigenetics of human disease. Nature 2007;447:433-40. [Crossref] [PubMed]
7. Nashine S, Nesburn AB, Kuppermann BD, et al. Age-related macular degeneration (AMD) mitochondria modulate epigenetic mechanisms in retinal pigment epithelial cells. Exp Eye Res 2019;189:107701. [Crossref] [PubMed]
8. Caputo V, Strafella C, Termine A, et al. Epigenomic signatures in age-related macular degeneration: Focus on their role as disease modifiers and therapeutic targets. Eur J Ophthalmol 2021;31:2856-67. [Crossref] [PubMed]
9. Farkas MH, DeAngelis MM. Age-Related Macular Degeneration: From Epigenetics to Therapeutic Implications. Adv Exp Med Biol 2021;1256:221-35. [Crossref] [PubMed]
10. Oliver VF, Jaffe AE, Song J, et al. Differential DNA methylation identified in the blood and retina of AMD patients. Epigenetics 2015;10:698-707. [Crossref] [PubMed]
11. Porter LF, Saptarshi N, Fang Y, et al. Whole-genome methylation profiling of the retinal pigment epithelium of individuals with age-related macular degeneration reveals differential methylation of the SKI, GTF2H4, and TNXB genes. Clin Epigenetics 2019;11:6. [Crossref] [PubMed]
12. Xu Z, Ruan Z, Huang X, et al. Identification of aberrantly methylated differentially expressed genes in age-related macular degeneration. Medicine (Baltimore) 2019;98:e15083. [Crossref] [PubMed]
13. Xiao H, Chen P, Zeng G, et al. Three novel hub genes and their clinical significance in clear cell renal cell carcinoma. J Cancer 2019;10:6779-91. [Crossref] [PubMed]
14. Zeng D, He S, Ma C, et al. Co-Expression Network Analysis Revealed That the ATP5G1 Gene Is Associated With Major Depressive Disorder. Front Genet 2019;10:703. [Crossref] [PubMed]
15. Shen Y, Li M, Liu K, et al. Integrated bioinformatics analysis of aberrantly-methylated differentially-expressed genes and pathways in age-related macular degeneration. BMC Ophthalmol 2020;20:119. [Crossref] [PubMed]
16. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
17. Irizarry RA, Hobbs B, Collin F, et al. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 2003;4:249-64. [Crossref] [PubMed]
18. Dedeurwaerder S, Defrance M, Calonne E, et al. Evaluation of the Infinium Methylation 450K technology. Epigenomics 2011;3:771-84. [Crossref] [PubMed]
19. Pidsley R, Y, Wong CC, Volta M, et al. A data-driven approach to preprocessing Illumina 450K methylation array data. BMC Genomics 2013;14:293. [Crossref] [PubMed]
20. von Mering C, Huynen M, Jaeggi D, et al. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res 2003;31:258-61. [Crossref] [PubMed]
21. Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498-504. [Crossref] [PubMed]
22. Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 2000;25:25-9. [Crossref] [PubMed]
23. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. [Crossref] [PubMed]
24. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
25. Yip AM, Horvath S. Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinformatics 2007;8:22. [Crossref] [PubMed]
26. Li Y, Guo XB, Wang JS, et al. Function of fibroblast growth factor 2 in gastric cancer occurrence and prognosis. Mol Med Rep 2020;21:575-82. [PubMed]
27. Cosker KE, Fenstermacher SJ, Pazyra-Murphy MF, et al. The RNA-binding protein SFPQ orchestrates an RNA regulon to promote axon viability. Nat Neurosci 2016;19:690-6. [Crossref] [PubMed]
28. Ersoy L, Schick T, de Graft D, et al. Extramacular drusen are highly associated with age-related macular degeneration, but not with CFH and ARMS2 genotypes. Br J Ophthalmol 2016;100:1047-51. [Crossref] [PubMed]
29. Takagi S, Mandai M, Hirami Y, et al. Frequencies of human leukocyte antigen alleles and haplotypes among Japanese patients with age-related macular degeneration. Jpn J Ophthalmol 2018;62:568-75. [Crossref] [PubMed]
30. SanGiovanni JP, SanGiovanni PM, Sapieha P, et al. miRNAs, single nucleotide polymorphisms (SNPs) and age-related macular degeneration (AMD). Clin Chem Lab Med 2017;55:763-75. [Crossref] [PubMed]
31. Regueiro C, Rodriguez-Rodriguez L, Triguero-Martinez A, et al. Specific Association of HLA-DRB1*03 With Anti-Carbamylated Protein Antibodies in Patients With Rheumatoid Arthritis. Arthritis Rheumatol 2019;71:331-9. [Crossref] [PubMed]
32. Xu Y, Mo N, Jiang Z, et al. Human leukocyte antigen (HLA)-DRB1 allele polymorphisms and systemic sclerosis. Mod Rheumatol 2019;29:984-91. [Crossref] [PubMed]
33. Liu X, Chen L, Fan Y, et al. IFITM3 promotes bone metastasis of prostate cancer cells by mediating activation of the TGF-β signaling pathway. Cell Death Dis 2019;10:517. [Crossref] [PubMed]
34. Yang Q, Deng Y, Xu Y, et al. Knockdown of SSATX, an alternative splicing variant of the SAT1 gene, promotes melanoma progression. Gene 2019;716:144010. [Crossref] [PubMed]
35. Park JH, Ju YH, Choi JW, et al. Newly developed reversible MAO-B inhibitor circumvents the shortcomings of irreversible inhibitors in Alzheimer's disease. Sci Adv 2019;5:eaav0316. [Crossref] [PubMed]
36. Sun Y, Zhang X, Gao H, et al. Expression of microRNA-514a-5p and its biological function in experimental pulmonary thromboembolism. Am J Transl Res 2019;11:5514-30. [PubMed]
37. Nie E, Miao F, Jin X, et al. Fstl1/DIP2A/MGMT signaling pathway plays important roles in temozolomide resistance in glioblastoma. Oncogene 2019;38:2706-21. [Crossref] [PubMed]
38. De Andrade WP, Da Conceição Braga L, Gonçales NG, et al. HSPA1A, HSPA1L and TRAP1 heat shock genes may be associated with prognosis in ovarian epithelial cancer. Oncol Lett 2020;19:359-67. [PubMed]

(English Language Editor: J. Jones)