RNA sequence analysis identified bone morphogenetic protein-2 (BMP2) as a biomarker underlying form deprivation myopia.

Purpose: Form deprivation myopia (FDM) is an urgent public issue characterized by pathological changes, but the underlying mechanism remained unclear. The aim was to investigate bone morphogenetic proteins (BMPs) utilizing the pathogenesis of FDM. Material and methods: Gene expression omnibus (GEO) database was used to analyze one mRNA profile (GSE89325) of FDM. Sixteen retina samples (8 FDM and 8 controls) were randomly divided into seven groups for differential gene expression analysis in R. software. The gene pathway and protein-protein interaction (PPI) analysis were performed by the DAVID and STRING databases. Cytoscape was used to draw the PPI network. The gene ontology (GO) enrichment and Kyoto encyclopedia of genes and Genomes (KEGG) analysis were determined to achieve gene annotation and visualization.
Results: A total of 18420 differentially expressed genes (DEGs) were identified associated with FDM. The only non-significant gene (BEND6) was separately analyzed between two groups. Thirteen hub genes were discovered, ACVR1, ACVR2A, ACVR2B, RGMB, BMPR2, BMPR1A, BMP2, BMPR1B, CHRD, PTH, PTH1R, PTHLH, and WNT9A. The expression alteration in FDM were mainly enriched in cytokine-cytokine, and neuroactive ligand receptor interaction pathways. BMP2 was the key gene in myopia progression. Conclusions: Of clinical perspective, our findings reveal that expression of BMP2 as an underlying mechanism of FDM, providing an insight for therapeutic interventions.

Introduction

The growth factor role of bone morphogenetic proteins (BMPs) has become an important area in numerous cellular functions, including chondrogenesis, cell proliferation, and extracellular matrix synthesis [1,2]. BMPs are crucially known for their salient sector in embryogenesis and development. Nevertheless, BMPs have been generally reported to promote early eye morphogenesis, collagen degeneration, and induce chondrogenic differentiation in human scleral fibroblasts [3]. Induction in chondrogenic differentiation has been regarded to be associated with some diseases, such as rheumatic arthritis and high myopia [[4], [5], [6]]. Apart from their efficiency in bone formation and repair, BMPs have other biological effects. For example, BMP2 is essential for retinal development. Both BMP2 and BMP4 are implicated in the formation of cardiac septa, and their deficiency might cause heart defects in mice. Moreover, BMP7 contributes to the deterioration of the kidney and the heart [[7], [8], [9]].

Myopia is a major cause of visual impairment among children. It has been related to ocular complications, such as macular degeneration, cataract, retinal detachment, and glaucoma. These comorbidities could increase the risk of irreversible visual loss [10,11]. The prevalence of myopia and high myopia is constantly rising globally and estimated to be 50% by the year 2050 [12,13]. Furthermore, early-onset myopic children have higher progression rates of high myopia [14,15]. However, the effect of BMP2 on the pathogenesis of myopia remains ambiguous. With the proliferation of low dose atropine [16], treatments for myopia have been largely evolved. Unfortunately, because low dose atropine induces various significant adverse events [17] and does not address the progression of high myopia, that is, there is an urgent need for new therapeutic methods. Therefore, a better understanding of the pathological mechanism for myopia progression will benefit the refreshment of clinical treatments.

Small mRNA suggests a role for post-transcriptional regulation that a light-dependent retina-to sclera signaling cascade and delineates possible pathobiological molecular drivers causing refractive errors [18]. Identification of differentially expressed mRNA and genes could be recognized through computational methods [19].

In this study, we identified the potential of BMPs as well as the difference in myopic retina between early and late timepoint recovery from FDM. Further analysis of the differentially expressed genes (DEGs) by functional annotation was performed to acknowledge the molecular basis and relevant pathways. Eventually, protein-protein interaction (PPI) network was analyzed to identify the key genes. We sought to explore the role of BMPs in FDM modulation to find a novel therapeutic target for myopia.

Material and methods

Data source and processing

We downloaded RNA sequencing dataset GSE89325 [20] as well as the corresponding characteristics of myopic retina and controls from the Gene Expression Omnibus (GEO) website (http://www.ncbi.nlm.nih.gov/geo/) in October 2021. The profile of GSE89325 was based on the platform of GPL21476 Illumina HiSeq 1500 (Gallus gallus). The RNA sequencing data were sorted by controls and FDM male chicks (Leghorn/New Hampshire), which were left untreated, and treated at different timepoints (0 h, 6 h, 24 h).

The expression levels were initially standardized and filtered using Fast-QC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/), multi-QC with matching reads beyond 80% [21], and Trim Galore (0.4.0), respectively. High quality sequence reads were aligned to chick genome (GalGal6) using Hisat2 [22] to obtain the counts for each sample.

Identification of DEGs

Our differential expression analysis based on RNA sequencing data was interpreted using the DEseq2 package (http://www.bioconductor.org/packages/release/bioc/html/DEseq2.html) in the R statistical environment (version 3.5.1) [23] to identify DEGs. The timepoint referred to occlusion (control vs myopia, 0 h vs control) and the 24 h recovery period (6 h vs 0 h, 24 h vs 6 h, and 24 h vs 0 h). Expression comparisons were conducted by moderated t-test method. The thresholds were |log 2fold change|>2, and false discovery rate (FDR) < 5% as the criteria for DEGs selection.

Gene annotation analysis

For those RNA target genes, Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) (http://www.genome.jp/kegg/pathway.html) analyses were performed on gene annotation. Pathway analysis was conducted through the Database for visualization, and integrated discovery (DAVID) (https://david.ncifcrf.gov/) [24,25]. Chick Ensemble Gene IDs were mapped to human gene symbol using BioMart (Ensemble release 101). While mapping detected no existing gene, they were connected to human gene orthologues (Homo sapiens GRCh38.p7).

For DEGs in myopic retina, PPI analysis was performed through the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.string-db.org) and produced by Cytoscape software with a confidence score of >1.0 [26].

Results

Data quality assessment

We searched 16 samples (8 myopia and 8 controls) from GEO database and then transferred into SRR documents using R software. Finally, a total of 18420 DEGs were identified between myopia and control groups. Through fast-QC and multi-QC assessment, we found 80% non-reproductive index sequence counts. The number of unique reads mostly centered on 6 reads, which demonstrated that RNA sequencing data were confidential. The plot was displayed in Fig. 1. We compared seven groups between controls and chicks recovery from FDM (0hrs, 6hrs, 24hrs).

Fig. 1

Plot of quality control for 15 SRR documents using multi-QC. Duplicate reads were identified using trim-gallore. QC, quality control.

Total analyses (FDM versus controls)

We compared control and FDM groups as shown in Fig. 2. The volcano plot indicated that only 1 gene (BEND 6) was detected to be differentially expressed (Fig. 2a). There was no statistically significant difference (p > 0.05). We found that 16 mRNAs expressions were obviously different from others. Among those mRNAs that 8 were upregulated with size Factor of 1.2 as shown in the heatmap (Fig. 2b).

Fig. 2

Differentially expressed genes and mRNAs in myopic retina. The volcano plot (a) showed that a total of 1 non-significant gene was detected to be differentially expressed. The heat map (b) detected that 16 mRNAs expressed differentially between myopic retina and controls, and 15 of them were upregulated.

Comparisons between early and late FDM versus controls

We compared the profiles of myopic and normal chicks, respectively. A total of 131 genes were detected as DEGs among seven groups. The volcano plots were displayed in Fig. 3. The volcano plot (a) showed that a total of 3 genes were detected to be differentially expressed between control and FDM 0 h; 2 were up-regulated and 1 was down-regulated. The volcano plot (b) displayed that 13 genes (up-regulated 8, down-regulated 5) differentially expressed between control and FDM 6 h. The volcano plot (c) indicated that 23 genes (up-regulated 7, down-regulated 15) differentially expressed between control and FDM 24 h. In early stage FDM, there were 23 DEGs about FDM 0 h vs FDM 6 h (d). Conversely, in late stage FDM, 26 and 42 DEGs were identified about FDM 0hrs vs FDM 24 h (e), and FDM 6hrs vs FDM 24 h, respectively (f). Altogether, only WNT9A was down-regulated in response to both FDM 0hrs vs FDM 6 h and FDM 0hrs vs FDM 24 h.

Fig. 3

Differentially expressed genes in early and long term FDM retina. The volcano plot (a) showed that a total of 3 genes were detected to be differentially expressed between control and FDM 0 h; 2 were upregulated and 1 downregulated. The volcano plot (b) displayed that 13 genes expressed differentially between control and FDM 6 h, and 8 of them were upregulated. The volcano plot (c) displayed that 23 genes expressed differentially between control and FDM 24 h, and 7 of them were upregulated. In early FDM retina, 23 genes were differentially expressed in response to (FDM 0hrs vs FDM 6 h) (d), 26 genes in response to (FDM 0hrs vs FDM 24 h) (e), and 42 genes in response to (FDM 6hrs vs FDM 24 h) (f). Altogether, only WNT9A was downregulated in response to both FDM 0hrs vs FDM 6 h and FDM 0hrs vs FDM 24 h. FDM 0hrs: 0 h recovery from FD. FDM 6hrs: 6 h recovery from FD; FDM 24hrs: 24 h recovery from FD.

GO and KEGG analysis

We conducted GO analysis of the DEGs in seven groups. Meanwhile, several comparisons were obtained among 0 h versus 6 h, 0 h versus 24 h, and 0 h versus control. A total of 44 biological processes (BP), 12 molecular functioning (MF), and 6 cellular components (CC) were found through the STRING database, and were mainly metabolism associated receptors, including signaling receptor activity, cell surface receptor signaling pathway, activin-activated receptor activity, and activin receptor complex. Moreover, the detected 3 KEGG pathways contributed to cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, and TGF-beta signaling pathway (Fig. 4). The top 6 GO and KEGG pathways were shown in Table 1. Cell surface receptor signaling pathway was one of the top 6 BPs, and the remaining was BMP signaling pathway, positive regulation of developmental process, regulation of cell differentiation, positive regulation of angiogenesis, animal organ morphogenesis, axis specification, regulation of response to stimulus and regulation of animal organ morphogenesis.

Fig. 4

Column plots showing all significant pathways from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway comparisons during the induction of myopia and recovery from FDM. The significantly enriched pathways are visualized for (a) 0 h versus 6 h, (b) 0 h versus 24 h, and (c) 0 h versus Control. The size of each column is proportional to the number of core genes within the pathway.

Table 1

The top 6 Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways of the differentially expressed mRNA targets between control and FDM 6 h.

Category	GO term	Description	Count	Genes
BP	GO:0051239	regulation of multicellular organismal process	8	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:0030154	cell differentiation	8	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:0007165	signal transduction	8	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:0048731	system development	8	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:2000026	regulation of multicellular organismal development	7	WNT9A,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:0007166	cell surface receptor signaling pathway	7	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,ACVR2A,RGMA
MF	GO:0005488	binding	7	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A
	GO:0005515	protein binding	6	WNT9A,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A
	GO:0038023	signaling receptor activity	5	ACVR2B,BMPR1B,ACVR1,ACVR2A,RGMA
	GO:0004675	transmembrane receptor protein serine/threonine kinase activity	4	ACVR2B,BMPR1B,ACVR1,ACVR2A
	GO:0017002	activin-activated receptor activity	3	ACVR2B,ACVR1,ACVR2A
	GO:0048185	activin binding	3	ACVR2B,ACVR1,ACVR2A
CC	GO:0110165	cellular anatomical entity	8	WNT9A,ACVR2B,BMP2,BMPR1B,ACVR1,PTHLH, ACVR2A,RGMA
	GO:0005886	plasma membrane	5	ACVR2B,BMPR1B,ACVR1,ACVR2A,RGMA
	GO:0031224	intrinsic component of membrane	5	ACVR2B,BMPR1B,ACVR1,ACVR2A,RGMA
	GO:0043235	receptor complex	4	ACVR2B,BMPR1B,ACVR1,ACVR2A
	GO:0005887	integral component of plasma membrane	4	ACVR2B,BMPR1B,ACVR1,ACVR2A
	GO:0048179	activin receptor complex	3	ACVR2B,ACVR1,ACVR2A
KEGG pathways	gga04350	TGF-beta signaling pathway	10	BMPR1A,ACVR2B,BMPR2,BMP2,BMPR1B,ACVR1,RGMB, ACVR2A,CHRD, RGMA
	gga04060	Cytokine-cytokine receptor interaction	7	BMPR1A,ACVR2B,BMPR2,BMP2,BMPR1B,ACVR1,ACVR2A
	gga04080	Neuroactive ligand-receptor interaction	4	PTGFR/F2RL1/PRLL/EDNRB

PPI analysis was performed on 13 DEGs as the following activin A receptor type 1 (ACVR1), activin A receptor type 2A (ACVR2A), activin A receptor type 2B (ACVR2B), repulsive guidance molecule BMP co-receptor b (RGMB), bone morphogenetic protein receptor type 2 (BMPR2), bone morphogenetic protein receptor type 1A (BMPR1A), bone morphogenetic protein 2 (BMP2), bone morphogenetic protein receptor type 1B (BMPR1B), chordin (CHRD), parathyroid hormone (PTH), parathyroid hormone 1 receptor (PTH1R), parathyroid hormone like hormone (PTHLH), and wnt family member 9A (WNT9A). To be specific, the size and color were influenced by the degree and the combined score dictating the edge size. Hence, low value corresponded to small size and light color. The result was shown in Fig. 5. BMP2 was identified as the key DEG between controls and recovery from FDM in our study.

Fig. 5

Protein-protein interaction analysis of differentially expressed genes in response to control vs FDM 0hrs and bone morphogenetic protein 2 (BMP2) was identified as the key gene.

Discussion

In this study, we sought to identify DEGs in retina between early and late stage FDM and then further analyzed functional annotation and potential biomarkers as a reliable therapy for myopia progression. We investigated that BEND6 was the only DEG although it was not statistically significant. This finding was inconsistent with the previous study that identified BEND6 causing reciprocal effects on neural stem cell renewal and neurogenesis [27]. Thus, we inferred that BEND6 may inhibit neuron genesis and signal connection in cells though considering no significant statistical meaning.

BMPs have a profound impact in bone formation and repair. Recent studies have underestimated their vital role in normal organs proliferation and vascular homeostasis in chronic diseases. Elevated expression amounts of BMPs have been associated with the commodities of diabetes mellitus. Increasing evidence from recent study demonstrates the involvement of BMP signaling in retinal inflammation, hyperpermeability and pathological neovascularization in diabetic retinopathy [28]. Previous studies have found that cellular signal transduction pathway and ion ligand activate pathway were related to proliferative vitreoretinopathy, proliferative diabetic retinopathy and wet age-related macular degeneration in zebrafish and mouse retina [29,30]. However, general pathway mechanism corresponding to FDM was still far from investigated.

Numerous mechanisms were detected to be associated with FDM through functional annotation of the DEGs. Among them, some have been reported to be related to myopia, including immunity [31,32], inflammation [[32], [33], [34]], and BMP2 upregulation [35,36]. We further analyzed the efficiency of BMP2 dysregulation on myopic retina by comparing early and late recovery from FDM. Because BMP2 dysregulation plays a pathogenic role in promoting complex growth-modulatory signal pathways. In early myopia, scleral fibroblast is a risk factor of disease progression, thus, our study reporting the influence of BMP2 on retina is vital to understand the scleral fibroblast function in myopia. In addition, we proposed that BMP2 was up-regulated in FDM 0hrs, and 4 DEGs were detected to skew myopic retina into a disease-promoting phenotype through scleral fibroblasts metabolism. BMP2 has been revealed to be up-regulated in FDM [36] and breast cancers [37]. Besides, it has been reported that BMP2 was involved with the adhesion of monocytes to fibronectin and endothelial cells [38], suggesting that BMP2 may contribute to the inflammatory system. However, the role of BMP2 in FDM and myopia pathogenesis has never been reported. Among the 14 mRNA target genes, BMP2, interacting with remaining six genes, was the most significant up-regulated gene. Meanwhile, our study observed that BMP2 disrupted the endothelial barrier function, induced endothelial-leukocyte adhesion and up-regulated retinal vascular endothelial growth factor levels, inflammatory markers and cytokines such as ICAM-1, IL-6 and IL-8 [28]. BMP2 has been found to be enriched in Smad pathway and up-regulated in myopia by Wang et al. [39]. They inferred that BMP2 may boost myopia by elevating the proliferation and differentiation of Human sclera fibroblast (HSF), as well as to benefit extracellular matrix synthesis potentially through classical pathway, which coincides with our results.

To the best of our knowledge, this is the first study to explore the DEGs by RNA sequencing analysis aiming to obtain deeply understanding of the pathophysiological mechanism of myopia. Besides, our study is of significance for other pathway-related diseases in adolescence such as refractive error and anisometropia. It should be noted that the present study has limitations. First, the relatively small sample size of this study may limit our ability. Therefore, we established many comparisons in multiple aspects to avoid possible bias. Secondly, our study was initially performed based on animal retina not human subjects. The main reason was little clinical data for human could be acquired in the present database. Last but not least, there is a lack of relevant experiments to identify the accuracy of bioinformatic results. To avoid potential risks, we implemented supreme criteria while selecting samples by setting |log 2fold change|>2, and FDR <5% in the GEO database. Further studies should be conducted on BMP2 to prove valuable functions in myopia prevention, as modulation of HSF may be a novel strategy for indulging myopia progression.

Conclusion

Taken together, the findings of this study suggest that cellular signal transduction plays a significant role in both myopic retina and sclera fibroblast about occurrence and progression of myopia. From the perspective of myopia prevention, it is an urgent method to inhibit BMP2 expression as well as proliferation of sclera fibroblast, which would induce risks of high myopia in myopic subjects, especially young children. To explore a supplementary bioinformatical signal pathway in BMP2 provides an original treatment of myopia progression that is based on molecular mechanism of myopia in clinical practice. In addition, BMP2 potentially represents a novel marker for myopia and may modulate extracellular matrix homeostasis in FDM.

Funding statement

This research did not receive a specific grant from any funding agency in the public, commercial, or not-for-profit sectors.

Declaration of competing interest

The authors declare no conflicts of interest regarding the publication of this article.