Network Integration Analysis and Immune Infiltration Analysis Reveal Potential Biomarkers for Primary Open-Angle Glaucoma

Primary open-angle glaucoma (POAG) is a progressive optic neuropathy and its damage to vision is irreversible. Therefore, early diagnosis assisted by biomarkers is essential. Although there were multiple researches on the identification of POAG biomarkers, few studies systematically revealed the transcriptome dysregulation mechanism of POAG from the perspective of pre- and post-transcription of genes. Here, we have collected multiple sets of POAG’s aqueous humor (AH) tissue transcription profiles covering long non-coding RNA (lncRNA), mRNA and mircoRNA (miRNA). Through differential expression analysis, we identified thousands of significant differentially expressed genes (DEGs) between the AH tissue of POAG and non-glaucoma. Further, the DEGs were used to construct a competing endogenous RNA (ceRNA) regulatory network and 1,653 qualified lncRNA-miRNA-mRNA regulatory units were identified. Two ceRNA regulatory subnets were identified based on the random walk algorithm and revealed to be involved in the regulation of multiple complex diseases. At the pre-transcriptional regulation level, a transcriptional regulatory network was constructed and three transcription factors (FOS, ATF4, and RELB) were identified to regulate the expression of multiple genes and participate in the regulation of T cells. Moreover, we revealed the immune desert status of AH tissue for POAG patients based on immune infiltration analysis and identified a specific AL590666.2-hsa−miR−339−5p-UROD axis can be used as a biomarker of POAG. Taken together, the identification of regulatory mechanisms and biomarkers will contribute to the individualized diagnosis and treatment for POAG.

Introduction

Glaucoma is the main cause of irreversible blindness, which includes several subtypes such as primary, secondary, angle-closure glaucoma and open-angle glaucoma (Weinreb and Khaw, 2004; Youngblood et al., 2019). Among them, primary open-angle glaucoma (POAG) is the most common. The clinical manifestations of POAG include optic nerve damage and loss of retinal ganglion cells, and high blood pressure and increased intraocular pressure are risk factors for POAG. Since the symptoms of POAG appear at a relatively late stage and the blindness caused by it is irreversible, early diagnosis is necessary (Weinreb et al., 2014). The identification of biomarkers is helpful for the early diagnosis of POAG patients (Kokotas et al., 2012). Although there were many studies on the identification of biomarkers of POAG, the limitations of the data have caused the limitations of the experimental results. For example, although Liu et al. found that hsa-miR-210-3p can be used as a biomarker of POAG from peripheral blood (Liu et al., 2019), the gene transcript also includes mRNA and long non-coding RNA (lncRNA) and did not reveal the pathogenesis of hsa-miR-210-3p.

With the progress of scientific research, the function of non-coding RNA has been unveiled. LncRNA and miRNA, as the two main types of non-coding RNA, have been shown to play an important role in the chromatin reprograming (Anastasiadou et al., 2018) and regulation of gene transcription through the competing endogenous RNAs (ceRNA) regulatory mechanism (Salmena et al., 2011). In the ceRNA network, mRNA and lncRNA act as a miRNA sponge to participate in ceRNA regulation determined by miRNA response elements (MREs) (Zhang et al., 2021). The ceRNA regulatory mechanism plays a role at the post-transcriptional level and is an important method to explain the dysregulation of transcript expression in diseases. For POAG, the construction of ceRNA regulatory network will assist in revealing its pathogenesis.

Over the past decade, the immune microenvironment has been a hot area of biological research, which includes immune infiltration, antigen presentation, immune cell exhaustion and immune cell communication. The immune microenvironment is composed of a variety of lymphocytes, such as T cells, B cells and macrophages, etc. Previous studies have shown that the neuroinflammatory response in POAG patients was thought to be caused by a defective immune response (Vernazza et al., 2020). For example, M1 polarization of macrophages enhances the antigen presentation ability and tissue inflammatory response (Yunna et al., 2020). Therefore, it is necessary to reveal the immune landscape of POAG to reveal its neuroinflammatory response mechanism.

In this study, we collected multiple sets of transcription profiles of aqueous humor (AH) tissues for POAG patients. Through the integrated analysis of the ceRNA competition network and the transcriptional regulatory network, we revealed the mechanism of the transcriptome dysregulation of POAG and the physiological functions that it affects. Immune infiltration analysis revealed the immune landscape of POAG. Additionally, potential biomarkers of POAG were identified based on machine learning algorithms.

Methods

Data Acquisition and Pre-Processing

The mRNA and long non-coding RNA (lncRNA) expression profiles of primary open-angle glaucoma (POAG) were downloaded from the Gene Expression Omnibus (GEO) database (accession number: GSE101727 (Xie et al., 2019), platform: GPL21827, Agilent-079487 Arraystar Human LncRNA miarray V4, Table 1). The raw annotation of GPL21827 only supported the sequence data format and not the gene symbol. Therefore, we mapped the sequences of GPL21287 probes to the human genome annotation file release GRCH37 in GENCODE (Harrow et al., 2012) using the R package “Rsubread” (Liao et al., 2019). Next, the average of standardized signal intensities was used to indicate mRNA/lncRNA expression intensity when multiple probes were mapped to the same mRNA/lncRNA. The miRNA expression profile was also obtained from the GEO database (accession number: GSE105269 (Drewry et al., 2018), platform: GPL24158, NanoString nCounter Human v3 miRNA Assay). The normalized miRNA expression matrix was used directly for the analyses (Table 1).

TABLE 1

TABLE 1. Expression datasets in the study.

Differential Expression Analysis of mRNA, lncRNA, and miRNA

Differentially expressed RNAs (mRNAs, lncRNAs, and miRNAs) were identified using the Linear Models for Miarray Data (Limma, v3.44.3) package in R (Ritchie et al., 2015). We considered the RNAs with |log₂FC| > 1.5 and p-value < 0.01 as the differentially expressed RNAs (DEmRNAs, DElncRNAs, and DEmiRNAs).

Constructing the POAG Associated Competing Endogenous RNA Network

To identify the POAG associated ceRNA relationships, we first recognized the candidate targets of DEmiRNAs based on the experimentally validated miRNA interaction relationships in lncACTdb v2.0 (http://www.bio-bigdata.net/LncACTdb/) (Wang et al., 2019), mirtarbase v2020 (http://miRTarBase.cuhk.edu.cn/) (Chou et al., 2018), and starbase v3.0 (https://starbase.sysu.edu.cn/) (Li et al., 2014). Next, according to the differentially expressed levels, the opposite changing trends between the expression levels of DEmiRNA-DEmRNA/DElncRNA pairs were retained in the AH (down-regulated miRNAs and up-regulated mRNAs/lncRNAs or up-regulated miRNAs and down-regulated mRNAs/lncRNAs). Furthermore, we calculated the rho between the expression levels of DElncRNAs and DEmRNAs. The raw p-values $\left(P_r\right)$ were adjusted by multiple hypotheses using a permutation method. For each mRNA, the expression value was held consistently, and 1,000 random lncRNAs were used to perform the same Spearman’s correlation test, generating a set of 1,000 permutation p-values $\left(P_p\right)$. Finally, an empirical p-value $\left(P_e\right)$ was corrected using the following formula (which introduces a pseudo-count of 1), i.e.

$P_e=\frac{\mathit{num}\left(P_p\le P_r\right)+1}{1001}(1)$

The mRNA-lncRNA pairs with the $rho>0.6$ and $P_e<0.01$ were treated as the co-expressed mRNA-lncRNA pairs. Finally, we constructed the ceRNA triplets relationships in POAG by integrating the miRNA-mRNA/lncRNA pairs and the co-expressed mRNA-lncRNA pairs (Wang et al., 2015). ceRNA network was visualized using the Cytoscape (Shannon et al., 2003).

Network-Based Prioritization of POAG-Related ceRNA Relationships Discovery

To identify the hub nodes in our ceRNA network, we employed the random walk with restart (Köhler et al., 2008). The POAG-related genes contained in the DisGeNet (Piñero et al., 2020) were considered as the seed genes. Performed random walk on the ceRNA network, with a restart probability of 0.7 using the function random walk in the R package RWOAG (Köhler et al., 2008). The nodes with top 30 visitation probabilities were treated as the hub nodes of the network. The ceRNA triplets consisting of hub nodes were considered as the critical ceRNAs relationships.

Construction of Transcriptional Regulatory Network

First, the immunosuppressive-related genes were collected from DisGeNET (Piñero et al., 2017) (http://www.disgenet.org) and HisgAtlas v1.0 (Liu et al., 2017) (http://biokb.ncpsb.org/HisgAtlas/). In addition, we searched for the keyword “immunosuppressive agents” in the Drugbank (Wishart et al., 2018) database (https://www.drugbank.ca/) and obtained 311 immunosuppressive-related drugs. In total, the 1,332 immunosuppressant-related genes were obtained from the above three databases. Next, the immunosuppressive-related genes in differentially expressed protein coding genes (PCGs) were extracted for the construction of transcriptional regulatory network. Moreover, the regulation data of transcription factors (TF) and target gene for human were downloaded from the TRRUST v2.0 (Han et al., 2018) (https://www.grnpedia.org/trrust/) and ORTI (Vafaee et al., 2016) databases (http://orti.sydney.edu.au/about.html). The TF target gene relationship pairs related to the immunosuppressive-related DEmRNAs were extracted. Further, the Spearman’s correlation coefficient (rho) between the genes of each pair was calculated and the cutoff of the p-value and rho were set to 0.05 and 0.5. Then, we constructed the TF-target network using Cytoscape software. We then analyzed the topological properties of the network and extracted the top 3 genes of degree as key drive factors.

Functional Enrichment Analysis

To annotated the potential biology functions of differentially expressed genes and ceRNA triplets, we performed functional enrichment analysis on the mRNAs using Metascape (http://metascape.org/gp/index.html) (Zhou et al., 2019). For the mRNA list, pathway and process enrichment analysis have been carried out with the following ontology sources: KEGG Pathway, GO Biological Processes, Reactome Gene Sets, and Canonical Pathways.

Immune Infiltration Analysis

The pre-processed expression matrix of PCGs for the GSE101727 series was used for immune infiltration analysis using CIBERSORT (Newman et al., 2015). CIBERSORT is a method to characterize the cellular composition of complex tissues from gene expression profiles.

Statistical Analysis

The ROC curves were performed using the R package pROC. The gene sets enrichment analysis using the Fisher’s exact test. All statistical analysis was performed using the R (v 3.6.2).

Results

Differential Expression Analysis Depicts the Transcriptional Features of POAG

In the regulation of gene expression, transcription is an initial step and one of the most critical steps (Prieto and McStay, 2005). To explore the changes in gene expression of POAG patients at the transcriptional level, the limma algorithm was used to identify genes that were significant differentially expressed in the AH of POAG compared to non-glaucoma. For GSE101727 series, 789 mRNAs were recognized to be down-regulated and 1,487 mRNAs were recognized to be up-regulated in AH (Figure 1A). In addition to considering the expression of protein-coding gene (PCG), the expression of non-coding gene that has been proven to play an important role in the activities of cell (Bridges et al., 2021) was also our focus. Further, there were 576 down-regulated and 614 up-regulated lncRNAs in AH (Figure 1B), which identified in GSE101727 series. Differentially expressed genes were used as features to identify non-glaucoma and POAG samples. We found that non-glaucoma samples and POAG samples can be distinguished and there are significant differences between the two groups (Figure 1C), indicating that the identified DEmRNAs and DElncRNAs can be used as the signature of POAG patients. In the ceRNA competition mechanism, miRNA is a crucial part (Smillie et al., 2018). Therefore, we specifically collected GSE105269 series data to identify DEmiRNA in AH of POAG. We found that 8 miRNAs (3 down-regulated and 5 up-regulated) were significantly differently expressed in AH of POAG (Figure 1D). Moreover, it is necessary to explore which physiological mechanisms these differentially expressed genes affect. The PCGs up-regulated in AH of POAG were used for functional enrichment analysis by Metascape tool. We found that the up-regulated PCGs are significantly enriched in protein synthesis and immune regulation (Figure 1E). Among them, the function of the top enrichment is the regulation of the expression of SLITs and ROBOs, which has been shown to be involved in the migration and positioning of neuronal precursor cells and the growth of neuronal axons (Tong et al., 2019). All these indicate that the neurons and immune microenvironment of AH in POAG patients have been changed.

FIGURE 1

FIGURE 1. The differentially expressed genes in AH of POAG. (A,B) The DEmRNAs and DElncRNAs between the two groups of AH tissue for POAG (GSE101727) and non-glaucoma are displayed by volcano map. The vertical dotted line is 1.5, and the horizontal dotted line is 2. These marked genes are identified in the following hub-subnets of ceRNA network. (C) The expression profile of differentially expressed genes (mRNA and lncRNA) is displayed with heat map. The column label shows the sample type. (D) The expression profile of the DEmiRNAs between the two groups of AH tissue for POAG (GSE105269) is displayed with heat map. The column label shows the sample type. (E) Functional enrichment results of DEmRNAs calculated by Metascape tool are displayed by bar plot. The darker the color, the more significant the enrichment function.

The ceRNA Network Reveal the Mechanism of Gene Expression Variations

The ceRNA regulatory mechanism plays an important role in the post-transcriptional regulation of genes. Using the ceRNA network to reveal the regulatory mechanisms of differentially expressed genes in AH tissue was conducive to the pathogenesis of POAG. Through the screening of genes involved in ceRNA regulation (see methods), we have identified 4 miRNAs that can bind to 13 lncRNAs and regulate their expression. Furthermore, the 4 miRNAs can regulate the expression of 333 mRNAs and then constitute 1,653 lncRNA-miRNA-mRNA regulatory units (Figure 2A). Through functional enrichment analysis, we found that the genes involved in ceRNA regulation are significantly enriched in cellular responses to stress, regulation of mRNA metabolic process and mRNA catabolic process (Figure 2B), indicating that the ceRNA mechanism in AH tissue of POAG could affect specific physiological mechanisms by regulating gene expression. Further, the hub nodes in the ceRNA network were identified by restarting the random walk algorithm (see methods). The ceRNA triad composed of hub nodes constitutes two ceRNA subnets. The ceRNA subnet_1 consists of 7 lncRNAs, 2 miRNAs and 10 mRNAs (Figure 2C). Among them, hsa-miR-21-5p as a ceRNA has been proven to play an important role in multiple complex diseases (Xiong et al., 2020; Jiang et al., 2020). The ceRNA subnet_2 consists of 5 lncRNAs, 1 miRNA and 5 mRNAs (Figure 2D). We found that FBN1 regulated by hsa-miR-339-5p is the causative gene of a variety of genetic diseases, including fibrinopathy and Marfan syndrome (Sakai et al., 2016). The genetic polymorphism of HSPA1B was closely related to the disorder of neuroregulation (Bosnjak Kuharic et al., 2020). Taken together, these suggest that the ceRNA regulatory mechanism plays an important role in AH tissue of POAG.

FIGURE 2

FIGURE 2. Construction of ceRNA regulatory network. (A) The POAG-related ceRNA triplets’ network in AH tissue. The yellow nodes represent the miRNAs. The red nodes represent the lncRNAs. The purple nodes represent the mRNAs. (B) Fuctional enrichment results of mRNAs in ceRNA networks are displayed by the bar plot. Three different colors represent three function sets. (C,D) The hub POAG-related ceRNA triplets’ sub-network in AH tissue. The yellow nodes represent the miRNAs. The red nodes represent the lncRNAs. The purple nodes represent the mRNAs.

Key Factors Driving the Progress of POAG

The TFs regulate the initiation and intensity of transcription of specific genes, which is an important driving factor in life activities (Lambert et al., 2018). To identify the driving factors that play an important role in POAG, a transcriptional regulatory network was constructed based on differentially expressed PCGs in GSE101727 series. By combining the previous research data and the correlation analysis of gene expression, 181 TF-target gene units were identified and the transcriptional regulatory network was constructed (Figure 3A). The network contained 73 TFs and 116 target genes. Further, functional enrichment analysis was used to explore the physiological mechanism involved in this transcriptional regulatory network. We found that the genes in this transcriptional regulatory network are significantly enriched in the regulation of myeloid cell differentiation and cell proliferation (Figure 3B), indicating that the specific expression of TF drives the expression of target genes and affects the immune microenvironment of POAG. Moreover, we identified the top 3 TFs (FOS, ATF4, and RELB) of degree as a key driver in the transcriptional regulatory network (Figures 3C‐E). The FOS and RELB were significantly down-regulated and ATF4 was significantly up-regulated (Figure 3D) in AH of POAG. Studies have shown that RELB plays a key role in the development of T cells and controls the proliferation of T cells, indicating that changes in RELB expression may be related to variations of the immune microenvironment for POAG (Zhou et al., 2020). Besides, we found that CDKN1A has the highest correlation with RELB at the transcript level (Figure 3F), and CDKN1A encodes cyclin which is regulated by kinase inhibitors (El-Deiry, 2016), indicating that RELB may control the proliferation of T cells by regulating the expression of CDKN1A. All these indicate that there are several driving factors that play an important role in the changes in the physiological mechanism of POAG.

FIGURE 3

FIGURE 3. Construction of transcriptional regulatory network. (A) The POAG-related transcriptional regulatory network in AH. The square represents the TF and target gene, the diamond represents TF, and the circle represents the target gene. Red color means up-regulation of gene expression, blue color means down-regulation. (B) Function enrichment results of mRNAs in transcriptional regulatory networks are displayed by the circle plot. The right side of the panel represents the function, and the left side represents the genes enriched in the function. (C–E) The expression of 3 TF (FOS, ATF4, and RELB) between AH tissue of POAG and non-glaucoma is shown by boxplot. (F) The correlation between the 3 TF (FOS, ATF4, and RELB) and its target genes. The color changes with the variation of the correlation coefficient.

Immune Infiltration Characteristics of POAG

The dynamics of the immune microenvironment is an important feature of the occurrence and development of diseases (Makowski et al., 2020). Identifying the immune characteristics is conducive to enriching the exploration of the pathogenesis of POAG. Therefore, the CIBERSORT tool was used to calculate the immune cell composition of each AH samples of POAG and non-glaucoma samples through the deconvolution algorithm. After preprocessing the immune cell fraction matrices, the consensus clustering algorithm was used to identify the distance between samples including AH of POAG and non-glaucoma samples. We found that POAG and non-glaucoma individuals can be distinguished by immune cell components and POAG individuals are in an immune desert state (Figure 4A). Further, the statistical test was used in the analysis of the difference between the immune cell components of the AH of POAG and non-glaucoma sample. From the perspective of the fold changes of immune cell components, the CD8⁺ T cell, CD4⁺ memory T cell, monocytes, macrophages, M1 and dendritic cell components of POAG individuals are significantly different from those of non-glaucoma individuals (Figure 4B). Besides, the Wilcoxon rank sum test was used to test the significance of differences in immune cell components between the two groups. We found that monocytes, γδ T cells, Tregs, CD8⁺ T cells and memory B cell components are significantly different between POAG and non-glaucoma individuals (Figure 4C). Monocytes, as a kind of myeloid cells, play an important role in presenting antigens in the organism and their fraction was significantly down-regulated in POAG, which may be an important sign of POAG patients. POAG is a neurodegenerative disease and neuroinflammation occurs during its pathogenesis (Weinreb and Khaw, 2004; Evangelho et al., 2019), which may be related to the lack of Treg. Taken together, these suggest that the AH of POAG is in an immune desert state and the significant down-regulation of specific immune cell components can be used as the marker of POAG.

FIGURE 4

FIGURE 4. The landscape of immune cell infiltration for POAG. (A) The fraction of immune cells in the AH tissue of POAG and non-glaucoma is shown by heatmap. The column label shows the sample type. (B) The fold change of the fraction of immune cells between the AH tissue of POAG and non-glaucoma. The two vertical dashed lines represent 1 and -1 respectively. (C) The fraction of immune cells in the AH tissue of POAG and non-glaucoma is shown by boxplot. Wilcoxon rank sum test is used to calculate the statistical difference between the two groups.

Biomarkers of POAG

The identification of biomarkers of POAG is helpful for its clinical diagnosis and treatment. Therefore, we collected important genes in the ceRNA regulatory subnet and drive factors identified above. For the two ceRNA regulatory subnet, the genes in each lncRNA-miRNA-mRNA unit were used as features to distinguish POAG from non-glaucoma individuals and the ROC curve was used to evaluate the stability of the feature. The AL590666.2-hsa−miR−339−5p-UROD axis was recognized to be able to stably distinguish between POAG and non-glaucoma individuals. Among them, UROD has the highest AUC value of 0.98 compared to 0.77 of AL590666.2 and 0.78 of hsa−miR−339−5p (Figures 5A–C). Uroporphyrinogen decarboxylase encoded by UROD was an important element in hemoglobin synthesis, which is significantly up-regulated in POAG (Figure 5D). Studies have shown that patients with POAG have red blood cell backlog and high plasma specific viscosity (Wang et al., 2004; Xu et al., 2020), indicating that the up-regulation of UROD may be an important cause of blood deformation in POAG patients. Further, we found that AL590666.2 and UROD have a strong correlation (Figure 5E), suggesting that AL590666.2 may be an important biomarker of POAG. For the top three TFs identified above, the AUC values of ATF4, FOS, and RELB were 0.91, 0.91, and 0.74, respectively (Figure 5F). All these suggest that these genes can be used as biomarkers of POAG for clinical diagnosis and treatment.

FIGURE 5

FIGURE 5. Machine learning was used to identify POAG biomarkers. (A–C) The potential of has-miR-339-5p, AL590666.2, and UROD in identifying POAG diseases is depicted by the ROC curve. The AUC values are calculated. (D) The expression of has-miR-339-5p, AL590666.2, and UROD in AH tissues of POAG and non-glaucoma. Wilcoxon rank sum test is used to calculate statistical differences. (E) The correlation between AL590666.2 and UROD in expression. (F) Same as in (A–C) but for the ATF4, FOS, and UROD.

Discussion

In this work, we have integrated multiple sets of transcript data (mRNA, lncRNA, miRNA) and revealed important functional subnets and driving factors of POAG through ceRNA competition network and transcriptional regulatory network analysis. Through statistical testing, thousands of genes (DEmRNA, DElncRNA, DEmiRNA) differentially expressed in POAG’s AH tissue have been identified. We have identified 1,653 lncRNA-miRNA-mRNA regulatory units and two functional subnets in AH tissue, which will help reveal the pathogenesis of POAG. Further, the transcriptional regulatory network was constructed based on differentially expressed genes and 3 TFs were recognized to play an important role in the transcriptome disorder of the POAG’s AH tissue. We have used the CIBERSORT tool and transcription profile of AH tissue to reveal the immune landscape of POAG. We found that the components of immune cells in AH tissue of POAG were globally down-regulated compared to non-glaucoma. Additionally, a ceRNA regulatory axis (AL590666.2-hsa−miR−339−5p-UROD) and 3 TFs (ATF4, FOS, and RELB) have been identified as potential biomarkers for POAG patients.

POAG is the most common form of glaucoma disease, which is a disease of the optic nervous system and causes irreversible blindness (Li et al., 2016). The identification of POAG’s biomarkers is the direction of the efforts of many researchers. For example, Liu et al. identified hsa-miR-210-3p in peripheral blood as a biomarker of POAG based on miRNA expression profile (Liu et al., 2019). However, polygene dysregulation and interaction were the inherent causes of POAG. We have revealed the regulatory relationship of dysregulated genes and the pathogenesis of POAG through multi-network integration analysis. The hsa-miR-21-5p and FBN1 were the key genes in the ceRNA regulatory subnet that we have identified, which play an important role in multiple complex diseases.

For our identified driving factor ATF4, it has been confirmed in previous studies that it can cause glaucoma by promoting ER client protein load (Kasetti et al., 2020) and regulating trabecular meshwork remodeling (Zhao et al., 2020). In addition to common transcriptome analysis, Liu et al. identified F-box protein (FBOX) and vaccinia-associated kinase 2 (VRK2) that may interact with tumor protein p53 (TP53) to regulate apoptosis and play a negative role in POAG from the perspective of genetic lineage (Liu et al., 2012). Among them, FBOX and TP53 were target genes of the key TF ATF4 and FOS identified in this work. Moreover, several potential biomarkers of POAG were revealed through integrated network analysis in this work.

Among identified protein-coding genes that are significantly differentially expressed in the AH tissue of POAG, RELB and CDKN1A were a pair of important transcriptional regulatory units. The RELB has been confirmed in previous studies to regulate the proliferation of T cells (Zhou et al., 2020) and the expression of its target gene CDKN1A was closely related to the cell cycle (El-Deiry, 2016). Therefore, the down-regulation of CDKN1A expression mediated by RENL may be an important reason for the decreased level of T cell infiltration. Red blood cell backlog and high plasma specific viscosity were an important physiological manifestation of POAG. This trait may be related to the up-regulation of UROD expression.

Conclusion

In conclusion, we integrated multi-network analysis to identify important functional subnets and driving factors, which will help advance the research on the pathogenesis of PAOG. Immune infiltration analysis and feature recognition reveal the immune desert state and the biomarkers of POAG. Taken together, our research provides theoretical guidance for the clinical diagnosis and treatment for POAG.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: GSE101727 and GSE105269 from the Gene Expression Omnibus (GEO) database.

Authors’ Contributions

HS, XC, and LW conceived and designed the experiments. LW and TY performed analysis and wrote the manuscript. TY and XZ collected the data. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (82004425, 82074497) and the Research Foundation of Heilongjiang University of Chinese Medicine (2019MS10).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors gratefully thank the GEO for providing data for this work.

References

Anastasiadou, E., Jacob, L. S., and Slack, F. J. (2018). Non-coding RNA Networks in Cancer. Nat. Rev. Cancer 18, 5–18. doi:10.1038/nrc.2017.99

PubMed Abstract | CrossRef Full Text | Google Scholar

Bosnjak Kuharic, D., Bozina, N., Ganoci, L., Makaric, P., Kekin, I., Prpic, N., et al. (2020). Association of HSPA1B Genotypes with Psychopathology and Neurocognition in Patients with the First Episode of Psychosis: a Longitudinal 18-month Follow-Up Study. Pharmacogenomics J. 20, 638–646. doi:10.1038/s41397-020-0150-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bridges, M. C., Daulagala, A. C., and Kourtidis, A. (2021). LNCcation: lncRNA Localization and Function. J. Cel Biol 220, e202009045. doi:10.1083/jcb.202009045

CrossRef Full Text | Google Scholar

Chou, C.-H., Shrestha, S., Yang, C.-D., Chang, N.-W., Lin, Y.-L., Liao, K.-W., et al. (2018). miRTarBase Update 2018: a Resource for Experimentally Validated microRNA-Target Interactions. Nucleic Acids Res. 46, D296–D302. doi:10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Drewry, M. D., Challa, P., Kuchtey, J. G., Navarro, I., Helwa, I., Hu, Y., et al. (2018). Differentially Expressed microRNAs in the Aqueous Humor of Patients with Exfoliation Glaucoma or Primary Open-Angle Glaucoma. Hum. Mol. Genet. 27, 1263–1275. doi:10.1093/hmg/ddy040

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Deiry, W. S. (2016). p21(WAF1) Mediates Cell-Cycle Inhibition, Relevant to Cancer Suppression and Therapy. Cancer Res. 76, 5189–5191. doi:10.1158/0008-5472.can-16-2055

PubMed Abstract | CrossRef Full Text | Google Scholar

Evangelho, K., Mogilevskaya, M., Losada-Barragan, M., and Vargas-Sanchez, J. K. (2019). Pathophysiology of Primary Open-Angle Glaucoma from a Neuroinflammatory and Neurotoxicity Perspective: a Review of the Literature. Int. Ophthalmol. 39, 259–271. doi:10.1007/s10792-017-0795-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, H., Cho, J.-W., Lee, S., Yun, A., Kim, H., Bae, D., et al. (2018). TRRUST V2: an Expanded Reference Database of Human and Mouse Transcriptional Regulatory Interactions. Nucleic Acids Res. 46, D380–D386. doi:10.1093/nar/gkx1013

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrow, J., Frankish, A., Gonzalez, J. M., Tapanari, E., Diekhans, M., Kokocinski, F., et al. (2012). GENCODE: the Reference Human Genome Annotation for the ENCODE Project. Genome Res. 22, 1760–1774. doi:10.1101/gr.135350.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, J., Wang, X., and Lu, J. (2020). PWRN1 Suppressed Cancer Cell Proliferation and Migration in Glioblastoma by Inversely Regulating Hsa-miR-21-5p. Cmar 12, 5313–5322. doi:10.2147/cmar.s250166

PubMed Abstract | CrossRef Full Text | Google Scholar

Kasetti, R. B., Patel, P. D., Maddineni, P., Patil, S., Kiehlbauch, C., Millar, J. C., et al. (2020). ATF4 Leads to Glaucoma by Promoting Protein Synthesis and ER Client Protein Load. Nat. Commun. 11, 5594. doi:10.1038/s41467-020-19352-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Köhler, S., Bauer, S., Horn, D., and Robinson, P. N. (2008). Walking the Interactome for Prioritization of Candidate Disease Genes. Am. J. Hum. Genet. 82, 949–958. doi:10.1016/j.ajhg.2008.02.013

CrossRef Full Text | Google Scholar

Kokotas, H., Kroupis, C., Chiras, D., Grigoriadou, M., Lamnissou, K., Petersen, M. B., et al. (2012). Biomarkers in Primary Open Angle Glaucoma. Clin. Chem. Lab. Med. 50, 2107–2119. doi:10.1515/cclm-2012-0048

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, S. A., Jolma, A., Campitelli, L. F., Das, P. K., Yin, Y., Albu, M., et al. (2018). The Human Transcription Factors. Cell 172, 650–665. doi:10.1016/j.cell.2018.01.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks from Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42, D92–D97. doi:10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Lindsley, K., Rouse, B., Hong, H., Shi, Q., Friedman, D. S., et al. (2016). Comparative Effectiveness of First-Line Medications for Primary Open-Angle Glaucoma. Ophthalmology 123, 129–140. doi:10.1016/j.ophtha.2015.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2019). The R Package Rsubread Is Easier, Faster, Cheaper and Better for Alignment and Quantification of RNA Sequencing Reads. Nucleic Acids Res. 47, e47. doi:10.1093/nar/gkz114

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Xie, L., Ye, J., Liu, Y., and He, X. (2012). Screening of Candidate Genes for Primary Open Angle Glaucoma. Mol. Vis. 18, 2119–2126.

PubMed Abstract | Google Scholar

Liu, Y., He, M., Wang, D., Diao, L., Liu, J., Tang, L., et al. (2017). HisgAtlas 1.0: A Human Immunosuppression Gene Database. Oxford: Database, 2017.

Google Scholar

Liu, Y., Wang, Y., Chen, Y., Fang, X., Wen, T., Xiao, M., et al. (2019). Discovery and Validation of Circulating Hsa-miR-210-3p as a Potential Biomarker for Primary Open-Angle Glaucoma. Invest. Ophthalmol. Vis. Sci. 60, 2925–2934. doi:10.1167/iovs.19-26663

PubMed Abstract | CrossRef Full Text | Google Scholar

Makowski, L., Chaib, M., and Rathmell, J. C. (2020). Immunometabolism: From Basic Mechanisms to Translation. Immunol. Rev. 295, 5–14. doi:10.1111/imr.12858

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust Enumeration of Cell Subsets from Tissue Expression Profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

Piñero, J., Ramírez-Anguita, J. M., Saüch-Pitarch, J., Ronzano, F., Centeno, E., Sanz, F., et al. (2020). The DisGeNET Knowledge Platform for Disease Genomics: 2019 Update. Nucleic Acids Res. 48, D845–D855. doi:10.1093/nar/gkz1021

PubMed Abstract | CrossRef Full Text | Google Scholar

Piñero, J., Bravo, À., Queralt-Rosinach, N., Gutiérrez-Sacristán, A., Deu-Pons, J., Centeno, E., et al. (2017). DisGeNET: a Comprehensive Platform Integrating Information on Human Disease-Associated Genes and Variants. Nucleic Acids Res. 45, D833–D839. doi:10.1093/nar/gkw943

PubMed Abstract | CrossRef Full Text | Google Scholar

Prieto, J.-L., and McStay, B. (2005). Nucleolar Biogenesis: the First Small Steps. Biochem. Soc. Trans. 33, 1441–1443. doi:10.1042/bst0331441

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakai, L. Y., Keene, D. R., Renard, M., and De Backer, J. (2016). FBN1: The Disease-Causing Gene for Marfan Syndrome and Other Genetic Disorders. Gene 591, 279–291. doi:10.1016/j.gene.2016.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA Hypothesis: the Rosetta Stone of a Hidden RNA Language. Cell 146, 353–358. doi:10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Smillie, C. L., Sirey, T., and Ponting, C. P. (2018). Complexities of post-transcriptional Regulation and the Modeling of ceRNA Crosstalk. Crit. Rev. Biochem. Mol. Biol. 53, 231–245. doi:10.1080/10409238.2018.1447542

PubMed Abstract | CrossRef Full Text | Google Scholar

Tong, M., Jun, T., Nie, Y., Hao, J., and Fan, D. (2019). The Role of the Slit/Robo Signaling Pathway. J. Cancer 10, 2694–2705. doi:10.7150/jca.31877

PubMed Abstract | CrossRef Full Text | Google Scholar

Vafaee, F., Krycer, J. R., Ma, X., Burykin, T., James, D. E., and Kuncic, Z. (2016). ORTI: An Open-Access Repository of Transcriptional Interactions for Interrogating Mammalian Gene Expression Data. PLoS One 11, e0164535. doi:10.1371/journal.pone.0164535

PubMed Abstract | CrossRef Full Text | Google Scholar

Vernazza, S., Tirendi, S., Bassi, A. M., Traverso, C. E., and Saccà, S. C. (2020). Neuroinflammation in Primary Open-Angle Glaucoma. J. Clin. Med. 9, 3172. doi:10.3390/jcm9103172

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G., Medeiros, F. A., Barshop, B. A., and Weinreb, R. N. (2004). Total Plasma Homocysteine and Primary Open-Angle Glaucoma. Am. J. Ophthalmol. 137, 401–406. doi:10.1016/j.ajo.2003.09.041

CrossRef Full Text | Google Scholar

Wang, P., Li, X., Gao, Y., Guo, Q., Wang, Y., Fang, Y., et al. (2019). LncACTdb 2.0: an Updated Database of Experimentally Supported ceRNA Interactions Curated from Low- and High-Throughput Experiments. Nucleic Acids Res. 47, D121–D127. doi:10.1093/nar/gky1144

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, P., Ning, S., Zhang, Y., Li, R., Ye, J., Zhao, Z., et al. (2015). Identification of lncRNA-Associated Competing Triplets Reveals Global Patterns and Prognostic Markers for Cancer. Nucleic Acids Res. 43, 3478–3489. doi:10.1093/nar/gkv233

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinreb, R. N., Aung, T., and Medeiros, F. A. (2014). The Pathophysiology and Treatment of Glaucoma. JAMA 311, 1901–1911. doi:10.1001/jama.2014.3192

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinreb, R. N., and Khaw, P. T. (2004). Primary Open-Angle Glaucoma. The Lancet 363, 1711–1720. doi:10.1016/s0140-6736(04)16257-0

CrossRef Full Text | Google Scholar

Wishart, D. S., Feunang, Y. D., Guo, A. C., Lo, E. J., Marcu, A., Grant, J. R., et al. (2018). DrugBank 5.0: a Major Update to the DrugBank Database for 2018. Nucleic Acids Res. 46, D1074–D1082. doi:10.1093/nar/gkx1037

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, L., Mao, M., Wang, C., Zhang, L., Pan, Z., Shi, J., et al. (2019). Potential Biomarkers for Primary Open-Angle Glaucoma Identified by Long Noncoding RNA Profiling in the Aqueous Humor. Am. J. Pathol. 189, 739–752. doi:10.1016/j.ajpath.2018.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, Y., Tang, Y., Fan, F., Zeng, Y., Li, C., Zhou, G., et al. (2020). Exosomal Hsa-miR-21-5p Derived from Growth Hormone-Secreting Pituitary Adenoma Promotes Abnormal Bone Formation in Acromegaly. Translational Res. 215, 1–16. doi:10.1016/j.trsl.2019.07.013

CrossRef Full Text | Google Scholar

Xu, M., Li, S., Zhu, J., Luo, D., Song, W., and Zhou, M. (2020). Plasma Lipid Levels and Risk of Primary Open Angle Glaucoma: a Genetic Study Using Mendelian Randomization. BMC Ophthalmol. 20, 390. doi:10.1186/s12886-020-01661-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Youngblood, H., Hauser, M. A., and Liu, Y. (2019). Update on the Genetics of Primary Open-Angle Glaucoma. Exp. Eye Res. 188, 107795. doi:10.1016/j.exer.2019.107795

PubMed Abstract | CrossRef Full Text | Google Scholar

Yunna, C., Mengru, H., Lei, W., and Weidong, C. (2020). Macrophage M1/M2 Polarization. Eur. J. Pharmacol. 877, 173090. doi:10.1016/j.ejphar.2020.173090

CrossRef Full Text | Google Scholar

Zhang, Y., Han, P., Guo, Q., Hao, Y., Qi, Y., Xin, M., et al. (2021). Oncogenic Landscape of Somatic Mutations Perturbing Pan-Cancer lncRNA-ceRNA Regulation. Front. Cel Dev. Biol. 9, 658346. doi:10.3389/fcell.2021.658346

CrossRef Full Text | Google Scholar

Zhao, Y., Zhu, H., Yang, Y., Ye, Y., Yao, Y., Huang, X., et al. (2020). AQP1 Suppression by ATF4 Triggers Trabecular Meshwork Tissue Remodelling in ET‐1‐induced POAG. J. Cel Mol Med 24, 3469–3480. doi:10.1111/jcmm.15032

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, S., Wu, W., Wang, Z., Wang, Z., Su, Q., Li, X., et al. (2020). RelB Regulates the Homeostatic Proliferation but Not the Function of Tregs. BMC Immunol. 21, 37. doi:10.1186/s12865-020-00366-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10, 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar