Molecular Vision 2020; 26:1-13
Received 05 January 2019 | Accepted 31 January 2020 | Published 01 February 2020
Rui Tian, Lufei Wang, He Zou, Meijiao Song, Lu Liu, Hui Zhang
Department of Ophthalmology, the Second Hospital of Jilin University, Changchun 130000, Jilin Province, China
Correspondence to: Hui Zhang, Department of Ophthalmology, the Second Hospital of Jilin University, No. 218 Ziqiang Street, Nan Guan District, Changchun 130000, Jilin Province, China; Phone: +86-0430-81136545; FAX: +86-0430-81136545; email: email@example.com
Background: As a disorder occurs in the eyes, keratoconus (KC) is induced by the thinning of the corneal stroma. This study was designed to reveal the key long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and mRNAs involved in the mechanisms of KC.
Methods: Transcriptome RNA-seq data set GSE112155 was acquired from the Gene Expression Omnibus database, which contained 10 KC samples and 10 myopic control samples. Using the edgeR package, the differentially expressed (DE)-mRNAs between KC and control samples were screened. The DE-lncRNAs and DE-miRNAs in this data set were identified using the HUGO Gene Nomenclature Committee (HGNC). Using the pheatmap package, bidirectional hierarchical clustering of the DE-RNAs was conducted. Then, an enrichment analysis of the DE-mRNAs was performed using the DAVID tool. Moreover, a competitive endogenous RNA (ceRNA) regulatory network was built using the Cytoscape software. After KC-associated pathways were searched within the Comparative Toxicogenomics Database, a KC-associated ceRNA regulatory network was constructed.
Results: There were 282 DE-lncRNAs (192 upregulated and 90 downregulated), 40 DE-miRNAs (29 upregulated and 11 downregulated), and 910 DE-mRNAs (554 upregulated and 356 downregulated) between the KC and control samples. A total of 34 functional terms and 9 pathways were enriched for the DE-mRNAs. In addition, 6 mRNAs (including PPARG, HLA-B, COL4A1, and COL4A2), 5 miRNAs (including miR-181a), 9 lncRNAs (including XIST), and the XIST-miR-181a-COL4A1 axis were involved in the KC-associated ceRNA regulatory network.
Conclusions: PPARG, HLA-B, COL4A1, COL4A2, miR-181a, and XIST might be correlated with the development of KC. Further, the XIST-miR-181a-COL4A1 axis might be implicated in the pathogenesis of KC.
As a disease of the eye, keratoconus (KC) leads to the uninflammatory thinning of the corneal stroma . KC usually occurs in the transition period from childhood to adulthood, and it is more common in Asian populations . KC may induce the symptoms of double vision, nearsightedness, blurry vision, light sensitivity, and astigmatism . KC patients may need special contact lenses or even corneal transplantation [4,5].
KC is caused by genetic, hormonal, and environmental factors , and genetic and environmental factors contribute to the pathogenesis and development of KC. Previous studies have shown that the alteration of inflammatory factors and genetic molecules is associated with KC [7-9]. For instance, Pahuja et al.  and Shetty et al.  showed that the levels of inflammatory factors such as tumor-necrosis factor (TNF)-α, interleukin (IL)-6, and matrix metalloproteinase 9 (MMP-9) were elevated in the epithelium of KC patients. Long non-coding RNA (lncRNA), microRNAs (miRNAs), and several pathways—including transforming growth factor β (TGF-β), phosphatidylinositol 3-kinase (PI3K)/protein kinase B (Akt), and Wnt signaling pathways—may function in the development and progression of KC [7,9,11]. The expression of miR-184 in normal cornea samples is higher than that of miR-205, and miR-184 may act in cornea development and corneal diseases [12,13]. Calpastatin plays a role in genetic susceptibility to KC, and the differential modulation of calpain–calpastatin complex may influence the functional defect of the cornea . Downregulated β-actin may be correlated with decreased human antigen R (HuR) in the corneal stroma of KC patients, which may also serve as a risk factor for the occurrence and development of KC . However, the key RNAs involved in the pathogenesis of KC have not been comprehensively identified, and the pathogenesis of KC remains unclear.
In recent years, more and more evidence has shown that the mutual regulation models between lncRNA and miRNA and their downstream target genes are closely related to the occurrence and development of diseases [16,17]. As an important factor in post-transcriptional regulation, miRNA activity can be regulated by lncRNA through sponge adsorption . As competitive endogenous RNA (ceRNA), lncRNA competitively binds to miRNA to regulate the protein level of the coding gene and participates in the regulation of cell biologic behaviors . Therefore, in-depth research into the control mechanism of ceRNA will help us better understand the occurrence and development of diseases.
By analyzing the gene expression profile of KC, the differentially expressed (DE)-RNAs between KC samples and myopic control samples were screened. Following this, a ceRNA regulatory network was built to select the key RNAs affecting the development of KC from those in myopia. This study might help to reveal the ceRNA regulatory mechanisms correlated with the pathogenesis of KC.
Transcriptome RNA-seq data set GSE112155 was downloaded from the Gene Expression Omnibus (GEO) database. There were 20 cornea epithelial tissue samples in GSE112155, including 10 KC samples (10 males) and 10 myopic control samples (6 males and 4 females). Female control samples were included in our study to avoid a gender bias, as reported by You . The expression profile data related to read count level was normalized using R package preprocessCore (version 1.40.0, Berkeley, CA) .
The lncRNAs, miRNAs, and mRNAs in GSE112155 were annotated and identified in the HUGO Gene Nomenclature Committee (HGNC) database , which includes 3,979 lncRNAs, 1,932 miRNAs and 19,197 recognized protein-coding genes. Using R package edgeR (version 3.22.5) , the DE-lncRNAs, DE-miRNAs, and DE-mRNAs between the KC and control samples were analyzed. A |log2 fold change (FC)| of >1 and false discovery rate (FDR) of <0.05 were selected as the thresholds of significant differential expression.
For the screened DE-RNAs, bidirectional hierarchical clustering based on a correlation algorithm was performed using the R package pheatmap (version 1.0.8) . In addition, using the DAVID tool (version 6.8, Frederick, MD) , the DE-mRNAs were analyzed though the Gene Ontology (GO) enrichment analysis, which included the Biology Process (BP), Molecular Function (MF), and Cellular Component (CC) categories , and through the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis . A p value of <0.05 was taken as the threshold of enrichment significance.
As a kind of ceRNA within miRNA, lncRNA is involved in the expression regulation of target genes and plays an important role in the occurrence and development of diseases . The lncRNA–miRNA regulatory pairs in the miRecode (version 11, Gothenburg, Sweden)  and starBase (version 2.0, Guangzhou, China)  databases were merged, and the regulatory pairs between DE-lncRNAs and DE-miRNAs were then screened from them. Using the Cytoscape software (version 3.6.1, Bethesda, MD) , the lncRNA–miRNA regulatory network was visualized.
Using the starBase database , the target genes were predicted for the DE-miRNAs involved in the identified lncRNA–miRNA pairs. The starBase database integrates target prediction information from the targetScan, picTar, RNA22, PITA, and miRanda databases . In this study, miRNA-target regulatory pairs included in at least one of the five databases were selected. Subsequently, the DE-mRNAs were correlated with the predicted target genes, and only the miRNA-target regulatory pairs involved the DE-miRNAs and DE-mRNAs with opposite expression change directions were retained. Moreover, the miRNA–mRNA regulatory network was built using the Cytoscape software .
Along with the identified lncRNA–miRNA and miRNA–mRNA regulatory pairs, the ceRNA regulatory network was constructed using the Cytoscape software . In addition, DAVID tool  was used to enrich the KEGG pathways for the mRNAs implicated in the ceRNA regulatory network, with a p value of <0.05 as the significance threshold.
Using “keratoconus” as the keyword, the KEGG pathways directly correlated with KC were searched in the Comparative Toxicogenomics Database (Raleigh, NC) . To obtain the mRNAs participating in the KC-associated pathways, the searched pathways were compared with the pathways enriched for the mRNAs implicated in the ceRNA regulatory network. Finally, the regulatory pairs involving the mRNAs enriched in the KC-associated pathways were extracted from the ceRNA regulatory network to construct the KC-associated ceRNA regulatory network.
A total of 22,281 RNAs were detected in transcriptome RNA-seq data set GSE112155. Based on the HGNC database, 3,261 lncRNAs, 789 miRNAs, and 18,231 mRNAs were identified. The expression distribution curves of the identified lncRNAs, miRNAs, and mRNAs are shown in Figure 1A. There were 282 DE-lncRNAs (192 upregulated and 90 downregulated), 40 DE-miRNAs (29 upregulated and 11 downregulated), and 910 DE-mRNAs (554 upregulated and 356 downregulated) between the KC and control samples. A scatter diagram of the screened DE-RNAs is presented in Figure 1B.
A clustering heatmap showed that the expression values of the screened DE-RNAs could separate different types of samples (Figure 1C,D). The result suggested that the DE-RNAs had expression characteristics in the samples. Enrichment analysis showed that 13 GO_BP terms (such as ion transport), 10 GO_CC terms (such as extracellular region), 11 GO_MF terms (such as channel activity), and 9 KEGG pathways (such as neuroactive ligand-receptor interaction) were enriched for the DE-mRNAs (Table 1).
After 66 eligible lncRNA–miRNA pairs were obtained, the lncRNA–miRNA regulatory network (involving 31 lncRNAs and 8 miRNAs) was built (Figure 2). Afterwards, 116 miRNA–mRNA regulatory pairs were acquired to construct the miRNA–mRNA regulatory network (involving 7 miRNAs and 92 mRNAs; Figure 3).
Then, lncRNA–miRNA–mRNA regulatory pairs were then identified, and the ceRNA regulatory network was constructed (Figure 4). In the ceRNA regulatory network, there were 131 nodes (including 32 lncRNAs, 7 miRNAs, and 92 mRNAs) and 182 edges (including 66 lncRNA–miRNA pairs and 116 miRNA–mRNA pairs). Additionally, 8 KEGG pathways (such as focal adhesion) were enriched for the mRNAs implicated in the ceRNA regulatory network (Table 2).
Only 12 KEGG pathways directly correlated with KC were found in the Comparative Toxicogenomics Database. After comparing the searched pathways with the pathways enriched for the mRNAs implicated in the ceRNA regulatory network, two overlapping pathways (type I diabetes mellitus, and pathways in cancer) were obtained as KC-associated pathways. Subsequently, the KC-associated ceRNA regulatory network was built (Figure 5). Meanwhile, the RNAs correlated with KC were obtained, including 6 mRNAs (peroxisome proliferator-activated receptor gamma [PPARG]; human leukocyte antigen-B [HLA-B]; runt-related transcription factor 1 [RUNX1T1]; carboxypeptidase E [CPE]; collagen, type IV, alpha 1 [COL4A1]; and collagen, type IV, alpha 2 [COL4A2]), 5 miRNAs (miR-301a, miR-181a, miR-222, miR-98, and miR-128), and 9 lncRNAs (X-inactive specific transcript [XIST]; ST7-AS2; LINC00309; LINC00299; LINC00261; LINC00276; LINC00355; rhabdomyosarcoma 2-associated transcript [RMST]; LINC00520; and prostate cancer antigen 3 [PCA3]; Figure 6). The XIST-miR-181a-COL4A1 axis was particularly involved in the KC-associated ceRNA regulatory network.
Previous studies have suggested that sex and ethnicity may affect morbidity, gene expression, and episode age in KC [20,33,34]. Gender consistency in gene expression, however, will generally explain the pathogenesis of KC. In this study, a total of 282 DE-lncRNAs (192 upregulated and 90 downregulated), 40 DE-miRNAs (29 upregulated and 11 downregulated), and 910 DE-mRNAs (554 upregulated and 356 downregulated) were screened for KC, relative to control samples with no gender bias. For the DE-mRNAs, 13 GO_BP terms, 10 GO_CC terms, 11 GO_MF terms, and 9 KEGG pathways were enriched. In the KC-associated ceRNA regulatory network, there were 6 mRNAs (including PPARG, HLA-B, COL4A1, and COL4A2), 5 miRNAs (including miR-181a), and 9 lncRNAs (including XIST). The number of DE-mRNAs (910) between the KC and control samples was higher than 13 and lower than 1422 without and with gender bias, respectively, in You . The number of DE-mRNAs (910) with strict criteria (edgeR, |log2FC| >1, and FDR <0.05) suggested that the result of the gene expression profile was credible.
The elevated inflammatory factors, including IL-6, TNF-α, and IL-1 receptors, have been identified in KC epithelium, despite KC having been defined as a non-inflammatory condition [8,35-37]. Multiple studies have shown conflicting results about inflammation in KC. Among the reported inflammatory factors in KC, IL-6 and IL-6 receptor (IL-6R) play roles in glaucomatous optic nerve and retina damage, and their abnormal single nucleotide polymorphisms (SNPs) are involved in the development and progression of primary open-angle glaucoma (POAG) . In addition, IL-6 assumes important roles in herpes simplex virus (HSV) type I infection-induced corneal nerve degeneration  and in ocular inflammation and angiogenesis in the cornea . The roles of PPARG in inflammation have been widely reported [41-43]. PPARG regulates the redox balance in macrophages . A previous study explored the functions of PPARG agonist rosiglitazone on retinoblastoma cells and found that rosiglitazone plays an antitumor role via the suppression of cell growth, metastasis, and invasion and via the promotion of cell apoptosis . HLA-A26, HLA-B40, and HLA-DR9 are frequently found in older Japanese populations and may be related to KC in young people . HLA-G contributes to establishing immune tolerance in allograft, which may also help to maintain the immune-privileged status of the cornea . Our present study, however, identified the downregulation of IL-6, PPARG and HLA in epithelium from patients with KC. We speculated the participation of IL-6, PPARG and HLA in KC pathogenesis. However, there were conflicting results between our study and others reporting the elevation of IL-6 in KC epithelium compared with controls [8,35-37].
In addition to conflicting inflammatory conditions in KC, histopathological changes in collagen decomposition or fibrosis are associated with KC induction [47,48]. Transcription factor 8 (TCF8) plays a role in about half of posterior polymorphous corneal dystrophy (PPCD) cases, and its target, COL4A3, is critical in both Alport syndrome and PPCD [49,50]. Enhancement during collagen decomposition is responsible for the damage of pathological tissues in KC corneas, which is still retained in initial cultures of KC fibroblasts . In KC corneas, type I, III, and V collagens have no difference in distribution, while the distribution of type IV collagen is disruptive and excrescent in the corneal basement membrane [47,52]. KC corneas have decreased collagen protein levels, and collagen type IV functions as a candidate gene in the development of KC . Our present study suggested that COL4A1 and COL4A2 were dysregulated in the KC epithelium relative to controls, suggesting the crucial role of collagen decomposition in the progression of KC.
In addition to the DE-mRNAs, we also identified the DE-miRNAs and lncRNAs in KC samples compared with controls. The miR-181a, miR-21, and Smad signaling coordinately regulate the expression of TGF-β-induced gene (TGFBI) protein (TGFBIp) in corneal fibroblasts, and their pharmacologic modulation may be applied to treat TGFBI-correlated corneal dystrophy . The overexpression of miR-181b can be induced by hypoxia, which further promotes the angiogenesis of retinoblastoma cells by mediating GATA-binding protein 6 (GATA6) and programmed cell death 10 (PDCD10) . Upregulated lncRNA XIST is positively related to an advanced stage and late differentiation state of retinoblastoma, and XIST may accelerate retinoblastoma progression via regulating the miR-124/signal transducer and activator of transcription 3 (STAT3) axis . The XIST-miR-181a-COL4A1 axis was involved in the KC-associated ceRNA regulatory network, indicating that XIST and miR-181a might be correlated with the pathogenesis of KC through the XIST-miR-181a-COL4A1 axis.
In conclusion, 282 DE-lncRNAs, 40 DE-miRNAs, and 910 DE-mRNAs were identified between the KC and control samples. These DE-RNAs were identified without gender bias. Further, PPARG, HLA-B, COL4A1, COL4A2, miR-181a, and XIST might be involved in the development and progression of KC in both females and males. Moreover, the XIST-miR-181a-COL4A1 axis might function in the mechanisms of KC. However, the specific roles of these RNAs in KC should be further explored and supported by experiments.
This study is by the Norman Bethune Program of Jilin University (2015327); the Postdoctoral research project of Jilin province (801171060413); the Health service capacity improvement project (2017) of Jilin province (3D517EC63429); the youth program of National natural science foundation (2018) of China (81802998); 2019 basic Research Natural Science Foundation (20190201150JC). Conception and design of the research: R. Tian and H. Zhang. Acquisition, analysis and interpretation of data: R. Tian, H. Zou, L.F. Wang and L. Liu. Drafting the manuscript: R. Tian. Revision of manuscript for important intellectual content: H. Zhang, L. Liu and M.J. Song. All authors have read and approved the manuscript. Compliance with Ethical Standards: This article does not contain any studies with human participants or animals performed by any of the authors, the authors have no ethical conflicts to disclose.