Non‐coding RNAs and related molecules associated with form‐deprivation myopia in mice

Abstract The role of miRNAs and its regulatory mechanism in myopia are indeterminate. Our study aimed to investigate potential myopia‐associated non‐coding RNAs and related molecules by performing a comprehensive bioinformatic analysis of miRNA expression profile of mice with form‐deprivation myopia (FDM). Differentially expressed miRNAs in two raw microarray data sets (GSE58124 and GSE84220) from Gene Expression Omnibus (GEO) database were comprehensively analysed using GEO2R. Target genes were predicted using miRDB and enriched with Metascape online tool. Protein‐protein interaction (PPI) networks were constructed utilizing STRING and Cytoscape. Significant differentially expressed miRNAs were validated by real‐time polymerase chain reaction (qRT‐PCR) using RNA extracted from monocular FDM ocular tissues. As result, we identified three upregulated miRNAs (mmu‐miR‐1936, mmu‐miR‐338‐5p, and mmu‐miR‐673‐3p) significantly associated with myopia in the two microarray data sets (p < 0.05 and |Log (Fold Change) |>1). GO functional analysis suggested these three miRNAs were targeted in genes mostly enriched in morphogenesis and developmental growth of retinal tissues. Enrichment analysis revealed top eight transcription factors, including PAX6 and Smad3, related to myopia. Ten hub genes, including Rbx1, Fbxl3, Fbxo27, Fbxl7, Fbxo4, Cul3, Cul2, Klhl5, Fbxl16 and Klhl42, associated with ubiquitin conjugation were identified. qRT‐PCR confirmed the increased expression of mmu‐miR‐1936 and mmu‐miR‐338‐5p (p < 0.05), but no statistical difference was observed in mmu‐miR‐673‐3p expression in myopic retinas. Our findings indicated mmu‐miR‐1936, mmu‐miR‐338‐5p and mmu‐miR‐673‐3p upregulation may be associated with myopia development via post‐transcriptional gene regulation, and identified potential molecules that could be further explored in future studies of the mechanism in myopia.


| INTRODUC TI ON
Myopia, characterized by a mismatch between eye shape and refractive ability, has become the major cause of worldwide severe vision loss with its rapidly increasing incidence rates. 1 The worldwide prevalence of myopia is estimated to reach 49.8% by 2050, affecting nearly half the world's population. 2 Myopia, especially high myopia, can also lead to glaucoma, cataracts, retinal detachment, macular degeneration and other vision-threatening complications, and even blindness in severe cases. 3,4 Therefore, it is necessary to clarify the detailed mechanism of myopia to identify potential targets of treatment and ultimately reduce incidence rates. [5][6][7] Current theories have acknowledged the vital role of scleral remodelling in the development of myopia, initiated through an environmental stimulus, a genetic stimulus or a combination of the two. [8][9][10][11] During the process of emmetropization, abnormal visual signals in the retina produce unknown molecules that pass to the choroid and then initiate responses in the sclera, leading to an abnormal metabolic process of scleral remodelling, characterized by the decelerated synthesis and accelerated degradation of extracellular matrix (ECM) components. [12][13][14][15] Several studies indicated that the regulation of gene expression is key to the changes in scleral remodelling caused by visual signals. 16 microRNAs (miRNAs) are small non-coding RNA that have assumed a high priority in recent studies due to their ability to modulate gene expression by targeting key epigenetic enzymes. 17,18 A previous study found that miR-328 can reduce the expression level of Pax6 by binding to the 3'-UTR, resulting in an increased expression of MMP2, which can enhance the degradation of extracellular matrix (ECM). 19,20 As a recognized ECM regulator, miR-29a has been found to affect the expression of collagen I through Wnt/β-catenin signalling pathway. 21 Other miRNAs, such as miR-184, miR-29a and miR-let-7i, have also been associated with myopia. However, the regulatory mechanism of miRNAs in myopia has not been explicitly described in previous studies. 18 The purpose of this study was to identify differentially expressed miRNAs involved in the development of form-deprivation myopia using two microarray data sets (GSE58124 and GSE84220) obtained from Gene Expression Omnibus database (GEO). Deep bioinformatics analysis, including target gene prediction, enrichment analysis, protein-protein network (PPI) analysis and identification of hub genes, were performed and identified differentially miRNAs were further validated in monocular mouse tissues with form-deprivation myopia using real-time polymerase chain reaction (qRT-PCR). The flow diagram of the analyses is shown in Figure 1.

| Non-coding RNA microarray profile collection
Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih. gov/geo) is an online public genome database that stores high-throughput gene expression data sets. 22 miRNA expression profiles examined in our study were retrieved from GEO by searching for the following key terms: ('microRNA' OR 'miRNA') AND 'form-deprivation myopia'.
The inclusion criteria for selecting the data sets for assessing the regulation mechanism of retina in myopic sclera were as follows: (i) profiles containing miRNA expression levels of retina samples, (ii) the myopia modelling was form deprivation instead of lensinduction, and (iii) the organism set was retained, with mouse assuming precedence.
The process resulted in the selection of two gene expression profiles (GSE58124 and GSE84220), separately based on GPL22134 and GPL15880 platforms (GPL15880, Agilent-029152 Unrestricted Mouse miRNA v15.0 microarray). GSE58124 consisted of six whole eye, three retina and three sclera samples, whereas GSE84220 contained three retina and three sclera samples. Both chips contained an equal amount of normal tissues, with the other eye as a control.
All miRNA information was downloaded for subsequent analysis. All raw data were then analysed by using GEO2R tools to identify differentially expressed miRNAs between the myopic eye and contralateral control eye. The cut-off criteria were p < 0.05 and |Log (Fold Change) | >1. The GSE58124 and GSE84220 data set were analysed separately, and a heat map was constructed using the limma package in the R Studio (version 1.3.959). Volcano plots were constructed using Prism (version 8.4.0), and overlapping miRNAs from the two data sets were represented on a Venn diagram. The intersecting portion of the Venn diagram was defined as differentially expressed miR-NAs, DE-miRs.
Except for retinas, miRNA expression level data of other ocular tissues (sclera and whole eye samples) were available. All miRNA information was analysed in the two GEO sets using GEO2R; however, only DE-miRs screened out of the retina samples were selected for subsequent analyses.

| Target gene prediction and enrichment analysis
The potential targets of DE-miRs were predicted using the miRDB (http://www.mirdb.org/) databases. 23,24 Metascape (http://metas cape.org/gp/index.html#/main/step1) is a simple and powerful tool for gene function annotation, analysis and visualization; it integrates many authoritative data resources such as Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), UniProt and Drugbank. 25 In this study, transcription factors, tissue-/cell-specific expression, GO and KEGG pathway analyses of target gene prediction were performed via Metascape.

| Protein-protein interaction (PPI) network creation and modular analysis
All predicted target genes of the DE-miRs were uploaded to STRING (version 11.0) to obtain the PPI pairs and subsequently imported into NetworkAnalyst (https://www.netwo rkana lyst.ca) for visualization. [26][27][28][29] Organism was limited to 'Mus musculus', and a minimum required interaction score >0.400 (medium confidence) was considered significant. Cytoscape (version 3.8.0) is an open source software platform for visualizing complex networks and integrating these with any type of attribute data. 29 After the PPI relation pairs were obtained, functional modules were figured out by cluster analysis using the Molecular Complex Detection (MCODE) plug-in of Cytoscape. Analysis parameters were as follows: clusters were found in whole network and 'haircut' was selected instead of 'fluff', degree cut was 2, node density cut-off was 0.1, node score cut-off was 0.2, k-core was 2, and maximum depth was 100.

| Identification of hub genes
After PPI network construction and cluster analysis based on MOCDE, hub genes were selected using another useful plug-in of Cytoscape, cytoHubba. cytoHubba is used to predict and explore important nodes and subnetworks in a given network based on several topological algorithms. 30 The constructed PPI network was uploaded to the cytoHubba plug-in of Cytoscape as a target network.
Subsequently, node scores were calculated by all twelve methods, and the top ten genes with the highest degree of connectivity (ranked by MCC) were selected as hub genes of the predicted DE-miR target genes. A miRNA-hub gene network was then constructed using Cytoscape. Biotech, Shanghai, China). The relative miRNA level was normalized to the internal control U6 using the 2 −∆∆CT method. p value < 0.05 was considered statistically significant.

| Identification of differentially expressed miRNAs
All raw data were uploaded to GEO2R to evaluate the differentially expressed miRNAs using the criterion: p value < 0.05 and |Log (Fold Change) | >1. From retina tissues in the data sets, sixteen differential expressed miRNAs (5 upregulated and 11 downregulated) were screened out from the GSE84220 data sets. A total of 115 miRNAs F I G U R E 2 Heat maps of differentially expressed miRNAs between a treated eye (myopia retina) and contralateral untreated eye (control retina) in two microarray profiles. (A) Heat map of GSE84220 data set; (B) heat map for GSE58124 data set (only 49 unique miRNAs are shown of a total of 115 microRNAs) F I G U R E 3 Volcano plots and Venn diagram of this study. (A, B) Volcano plots of the distribution of all differentially expressed miRNAs that were selected with a p value < 0.05 and |Log (Fold Change) | >1. Upregulated, downregulated and not significantly changed miRNAs are marked, respectively, as red, green and black dots. From the retina tissues in the data set, a total of 115 miRNAs were differentially expressed in the GSE58124 data sets, and all were upregulated. Sixteen differentially expressed miRNAs (5 upregulated and 11 downregulated) were screened out from the GSE84220 data sets. (C) Venn diagram of co-expressed miRNAs within GSE58124 and GSE84220, with consensus defined as DE-miRs. (D) mmu-miR-1936, mmu-miR-338-5p and mmu-miR-673-3p were screened out as DE-miRs; their specific p and logFC values are listed in the two data sets (49 unique miRNAs) were differential expressed in the GSE58124 data sets, and all were upregulated. The differentially expressed unique miRNAs of the two data sets are displayed by separate heat maps in Figure 2. Volcano plots ( Figure 3A,B) and a Venn diagram ( Figure 3C) were used to highlight the differential expressed miRNAs and the overlap of the two data sets respectively. Overall, we identified three differentially upregulated miRNAs, defined as DE-miRs: mmu-miR-1936, mmu-miR-338-5p and mmu-miR-673-3p ( Figure 3D).
In scleral tissues, no differentially expressed miRNAs were found in both data sets. However, in whole eye samples, we identified 39 upregulated miRNAs (Table S1), although the DE-miRs of these samples were not further analysed.

| Target gene prediction and enrichment analysis
A total of 1340 unique target genes were predicted by an online bioinformatics tool miRDB as presented in Table S2. GO Figure 4D).

| Protein-protein interaction (PPI) network analysis and modular analysis
The PPI network is the relationship network of biomolecules, which plays an important role in life movement. Identifying the functional modules in a complex interactome can therefore help us to understand the pathogenesis of myopia. STRING database was used to assess functional associations, and a PPI network was constructed to include 1330 nodes and 5644 edges, which was visualized using NetworkAnalyst ( Figure 5A). In the protein network graph, each node represents a protein and the edge represents a connection between two proteins. A total of 36 clusters were screened out using the MCODE app, and the top two clusters are presented in

| Hub gene analysis and miRNAs-hub gene network construction
Hub gene analysis identified ten genes that could be regulated by the three differentiated miRNAs ( Figure 6A Figure 6A, and the shortest path is shown in Figure 6B. The integrated result of all the twelve methods used to calculate scores in the cytoHubba plug-in of Cytoscape is shown in Table S3. The miRNA-hub gene network is displayed in Figure 6C. mmu-miR-338-5p potentially could target 8 (Rbx1, Fbxl3, Fbxo27, Cul3, Cul2, Klhl42, Klhl5 and Fbxl16) of the10 hub genes, whereas Fbxo4 and Fbxl7 were identified as potential target of mmu-miR-1936.
Interestingly, we found that the 10 hub genes were all related to ubiquitination. IUUCD (version 2.0), a family-based database of ubiquitin and ubiquitin-like conjugation, was used to obtain the latest information of functional roles of the 10 genes. 33 The placings, score, gene name, protein name and functional description are shown in Table S4.

| Upregulation of mmu-miR-1936 and mmu-miR-338-5p in retinal tissues of mice with formdeprivation myopia
We further verified the expression level of miRNAs using qRT-PCR analysis. As shown in Figure 7, the expression levels of mmu-miR-1936 and mmu-miR-338-5p in retinal tissues of myopic eyes were statistically significantly higher than those of controls, but there was no statistical difference in mmu-miR-673-3p levels.

| DISCUSS ION
In present study, we identified three upregulated miRNAs (mmu-miR-1936, mmu-miR-338-5p and mmu-miR-673-3p) in the retina of mice with form-deprivation myopia based on the results of overlapping two microarray data sets (GSE58124 and GSE84220) from the Gene Expression Omnibus (GEO) database. We verified the abnormal expression of two of the three miRNAs using qRT-PCR and subsequently performed a comprehensive bioinformatic analysis, providing a framework for understanding the molecular mechanisms underlying the role of these miRNAs in formdeprivation myopia.
In our study, we found three differential expressed miRNAs (mmu-miR-1936, mmu-miR-338-5p and mmu-miR-673-3p) associated with myopia. Currently, there are few studies on these three miRNAs, most of which have focused on cancer. [34][35][36] Previous studies found that miR-338-5p promoted glioma cell invasion by regulating matrix metalloproteinases 2 (MMP2), which has been implicated in the development of myopia. 37,38 Given the fact that miRNAs can mediate ECM remodelling via MMP2, a process that drives tumour invasion, it may be possible that these miRNAs influence myopic scleral ECM remodelling in a similar way during visual manipulation in form-deprivation myopia. We surmise that the retina is involved in the regulation of myopic sclera through DE-miRs; however, this process needs to be confirmed by further research.  [43][44][45] There are a few limitations in the current study. First, to determine actual differentially expressed miRNAs in myopic retina, a larger sample size would be needed for external verification. Moreover, changes in myopic retina miRNAs/target mRNAs were analysed in the present study. However, there was no direct evidence that the changes in biochemical substances or gene expression in the retina led to the change in RPE layer miRNAs, with consequent action on the sclera to cause myopia. It is worth mentioning that miRNAs can be packaged in structures called extracellular vesicles. These vesicles, including exosomes and microcapsules, are cell-derived membranous structures that are transferred between cells to establish intercellular communication. 46,47 Future studies can further clarify the mechanism of visual manipulation-induced myopia by labelling and tracing retinal exosomes that contain miRNAs.
In summary, this study identified key miRNA signalling factors and their interactions of myopia, expanding our current knowledge of the patterning mechanisms active throughout myopia formation and bringing new insights into the molecular mechanisms underlying form-deprivation myopia.

ACK N OWLED G M ENTS
This work was supported by the National Natural Science Foundation of China (Grant numbers 81800871).

CO N FLI C T O F I NTE R E S T
The authors declare no competing financial interests.

R E FE R E N C E S
F I G U R E 7 Validation of the relative expressed level of mmu-miR-1936, mmu-miR-338-5p and mmu-miR-673-3p. *p < 0.05. The fold changes were normalized to small nuclear RNA U6 and calculated using the 2 −∆∆CT method. The validation results show that the expression levels from myopia of mmu-miR-1936 and mmu-miR-338-5p in retinal tissues were statistically significantly higher than those of controls