Integrative analysis of miRNA–mRNA network in high altitude retinopathy by bioinformatics analysis

High-altitude retinopathy (HAR) is one clinical entity of acute high-altitude illness (AAI), which also includes acute mountain sickness (AMS), high-altitude pulmonary edema (HAPE) and high-altitude cerebral edema (HACE) [1,2]. As retinopathy is caused by the inability to adapt to acute oxygen deficiency at high altitudes, HAR mainly manifests as leakage of the peripheral retinal vessels, retinal hemorrhages, tortuous retinal vessels, optic disc edema and macular edema [3,4]. With the development of the economy and tourism in the plateau area, the incidence of HAR is up to 79% and is increasing [5]. However, the pathogenesis of HAR remains unclear, and there is a lack of effective prevention and treatment measures.

Current studies suggest that hypoxia is the main reason for the pathophysiological change in high-altitude situations [6,7]. However, under the same hypoxic conditions, some high-altitude residents will remain in good condition, while others may have AAI [8]. Several genes, including ACE, EDN1, ACYP2, RTEL1, and VEGF, are related to hypoxia [9–11]. In hypoxic environments, transcription of various genes, such as endothelial PAS domain-containing protein 1 (EPAS1) and prolyl hydroxylase domain-containing protein 2 (PHD2), is initiated by hypoxia-related pathways [12]. Also, genetic susceptibility has been reported as one of the major determinants of HAPE and AMS. The oxidative stress-related genes CYBA and GSTP1, contributing to endothelial damage under hypobaric hypoxia, are potential candidate genes for HAPE [13–17]. However, there are a few genetic studies on HAR.

MicroRNAs (miRNAs) can inhibit mRNA translation or induce degradation of target mRNAs by binding to complementary sequences in target regions [18,19]. Huang et al. revealed that saliva miR-134-3p and miR-15b-5p could be used as noninvasive biomarkers to predict AMS individuals exposed to high altitude in advance [19]. In a hypoxic atmosphere, the expression levels of miR-16, 20b, 22, 206, and 17/92 are reduced, which can inhibit ion channels, increase pulmonary artery pressure, and cause vascular dysfunction and the loss of cell integrity, promoting the occurrence of HAPE, and miR-155, 23b, 26a can inhibit TGF signals and promote pulmonary pressure, miR-210 can inhibit mitochondrial function. These may promote the occurrence of HAPE [20]. However, the differential expression of miRNAs has not been reported in HAR.

In the present study, we aimed to explore HAR-related miRNAs and mRNAs via analyzing the GSE90500 and GSE75665 datasets from the Gene Expression Omnibus (GEO). Gene Ontology (GO) functional annotations, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, potential target gene analysis, and construction of potential miRNA–mRNA network were also conducted. Besides, the HAR cell model was established to validate the results of the bioinformatics analysis via quantitative RT-PCR (qRT-PCR). The present study may help elucidate the pathogenic mechanism and identify valuable diagnostic biomarkers for HAR.

The miRNA expression profile GSE90500 and the mRNA expression profile GSE75665 were obtained from the GEO database. The RNA-seq data of GSE90500 were based on the platform of GPL18058 (Exiqon miRCURY LNA microRNA array, 7th generation) and contained 13 patients with AMS and 9 non-AMS volunteers. The microarray data of GSE75665 were based on the platform GPL11154 (Illumina HiSeq 2000) and contained 5 patients with AMS and 5 non-AMS volunteers. Figure 1 shows the workflow for the study.

The GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) based on the R software ‘LIMMA’ package was used to identify differentially expressed miRNA (DE-miRNA) between AMS and non-AMS blood samples of GSE90500 with thresholds of |logFC| > 2 and P-value <0.01. A volcano plot was generated to visualize the significant miRNA expression changes using GraphPad Prism 8.0 (GraphPad, San Diego, CA, USA). For GSE75665, differentially expressed mRNAs (DE-mRNAs) related to AMS were obtained from the study of Liu et al. [21].

We used three prediction tools to predict the target mRNAs of DE-miRNA: TargetScan (http://www.targetscan.org/), miRDB (http://mirdb.org/), and DIANA-microT (http://www.microrna.gr/microT). The genes predicted by at least two programs were chosen as the targets of DE-miRNA. We selected the overlapping genes between targets and DE-mRNAs as target differentially expressed genes (DEGs). A heatmap was performed using the R package ‘Pheatmap’. Then, using Cytoscape software [22], the miRNA–mRNA regulatory network was constructed.

GO annotation and KEGG pathway enrichment analyses were performed by utilizing the Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.abcc.Ncifcrf.gov/) online software [23]. P<0.05 was set as the cut-off value. The GO analysis included three categories: biological process (BP), cellular component (CC), and molecular function (MF). Heatmaps of DEGs significantly enriched KEGG pathways was performed by the R package ‘Pheatmap’.

By putting the target DEGs in the Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) database [24], the interaction relationships among DEGs at the protein level were screened. Then, the protein–protein interaction (PPI) networks for up- and down-regulated target DEGs were constructed under the criterion of a combined score > 0.4. Cytoscape software was used to visualize the networks. The connectivity degrees were calculated through network statistical methods.

Human retinal microvascular endothelial cells (HRMECs) were purchased from Cell Systems Corporations (catalog no. ACBRI181; Kirkland, WA, USA) [25]. The HRMECs were cultured in M199 medium supplemented with 20% fetal bovine serum, 3 ng/ml FGF-basic, 10 units/ml heparin, and 1% streptomycin/penicillin at 37°C. To establish the hypoxia model, cell culture was performed in a hypoxia chamber filled with an anaerobic gas mixture of 94% N₂, 5% CO₂, and 1% O₂ for 12 h. Control groups were cultured in 95% air and 5% CO₂.

Total RNA was extracted from HRMECs by using TRIzol reagent. RNA was reverse-transcribed into complementary DNA using PrimeScript RT Master Mix or miRNA First-Strand Synthesis kits (TaKaRa, Kumamoto, Japan), following the manufacturer’s protocols. qRT-PCR was performed by TB Green Premix ExTaq (TaKaRa, Kumamoto, Japan). The oligonucleotide primers used for PCR amplification were purchased from BioSune Biotechnology (Shanghai, China) and are listed in Supplementary Table S1. ACTB and U6 served as internal control for gene and miRNA expression in analysis, respectively. The relative expression of the DE-miRNAs and the DE-mRNAs was calculated using the 2^−ΔΔC_t method.

The data in the present study are presented as the mean ± SD. The results were analyzed by the two-tailed Student’s t test, and P<0.05 was considered to show a statistical difference.

According to the analysis of Liu et al. for GSE75665 [21], a total of 807 genes were differentially expressed between the non-AMS and AMS groups, with 271 up-regulated and 535 down-regulated DE-mRNAs (Supplementary Table S2). During an analysis of the GSE90500 dataset, we obtained 16 DE-miRNAs (8 up-regulated and 8 down-regulated) by using the GEO2R online tool (P<0.01 and |logFC| > 2) (Table 1). These DE-miRNAs are indicated by a volcano plot (Figure 2). MiR-495 includes miR-495-3p and miR-495-5p, and we discussed them separately in the study, so there were 17 DE-miRNAs (9 up-regulated and 8 down-regulated).

The target mRNAs of the DE-miRNAs were predicted by the three programs. Besides, by comparing the targets with DE-mRNAs, we only selected the overlapping genes as target DEGs. As shown in Figure 3A and B, 122 out of 5652 predicted targets of up-regulated DE-miRNAs overlapped down-regulated DE-mRNAs, and 37 out of 3857 predicted targets of down-regulated DE-miRNAs overlapped up-regulated DE-mRNAs. After removing 1 duplicate, we got 157 target DEGs, and the expression of these DEGs was visualized in Figure 3C. Then, we constructed the miRNA–mRNA regulatory network (Figure 4), which consisted of 240 miRNA–mRNA pairs in total. Among them, there were 196 up-regulated miRNA–mRNA pairs and 44 down-regulated miRNA–mRNA pairs. Among the target DEGs, there were 29 transcription factors and 128 non-transcription factors. Table 2 listed the target DEGs of DE-miRNAs.

To further explore the biological function of the target DEGs, GO enrichment and KEGG pathway analysis were performed using DAVID. GO analysis results showed that acrosome assembly was the only significantly enriched BP for the up-regulated target DEGs. As shown in Figure 5A and Table 3, the down-regulated target DEGs were mainly enriched in the process of RNA expression, such as positive regulation of transcription from RNA polymerase II promoter, positive regulation of gene expression at BP level. On the CC level, the down-regulated DEGs were mainly enriched in the nucleus, cytoplasm, and plasma membrane. For the MF, the down-regulated DEGs were mainly enriched in sequence-specific DNA binding, transcriptional activator activity, and RNA polymerase II core promoter proximal region sequence-specific binding. In addition, KEGG pathway analyses indicated that down-regulated DEGs were significantly enriched in the Hippo signaling pathway (Figure 5B and Table 4), which is associated with the proliferation and migration of vascular endothelial cells [26,27] and can be deactivated by hypoxia [28]. In addition, as shown in Figure 6A–D, heatmaps illustrated the allocation of down-regulated DEGs significantly enriched four KEGG pathways.

Based on information from the STRING database, the PPI network was constructed and contained 59 nodes and 99 interactions (Figure 7). The details of the PPI network are described in Supplementary Table S3. In this network, with the criteria of filtering degree > 5, eight genes were defined as hub genes, including IL10, CDK1, FOS, HNRNPH1, IL7R, BCL11B, FCGR2A, and DDX3X (Table 5).

To validate the results of bioinformatics analysis, we chose the nine up-regulated DE-miRNAs and eight hub genes for qRT-PCR under normoxia and hypoxia in HRMECs. As shown in Figure 8A and B, miR-3177-3p, miR-369-3p, miR-603, miR-495-3p, miR-495-5p, miR-4791, and miR-424-5p were significantly up-regulated in hypoxia-cultured HRMECs compared with normoxia (P<0.05). FOS, IL10, and IL7R were strongly down-regulated under hypoxia conditions (P<0.05), which was also consistent with our bioinformatics results. However, the expression of DDX3X, BCL1B, HNRNPH1, and CDK1 was higher than that in normoxia-cultured HRMECs.

HAR is an ocular manifestation of AAI [1], and studies have found a statistically significant correlation between HAR and HACE [3]. Some researchers have proposed that AMS is a mild form of HACE [29,30]. Since the retina and optic nerve act as the directly visible part of the brain, belonging to the central nervous system due to its embryonic origin [31,32], we obtained the DE-miRNAs and DE-mRNAs of HAR by analyzing the AMS datasets and then explored the possible potential target genes of HAR.

MiRNAs can regulate many BPs by negatively regulating the expression of target genes and are thus important targets for studying the genetic susceptibility of diseases [19,33]. Through a series of bioinformatics analyses and experiments, we selected 16 DE-miRNAs and 157 target DEGs related to AMS from GSE90500 and GSE75665 and constructed an miRNA–mRNA network containing 240 relationship pairs. The top eight hub genes—IL7R, FOS, IL10, FCGR2A, DDX3X, CDK1, BCL11B and HNRNPH1—filtered from the PPI network may serve as candidate biomarkers of HAR. The expression of miR-3177-3p, miR-369-3p, miR-603, miR-495-3p, miR-495-5p, miR-4791, miR-424-5p, FOS, IL10, and IL7R were confirmed by qRT-PCR in our HAR cell model to have the same trend as that found in our bioinformatics results. These miRNAs and mRNAs may play an important role in the development of HAR.

FOS, consisting of c-Fos, FosB, Fra1, and Fra2 [34], was significantly down-regulated in our hypoxia-cultured HRMECs. In the lung tissues of hypoxic rats, the expression of c-Fos and FosB was significantly decreased [35], which is consistent with our results. c-Fos promotes the development of inflammatory diseases such as arthritis [36], but it also inhibits inflammation in myeloid and lymphoid cell lineages [37]. In studies on retinal development, it has been found that c-Fos induces inflammatory signals in photoreceptor cells, leading to an increase in VEGF, thus promoting the formation of new blood vessels [38]. Therefore, FOS may be associated with the leakage of the peripheral retinal vessels and retinal hemorrhages in HAR through the regulation of inflammatory constituents. The bioinformatics analysis showed that miR-603 and miR-4791 may down-regulate FOS, and miR-603 can promote the growth of glioma cells via the Wnt/β-catenin pathway [39]. The Wnt pathway mediates many complex biological and pathological processes of retinopathy [40], but there are no reports about miR-603 and this pathway in the development of HAR. Further research about the role of FOS and miR-603 in HAR is needed.

Another inflammation-related gene, IL10, could reduce tissue damage and promote immune tolerance [41]. It was down-regulated in our study, and maybe serve as a biomarker in HAR. Under hypoxic conditions, the decrease in IL10 indicates hypoxia reduces the ability to inhibit immunity and inflammation [42]. However, in the study of Kang et al., mice adapted to the low oxygen environment (10% O₂) had higher IL10 levels [43]. Therefore, the role of anti-inflammatory markers IL10 in HAR needs further investigation.

IL7R mediates the signaling of IL7, affecting the role of IL7 in regulating the development of T and B cells and the balance of the internal environment of T cells [44], and IL7R is involved in the pathogenesis of several autoimmune diseases, such as rheumatoid arthritis, multiple sclerosis, and systemic lupus erythematosus etc. [45] High altitude also could cause changes in different immune cells [46]. According to our bioinformatics analysis, we found that IL7R was regulated by miR-369-3p and miR-424-5p. It is noted that the induction of miR-369-3p expression inhibits the release of inflammatory cytokines [47]. Furthermore, up-regulated miR-424-5p in hypoxic vascular endothelial cells was shown to stabilize hypoxia-inducible factor-1α (HIF-1α) and fine-tune VEGF, promoting angiogenesis and preserving the integrity of the endothelial barrier [48]. Although the down-regulated IL7R and the up-regulated miR-369-3p and miR-424-5p may play an important role in the progression of HAR, the specific mechanism is not clear and remains to be clarified.

High-altitude retinal hemorrhages and vascular leakage may be caused by the retinal microcirculation disorder under hypoxia conditions at high altitude, but the specific mechanism is not clear. The expression of DDX3X was found to be correlated with overexpression of HIF-1α in breast cancer, indicating the oncogenic role of DDX3X [49]. CDK1 has been shown to promote tumor angiogenesis by stabilizing HIF-1α, and its disruption could inhibit retinal angiogenesis by inducing cell cycle arrest and apoptosis [50]. However, inhibition of miR-495 improves angiogenesis and blood flow recovery after ischemia [51]. Although the expression of DDX3X and CDK1 did not down-regulate in HRMECs as predicted, these genes had high degrees in the PPI network, indicating their potential roles in the development of HAR. It worth mentioning that very little was reported on miR-3177-3p, so the specific function of miR-3177-3p needs further research.

To our knowledge, our study is the first to explore the DE-miRNAs and DEGs of HAR and establish the miRNA–mRNA network through comprehensive bioinformatics analyses. However, there are still some limitations in the present study. First, the data of the GSE90500 and GSE75665 are collected from blood samples, which may be the reason why the expression of several hub genes in HRMECs was inconsistent with the bioinformatics analyses. Then, there are some differences between the analytical methods of GSE90500 and GSE75665, which reduces the reliability of this study to a certain extent. Third, the miRNA–mRNA relationships in HAR were based on target prediction and statistical evidence, which need experimental investigations. Thus, Further studies are required to verify the specific functions of the identified DE-miRNAs and hub genes in HAR.

In conclusion, the present study identified FOS, IL10, IL7R and seven miRNAs as candidate biomarkers of HAR, which may participate in the development of HAR through possible pathways. It has been proposed that HAR is related to individual sensitivity, promoting the effective prevention, and treatment of HAR in the future. Nevertheless, we carried out the systematic and comprehensive bioinformatics analyses to identify the new miRNA and genes of HAR and performed experimental validation in the HAR cell model, which is helpful for understanding the gene changes in HAR and provides novel candidate biomarkers for the diagnosis of HAR.

All data associated with the present paper are already included in the manuscript.

The authors declare that there are no competing interests associated with the manuscript.

This work was supported by the National Natural Science Foundation of China [grant numbers 81970811, 81770947]; the National Science and Technology Major Project of China [grant number 2017ZX09304010]; the National Key R&D Program of China [grant numbers 2016YFC0904800, 2019YFC0840607]; and the Tibet Natural Science Foundation of China [grant numbers XZ2017ZR-ZYZ09, XZ2018ZRG-95, XZ2018ZRG-100(Z)].

All authors contributed significantly to this work. Tong Su, Chufeng Cu, and Deji Draga were involved in all aspects of research, analyzed data and wrote the paper. Chuandi Zhou performed the experiments. Thashi Lhamo was involved in research design. Qinghua Qiu and Zhi Zheng contributed to research design, supervision of research, data analysis and revised the manuscript. All authors have read and approved the final manuscript.

AAI

acute high-altitude illness

AMS

acute mountain sickness

BP

biological process

CC

cellular component

DAVID

Database for Annotation, Visualization, and Integrated Discovery

DE-miRNA

differentially expressed miRNA

DE-mRNA

differentially expressed mRNA

DEG

differentially expressed gene

GEO

Gene Expression Omnibus

GO

Gene Ontology

HACE

high-altitude cerebral edema

HAPE

high-altitude pulmonary edema

HAR

high-altitude retinopathy

HIF-1α

hypoxia-inducible factor-1α

HRMEC

human retinal microvascular endothelial cell

KEGG

Kyoto Encyclopedia of Genes and Genomes

MF

molecular function

miRNA

microRNA

qRT-PCR

uantitative RT-PCR

STRING

Search Tool for the Retrieval of Interacting Genes