Role of the XIST-miR-181a-COL4A1 axis in the development and progression of keratoconus.

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.

Introduction

As a disease of the eye, keratoconus (KC) leads to the uninflammatory thinning of the corneal stroma [1]. KC usually occurs in the transition period from childhood to adulthood, and it is more common in Asian populations [2]. KC may induce the symptoms of double vision, nearsightedness, blurry vision, light sensitivity, and astigmatism [3]. KC patients may need special contact lenses or even corneal transplantation [4,5].

KC is caused by genetic, hormonal, and environmental factors [6], 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. [8] and Shetty et al. [10] 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 [14]. 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 [15]. 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 [18]. 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 [19]. 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.

Methods

Data preprocessing and differential expression analysis

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 [20]. The expression profile data related to read count level was normalized using R package preprocessCore (version 1.40.0, Berkeley, CA) [21].

The lncRNAs, miRNAs, and mRNAs in GSE112155 were annotated and identified in the HUGO Gene Nomenclature Committee (HGNC) database [22], which includes 3,979 lncRNAs, 1,932 miRNAs and 19,197 recognized protein-coding genes. Using R package edgeR (version 3.22.5) [23], 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.

Hierarchical clustering analysis and enrichment analysis

For the screened DE-RNAs, bidirectional hierarchical clustering based on a correlation algorithm was performed using the R package pheatmap (version 1.0.8) [24]. In addition, using the DAVID tool (version 6.8, Frederick, MD) [25], 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 [26], and through the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis [27]. A p value of <0.05 was taken as the threshold of enrichment significance.

Regulatory network analysis of ceRNA

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 [28]. The lncRNA–miRNA regulatory pairs in the miRecode (version 11, Gothenburg, Sweden) [29] and starBase (version 2.0, Guangzhou, China) [30] 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) [31], the lncRNA–miRNA regulatory network was visualized.

Using the starBase database [30], 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 [30]. 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 [31].

Along with the identified lncRNA–miRNA and miRNA–mRNA regulatory pairs, the ceRNA regulatory network was constructed using the Cytoscape software [31]. In addition, DAVID tool [25] 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.

Construction of KC-associated ceRNA regulatory network

Using “keratoconus” as the keyword, the KEGG pathways directly correlated with KC were searched in the Comparative Toxicogenomics Database (Raleigh, NC) [32]. 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.

Results

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.

Figure 1

The expression distribution curves, scatter diagram, clustering heatmap, and histogram of the differentially expressed RNAs (DE-RNAs). A: The expression distribution curves of the identified RNAs. The red, green, and blue curves represent long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and mRNAs, respectively. B: The scatter diagram of the DE-RNAs. The red, green, and blue dots represent significantly upregulated RNAs, significantly downregulated RNAs, and non–DE-RNAs, respectively. C: The clustering heatmap of the DE-RNAs. B: The histogram showing the proportion of upregulated and downregulated RNAs.

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).

Table 1

The results of enrichment analysis for the differentially expressed mRNAs.

Category	Term	Count	P-value	FDR
Biology Process	GO:0006811~ion transport	78	1.69E-11	3.01E-08
	GO:0007267~cell-cell signaling	62	1.49E-09	2.65E-06
	GO:0009611~response to wounding	55	1.30E-08	2.31E-05
	GO:0030001~metal ion transport	49	5.78E-08	1.03E-04
	GO:0006955~immune response	63	1.17E-07	2.08E-04
	GO:0006812~cation transport	54	1.34E-07	2.40E-04
	GO:0002526~acute inflammatory response	19	3.03E-07	5.39E-04
	GO:0006952~defense response	57	3.18E-07	5.67E-04
	GO:0006959~humoral immune response	17	3.43E-07	6.11E-04
	GO:0019226~transmission of nerve impulse	39	4.14E-07	7.37E-04
	GO:0007268~synaptic transmission	35	5.40E-07	9.62E-04
	GO:0006954~inflammatory response	36	1.44E-06	2.56E-03
	GO:0006956~complement activation	12	1.63E-06	2.90E-03
Cellular Component	GO:0005576~extracellular region	166	8.85E-13	1.22E-09
	GO:0005886~plasma membrane	264	2.56E-12	3.53E-09
	GO:0044421~extracellular region part	97	4.47E-12	6.16E-09
	GO:0044459~plasma membrane part	173	1.84E-11	2.54E-08
	GO:0031224~intrinsic to membrane	348	1.92E-11	2.65E-08
	GO:0005615~extracellular space	74	1.32E-10	1.83E-07
	GO:0016021~integral to membrane	334	2.43E-10	3.35E-07
	GO:0031226~intrinsic to plasma membrane	106	2.24E-09	3.10E-06
	GO:0005887~integral to plasma membrane	104	2.84E-09	3.92E-06
	GO:0034702~ion channel complex	27	6.91E-06	0.009529
Molecular Function	GO:0015267~channel activity	51	1.14E-10	1.74E-07
	GO:0005216~ion channel activity	49	1.17E-10	1.78E-07
	GO:0022803~passive transmembrane transporter activity	51	1.24E-10	1.88E-07
	GO:0022838~substrate specific channel activity	49	3.38E-10	5.15E-07
	GO:0022836~gated channel activity	40	4.97E-09	7.56E-06
	GO:0005261~cation channel activity	33	7.27E-07	1.11E-03
	GO:0005509~calcium ion binding	71	7.41E-06	1.13E-02
	GO:0046873~metal ion transmembrane transporter activity	34	1.19E-05	1.81E-02
	GO:0022834~ligand-gated channel activity	19	1.40E-05	2.13E-02
	GO:0015276~ligand-gated ion channel activity	19	1.40E-05	2.13E-02
	GO:0030594~neurotransmitter receptor activity	16	1.97E-05	2.99E-02
Pathway	hsa04080:Neuroactive ligand-receptor interaction	35	7.65E-08	1.07E-05
	hsa04610:Complement and coagulation cascades	15	5.01E-06	3.51E-04
	hsa04940:Type I diabetes mellitus	11	2.69E-05	1.26E-03
	hsa00982:Drug metabolism	13	3.97E-05	1.39E-03
	hsa00140:Steroid hormone biosynthesis	10	3.43E-04	9.55E-03
	hsa04060:Cytokine-cytokine receptor interaction	26	1.07E-03	2.47E-02
	hsa05330:Allograft rejection	8	1.64E-03	3.22E-02
	hsa04020:Calcium signaling pathway	19	2.53E-03	3.86E-02
	hsa00980:Metabolism of xenobiotics by cytochrome P450	10	2.51E-03	4.30E-02

Note: FDR, false discovery rate.

Analysis of the ceRNA regulatory network

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).

Figure 2

The long non-coding RNA (lncRNA)–microRNA (miRNA) regulatory network. The squares and diamonds represent lncRNAs and miRNAs, respectively. Blue and red represent upregulation and downregulation, respectively.

Figure 3

The microRNA (miRNA)–mRNA regulatory network. The diamonds and circles represent miRNAs and mRNAs, respectively. Blue and red represent upregulation and downregulation, respectively.

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).

Figure 4

The competing endogenous RNA (ceRNA) regulatory network. The squares, diamonds, and circles represent long non-coding RNAs (lncRNAs), microRNAs (miRNAs), and mRNAs, respectively. Blue and red represent upregulation and downregulation, respectively. Black lines and red lines represent, respectively, lncRNA–miRNA and miRNA–mRNA relationships.

Table 2

The pathways enriched for the mRNAs implicated in the competing endogenous RNA (ceRNA) regulatory network.

Term	Count	P-value	Genes
hsa04510:Focal adhesion	5	0.001561	COL4A2, COL4A1, CCND2, TNC, KDR
hsa04060:Cytokine-cytokine receptor interaction	5	0.003707	CXCL1, CCL22, CXCL5, KDR, CXCL10
hsa04062:Chemokine signaling pathway	4	0.006217	CXCL1, CCL22, CXCL5, CXCL10
hsa04512:ECM-receptor interaction	3	0.006331	COL4A2, COL4A1, TNC
*hsa04940:Type I diabetes mellitus	2	0.018766	CPE, HLA-B
*hsa05200:Pathways in cancer	4	0.021608	COL4A2, COL4A1, PPARG, RUNX1T1
hsa04623:Cytosolic DNA-sensing pathway	2	0.023854	IL33, CXCL10
hsa04270:Vascular smooth muscle contraction	2	0.042772	EDNRA, PRKG1

Note: “*” indicates the pathways overlapped with the keratoconus-correlated pathways searched from Comparative Toxicogenomics Database.

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.

Figure 5

The keratoconus (KC)-associated competing endogenous RNA (ceRNA) regulatory network. The squares, diamonds, and circles represent long non-coding RNAs (lncRNAs), miRNAs, and mRNAs, respectively. Blue and red represent upregulation and downregulation, respectively. Black lines and red lines represent, respectively, lncRNA–miRNA and miRNA–mRNA relationships. Purple triangles represent KC-associated pathways.

Figure 6

The expression levels of the RNAs correlated with keratoconus (KC). A: The expression levels of the mRNAs correlated with KC. B: The expression levels of the microRNAs (miRNAs) correlated with KC. C: The expression levels of the long non-coding RNAs (lncRNAs) correlated with KC. Black and white columns represent KC samples and control samples, respectively.

Discussion

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 [20]. 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) [38]. In addition, IL-6 assumes important roles in herpes simplex virus (HSV) type I infection-induced corneal nerve degeneration [39] and in ocular inflammation and angiogenesis in the cornea [40]. The roles of PPARG in inflammation have been widely reported [41-43]. PPARG regulates the redox balance in macrophages [41]. 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 [44]. HLA-A26, HLA-B40, and HLA-DR9 are frequently found in older Japanese populations and may be related to KC in young people [45]. HLA-G contributes to establishing immune tolerance in allograft, which may also help to maintain the immune-privileged status of the cornea [46]. 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 [51]. 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 [53]. 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 [54]. 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) [55]. 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 [56]. 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.