Hypoxia-inducible factor stabilisation-related lncRNAs in retinopathy of prematurity

Abstract

Long non-coding RNAs (lncRNAs) play an important role in the response to many diseases. The previous study reported the transcriptomes of mice that were cured of oxygen-induced retinopathy (OIR, retinopathy of prematurity (ROP) model) by hypoxia-inducible factor (HIF) stabilisation via HIF prolyl hydroxylase inhibition using the isoquinolone Roxadustat or the 2-oxoglutarateanalog dimethyloxalylglycine (DMOG). However, there is little understanding of how those genes are regulated. In the present study, 6918 known lncRNAs and 3654 novel lncRNAs were obtained, and a series of differentially expressed lncRNAs (DELncRNAs) were also identified. By cis- and trans-regulation analyses, the target genes of DELncRNAs were predicted. Functional analysis demonstrated that multiple genes were involved in the MAPK signalling pathway, adipocytokine signalling pathway was regulated by the DELncRNAs. By HIF-pathway analysis, two lncRNAs Gm12758 and Gm15283 were found that can regulate the HIF-pathway by targeting the Vegfa, Pgk1, Pfkl, Eno1, Eno1b and Aldoa genes. In conclusion, the present study provided a series of lncRNAs for further understanding and protecting the extremely premature infant from oxygen toxicity.

Impact Statement

• What is already known on this subject? Roxadustat can prevent oxygen-induced retinopathy (OIR) by two pathways: direct retinal hypoxia-inducible factor (HIF) stabilisation and induction of aerobic glycolysis or indirect hepatic HIF-1 stabilisation and increased serum angiokines. However, underlying the long non-coding RNAs (lncRNAs) that may regulate the HIF stabilisation-related genes have not been investigated thoroughly.

• What do the results of this study add? Six thousand nine hundred and eighteen known lncRNAs and 3654 novel lncRNAs were identified. GO and KEGG enrichment analysis showed that the MAPK signalling pathway and adipocytokine signalling pathway were regulated by the differentially expressed lncRNAs (DELncRNAs). Two lncRNAs Gm12758 and Gm15283 were found that may regulate the HIF-pathway by targeting the Vegfa, Pgk1, Pfkl, Eno1, Eno1b and Aldoa genes.

• What are the implications of these findings for clinical practice and/or further research? It provides a further rationale for protecting severe premature infants from oxygen poisoning.

Keywords:

• ncRNAs
• hypoxia-inducible factor
• Roxadustat
• DMOG
• target genes

Introduction

Recent evidence has demonstrated that more than 90% of the mammalian genome is transcribed into non-coding RNAs (NcRNAs), which cannot be considered junk or transcriptional noise due to their biological significance (Jathar et al. 2017). Most NcRNAs, such as rRNA, tRNA and snoRNA, are smaller than 200 nucleotides (nt). Those longer than 200 nt are called long non-coding RNAs (lncRNAs), which were predicted largely due to the efforts of the FANTOM, GENCODE and the ENCODE consortia (Khorkova et al. 2015). These RNAs play important roles in gene expression control in developmental processes, such as dosage compensation, genomic imprinting, cell differentiation and organogenesis (Wu et al. 2020, Cai et al. 2021). Many reports have shown that lncRNAs located in the nucleus regulate gene expression at the mRNA level in the cytoplasm (Beermann et al. 2016). Hung et al. (2011) found that 56 cell-cycle genes harbour lncRNAs and regulate gene expression in human cancers. The expression levels of these genes are regulated by specific oncogenic stimuli, stem cell differentiation or DNA damage (Fatica and Bozzoni 2014). Huarte et al. (2010) found that the inhibition of lncRNA-p21 affects the expression of hundreds of targets genes normally repressed by p53. Tissue differentiation inducing non‑protein coding RNA is a key lncRNA that binds to differentiation mRNAs, including FLG, LOR, ALOXE3, ALOX12B and ABCA12, to ensure their expression (Huarte et al. 2010). These examples show that lncRNAs play crucial roles in biological processes (BPs) by targeting mRNAs or other RNAs. Furthermore, lncRNAs have been reported to be associated with many diseases, including cancer, obesity, cardiovascular disease and visual impairment. An improved understanding of the role of lncRNAs in vascular disease may provide new pathophysiological insights, given recent findings on the expression, mechanism and function of lncRNAs implicated in a range of vascular disease states, ranging from those occurring in mice to human subjects (Simion et al. 2019). Cissé et al. (2018) expounded the characteristics of lncRNAs and their regulation in glaucoma and Biswas et al. (2019) highlighted some of the key roles of lncRNAs in diabetic retinopathy and inflammation.

In 1942, retrolental fibroplasia (RLF) was first reported in retinopathy of prematurity (ROP), which occurs in premature infants born under 36-week gestation with a low birth weight and long-term oxygen supplementation. Uncontrolled supplemental administration of oxygen was found to be the cause of RLF (Sola and Zuluaga 2013). Retinopathy of prematurity was called RLF in the 1940s (Terry 1942). Retinovascular growth attenuation and vascular obliteration created by hyperoxia cause ROP, resulting in approximately 100,000 new cases of childhood blindness each year (Hartnett and Penn 2012). Clinical trials have demonstrated that oxygen saturation reduces ROP, but increases mortality. A recent study found significant international practice variation and uncertainty pertaining to target oxygen saturation in premature infants (Terry 1942). Roxadustat (RXD) and dimethyloxalylglycine (DMOG) used in hypoxia-inducible factor (HIF) stabilisation are reported to have the potential to cure infants with ROP through oxygen tension response, and the cells were regulated by HIFs (Wang et al. 1995).

Hoppe et al. (2016) used organ systems pharmacology to define the transcriptomes of mice that were cured of oxygen-induced retinopathy (OIR) (OIR, ROP model) by HIF stabilisation via HIF prolyl hydroxylase inhibition using the isoquinolone RXD or its 2-oxoglutarate analogue DMOG. In the present study, the expression patterns of lncRNAs in each liver and retina sample, and their effect on the target genes were analysed. Six thousand nine hundred and eighteen known lncRNAs and 3654 novel lncRNAs were predicted, and many genes, especially involved in the MAPK signalling pathway, adipocytokine signalling pathway, and the ABC transport pathway that are regulated by lncRNA were observed. We aimed to elucidate the regulatory action of lncRNAs in HIF stabilisation mechanisms in vivo. The present study showed that the differentially expressed lncRNAs (DELncRNAs) may regulate the genes involved in the process of reducing oxygen toxicity. For instance, two lncRNAs Gm12758 and Gm15283 can regulate the HIF-pathway by targeting the Vegfa, Pgk1, Pfkl, Eno1, Eno1b and Aldoa genes, and these two lncRNAs are important for the regulation mechanism of combating oxygen toxicity in RXD treatment.

Methods

Data processing

All the RNA-Seq data of OIR mice treated with RXD and DMOG were downloaded from the Gene Expression Omnibus (GEO) database (accession no. GSE74170; NCBI tracking system no. 17567121) (Hoppe et al. 2016). Sequencing adaptors and low-quality sequences were trimmed using Trimmomatic (Bolger et al. 2014). The clean reads were aligned to the mouse reference genome GRCm38 using HISAT2 (Kim et al. 2015, 2019).

Identification of lncRNAs

To obtain qualified lncRNA candidates, we performed a preliminary filtering of transcripts and predicted the coding potential of the screened transcripts (Supplementary Figure 3). Four criteria were employed in the filtering step: (1) transcript length and exon number should be above 200 bp and 2, respectively; (2) transcripts should be covered by at least five reads in all samples; (3) transcripts aligned to known mRNAs and other NcRNAs (such as rRNA, tRNA, snoRNA and snRNA) of this species should be discarded; (4) the potential lncRNA, intronic lncRNA and anti-sense lncRNA were predicted according to class_code information (‘u’, ‘i’, ‘x’). Due to the lack of coding ability, the candidate lncRNA transcripts were further assessed for their coding potential using four methods, including coding-non-coding index (CNCI) (Sun et al. 2013), coding potential calculator (CPC) (Kong et al. 2007) and coding potential assessment tool (CPAT) (Wang et al. 2013), which were used to calculate the encoding and the non-coding capabilities of the transcripts. In addition, the protein domains in the HMM library were searched for each candidate sequence using the pfamscan (https://www.ebi.ac.uk/Tools/pfa/pfamscan/) tool to screen out sequences with known protein domains.

Differential expression lncRNAs analysis

The DESeq2 package in R (Love et al. 2014) was used for the lncRNA differential expression analysis. Differentially expressed lncRNAs (q≤ 0.05 and |log2 ratio| ≥1) were denoted as DElncRNAs. For DElncRNAs, cis- and trans-regulation analyses were performed. Cis-regulation refers to the regulatory action of lncRNAs on the adjacent protein-encoding genes in the genomic locus. The target genes for cis-regulation were predicted by screening protein-encoding genes adjacent to the lncRNA (upstream and downstream 50 kb). The trans-regulated genes were recognised based on the correlation coefficient of the lncRNA and the mRNA expression values (Pearson’s correlation coefficient ≥0.9, Pearson’s correlation coefficient ≤−0.9, respectively).

GO and KEGG enrichment analysis

The cis-acting regulatory genes and the trans-acting regulatory genes of DElncRNAs obtained from the retina and the liver in RXD and DMOG treatment were used to perform GO (Blake et al. 2015) and KEGG (Kanehisa 2000) enrichment analysis, for which the hypergeometric test was applied. An adjusted p value of .05 was used as the threshold for significant enrichment.

DElncRNAs and differentially expressed genes (DEGs) interaction network analysis

All DEGs were downloaded from the Supplementary Information of Hoppe et al. (2016). Common and unique DElncRNAs expressed in different organs, under RXD and DMOG treatment, were identified for further analysis. An adjusted p value cut-off (q≤ 0.05) was used to determine the significance of the pathways identified. DElncRNAs from commonly appearing DEGs, for which targeted regulation was identified upon RXD and DMOG treatment (downloaded from Hoppe, Yoon (Hoppe et al. 2016)), were selected for interactive network analysis. Interactive networks involving more than two DElncRNAs were explored using the Cytoscape software (Shannon et al. 2003).

HIF-pathway analysis

Genes pertaining to the HIF-1 signalling pathway in mouse (Mus musculus) (mmu04066) were downloaded from the KEGG PATHWAY database. A total of 112 genes were identified to be involved in the HIF-1 signalling pathway. Common and unique DElncRNAs (q≤ 0.05), for which targeted regulation was identified under RXD and DMOG treatment, were selected to draw an interactive network diagram.

Results

Data pre-processing and lncRNA identification

The RNA-seq data of 36 samples from OIR mice subjected to RXD and DMOG treatment were obtained from a previous study (Hoppe et al. 2016). Every data group, including Liver_PBS, Liver_DMOG, Liver_RXD, Retina_PBS, Retina_DMOG and Retina_RXD, contained six biological repeats. The groups Liver_PBS and Retina_PBS were used as control groups. After removing sequencing adaptors and low-quality sequences, clean RNA-seq data from each sample were aligned to the mouse reference genome using HiSAT2 (Kim et al. 2019). The average total reads, mapped reads and the mapping rate were 92.97 million, 86.97 million and 93.54%, respectively (Table 1), and there were no significant differences among the six groups.

Table 1. RNA-seq data summary.

#Library	Total reads	Mapped reads	Mapping rate
Liver_PBS	98,526,919	91,335,829	92.74%
Liver_DMOG	95,106,800	89,609,117	94.25%
Liver_RXD	87,799,495	82,661,434	94.10%
Retina_PBS	83,426,826	76,767,316	92.08%
Retina_DMOG	99,742,039	93,549,367	93.90%
Retina_RXD	93,244,413	87,898,398	94.14%
Average	92,974,415	86,970,244	93.54%

A total of 6918 known lncRNAs were identified in the 36 samples. A total of 4593, 4593, 4593 and 6052 novel lncRNAs were obtained using the CNCI, CPC and CPAT software, and the results of all four methods (Supplementary Figure 1C, 1D and 1F). As shown in supplementary Figures 1C and 1D, most lncRNAs had two exons and were shorter than 1000 bp. Compared to the lncRNAs, the mRNAs presented a larger exon number and a longer gene length (Supplementary Figures 1E and 1G). The fragments per kilobase of transcript (FPKM) of the transcripts corresponding to the mRNAs and lncRNAs in the 36 samples was calculated in order to explore the expression level of the identified lncRNAs (Supplementary Figure 4A, 4B and 4C); the results showed a lower expression level of lncRNA (Wilcoxon’s rank sum test; p < .001) than that of the protein-coding genes (Supplementary Figure 1B). An entropy-based metric relying on the Jensen-Shannon (JS) distance was employed to calculate the expression-specific score of each transcript. For each transcript, the maximum JS score of each RNA was considered to be either tissue-specific or treatment-specific. A large proportion of the lncRNAs exhibited a higher maximal JS score (Supplementary Figure 1A), indicating that the identified lncRNAs had a stronger tissue-specificity.

Differential expression lncRNA results

Compared to the control group (Liver_PBS group or Retina_PBS group), significantly (q≤ 0.05) DELncRNAs upon DMOG and RXD treatment identified in a data group were defined as DElncRNAs (Figure 1(A–D)). In the liver data groups, 17 known DElncRNAs were identified upon DMOG and RXD treatment, while 62 known DElncRNAs were identified only in the DMOG treatment group, and 25 known DElncRNAs were identified only in the RXD treatment group (Figure 1(E)). In the retina data groups, two downregulated known DElncRNAs were identified upon DMOG and RXD treatment, while 19 upregulated known DElncRNAs were identified only in the DMOG treatment group and 13 known DElncRNAs were identified only in the RXD treatment group (Table 2; Figure 1(G)). In our study, the number of DELncRNAs in retinal was less than that in the liver, while DMOG influenced more lncRNAs than RXD. The details of the DElncRNAs, including lncRNA FPKM, fold change and p values, are presented in Supplementary Table I (S1).

Figure 1. (A–D) Volcano plots showing the differentially expressed lncRNAs. (F) Venn diagram of the overlap number of DElncRNAs. (E, G) Histogram of differentially expressed known lncRNAs. DElncRNAs were defined as the lncRNAs with significant (q≤ 0.05) differential expression upon Roxadustat (RXD) and dimethyloxalylglycine (DMOG) treatment, when compared with the control groups (Liver_PBS group or Retina_PBS group).

Table 2. The numbers of known differentially expressed lncRNAs.

Organ	DElncRNAs	Total	Upregulated	Downregulated
Liver	Common	17	4	13
	Unique to DMOG	62	43	19
	Unique to RXD	25	13	12
Retina	Common	2	0	2
	Unique to DMOG	19	19	0
	Unique to RXD	13	6	7

DElncRNAs: differentially expressed lncRNAs.

Function enrichment analysis of the target genes of DElncRNAs

To explore the function of the DElncRNAs, GO and KEGG pathway analyses were performed on the target genes of the DELncRNAs. Among the liver groups, 439 BPs terms, 89 cellular components (CCs) terms, 72 molecular functions (MFs) terms and 76 KEGG pathways were found to be significantly enriched upon DMOG treatment. In the RXD treatment group, 334 BPs terms, 55 CC terms, 62 MF terms and 37 KEGG pathways were found to be significantly enriched (p ≤ .05). As shown in Figure 2(A,B) and Supplementary Table II(S2), several signal transport pathways, such as the MAPK signalling pathway, adipocytokine signalling pathway and the ABC transport pathway, were found to be enriched.

Figure 2. Function enrichment results of the target genes of DElncRNAs. (A) Top 10 enrichment KEGG pathways of target genes of DElncRNAs between liver DMOG treatment samples and liver PBS treatment samples. (B) Top 10 enrichment KEGG pathways of target genes of DElncRNAs between liver RXD treatment samples and liver PBS treatment samples. (C) Heatmap of enrichment q-value value of GO terms (biological processes). (D) Heatmap of enrichment FDR value of GO terms (cellular component GO terms).

Among the retina groups, the number of the enriched functions and pathways were found to be drastically reduced not only in the DMOG group but also in the RXD group. In the GO analysis, five CCs and one MF were found to be enriched significantly upon DMOG treatment, while 27 CCs were found to be enriched significantly upon RXD treatment. In the DMOG and the RXD groups, three significant KEGG pathways and one significant KEGG pathway were respectively found (Table 3). The details of the GO enrichment and the KEGG pathway enrichment results are shown in Supplementary Table II (S2).

Table 3. The numbers of GO terms and KEGG pathways of the target genes of differentially expressed lncRNAs.

Group	KEGG	GO
		BP	CC	MF
Liver_DMOG	76	439	89	72
Liver_RXD	37	334	53	62
Retinal_DMOG	3	0	5	1
Retinal_RXD	1	0	27	0

BP: biological processes; CC: cellular components; MF: molecular functions.

Among the liver groups, many GO terms related to stress response, such as regulation of response to stimulus and regulation of response to stress, were found to be enriched. Moreover, many GO terms related to signal transport, such as regulation of signalling, intracellular signal transduction, and regulation of signal transduction, were also found to be enriched (Figure 2(C)). However, in the retinal groups, these GO terms were not found to be enriched.

DElncRNAs and DEGs interaction networks

All DEGs data were downloaded from the Supplementary Information of the study by Hoppe et al. (2016). Correlation analysis of DElncRNAs and DEGs was performed for the six groups. Among the liver groups, the target genes of six DElncRNAs were found in both the DMOG and the RXD treatment groups, while target genes of 20 DElncRNAs and six DElncRNAs were found in the DMOG treatment and the RXD treatment groups, respectively. Among the retina groups, the target genes of five DElncRNAs and three DElncRNAs were found to be DEGs in the DMOG treatment and the RXD treatment groups, respectively. The details of these results are provided in Supplementary Table III (S3). As shown in supplementary Figure 2A and 2B, it can be found that many targeted genes of DElncRNAs were also differentially expressed, while in different tissues, they had different expression patterns. In different tissues, there were more DEGs and DElncRNAs in the liver than in the retina, regardless of whether RXD or DMOG was used for treatment. In both the liver and the retina, RXD treatment resulted in more DEGs and DElncRNAs. Most DEGs and DElncRNAs were not the same genes in different tissues.

Effects of lncRNA expression level on HIF-pathway

Data pertaining to the genes involved in the HIF-1 signalling pathway of Mus musculus (mmu04066) were downloaded from the KEGG PATHWAY database. Similar to the complete DElncRNAs and the DEG networks, fewer differential genes were found in the retina than in the liver, which also matched the number of DEGs in the study by Hoppe et al. (2016). In the liver, both DMOG and RXD stabilisation action were found to play crucial roles with a functional difference, while only RXD stabilisation action had significance in the retina (Figure 3). This result matched the conclusion obtained from the mRNA analysis in the study by Hoppe et al. (2016), which reported that RXD stabilisation action affects the HIF-1α or HIF-2α pathway in the retina and the liver, while DMOG plays a role only in the liver. Although in both liver and retina samples, only Gm12758 and Gm15283 regulated the genes involved in the HIF-pathway, more genes such as Pfkp, Eno1, Eno1b, Vhl, Slc2a1, Pfkl, Pgk1, Vegfa and Aldoa were regulated in different organ samples. Among these genes, Vhl, Eno1b, Eno1, Aldoa, Pgk1, Pfkl and Vegfa both appeared in RXD treatment. The details of these results are provided in Supplementary Table III (S3).

Figure 3. The interaction network of DElncRNAs and the target genes in the HIF-1 pathway. (A) In the liver; (B) in the retina. Genes involved in the HIF-1 signalling pathway of Mus musculus (mmu04066) were downloaded from the KEGG PATHWAY Database. The cis-acting regulatory genes or the trans-acting regulatory genes of DElncRNAs were defined as the target genes.

Discussion

The RNA-seq dataset was downloaded from the study by Hoppe et al. (2016), which analysed the comparative systems pharmacology of HIF stabilisation in relation to the prevention of ROP. In this study, the transcriptomes of mice that were cured of OIR upon RXD and DMOD treatment were analysed separately. It was found that the stabilisation action of RXD rescued the mice from retinal oxygen toxicity, whereas DMOG did not exhibit this effect. Our study design included a systematic analysis of related lncRNAs, which is more comprehensive and reveals the lncRNA regulation mechanism of the pathways found by Hoppe et al. (2016).

From our results, we found some DElncRNAs in the liver and the retina. Some of them were identified in both the DMOG and the RXD treatment groups, while others were identified in only one of the treatments (Table 2). These results showed that lncRNAs are involved in the regulation processes occurring in both RXD and DMOG treatment. DElncRNAs were more abundant in the liver than in the retina (104/34), while more DElncRNAs were identified in the DMOG treatment than in the RXD treatment, both in the liver (62/25) and in the retina (19/13). In the KEGG/GO enrichment results, a similar trend was observed. This result matched the conclusion of Hoppe et al. (2016) that both RXD and DMOG exhibit stabilisation action in the liver, while only RXD stabilisation action plays a crucial role in reducing oxygen toxicity in the retina.

Many genes such as Pfkp, Eno1, Eno1b, Vhl, Slc2a1, Pfkl, Pgk1, Vegfa and Aldoa are regulated by DElncRNAs in different organs. Among these genes, Vhl, Eno1b, Eno1, Aldoa, Pgk1, Pfkl and Vegfa were associated with the HIF-pathway in RXD treatment. These genes seem likely to have significance in ROP treatment, necessitating further functional research. The Vhl gene is regulated in von Hippel-Lindau (VHL) disease identified by benign and malignant tumours, including retinal and central nervous system hemangioblastomas (Hos et al. 2019). Eno1 was found to be upregulated in the retinal pigment epithelial cells in response to hypoxia (Zheng et al. 2016). Chen et al. (1987) found that Pgk1 mRNA in retinal cells decreased with the prolongation of post-mortem interval. Vegfa plays a key role in the neovascularisation and the neuroprotection of the retina, and differential gene expressions of Vegfa, Aldoa and Pgk1 were identified in the retinal choroid by RT-PCR (Solberg et al. 2010). Almasry et al. (2018) found a significant upregulation of Vegfa and prom1 immunoreactivity in diabetic retinas compared to controls.

In summary, our data matched the mRNA analysis results from Hoppe et al. (2016), which showed that DMOG stabilisation action involved pathways for retinovascular protection against OIR, while RXD played a crucial role both in the retina and the liver. This result demonstrated that the use of low-dose, intermittent HIF prolyl-hydroxylase inhibitor can be used to treat oxygen toxicity. Through co-expression analysis, we found that some DElncRNAs regulated key genes in the HIF-1 pathway, including Vegfa, Pgk1, Eno1 and Vhl, are involved in retinal activity and retinal diseases. In conclusion, DElncRNAs are regulatory factors involved in the related pathways in the process of reducing oxygen toxicity. For instance, Gm12758 and Gm15283 regulate the HIF-pathway in the retina and the liver, regulating the expression of Vegfa, Pgk1, Pfkl, Eno1, Eno1b and Aldoa, and these two DElncRNAs are crucial for the complete mechanism of combating oxygen toxicity in RXD treatment.

Ethical approval

This study does not involve moral and ethical issues.

Author contributions

DQC conceived and designed the experiments. MKD and ZHC contributed to paper writing. YMC and ZC contributed to the data analysis.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data were downloaded from the Gene Expression Omnibus (GEO) database, www.ncbi.nlm.nih.gov/geo (accession no. GSE74170; NCBI tracking system no. 17567121).

Additional information

Funding

This work was supported by Zhejiang Natural Science Youth Science Foundation [LQ19H040008], the Key Research and Development Plan of Zhejiang Provincial Department of Science and Technology [2018C03010] and the Natural Science Foundation of Zhejiang Province [LQ20H040003].