Differential Expression of Coding and Long Noncoding RNAs in Keratoconus-Affected Corneas

Abstract

Purpose: Keratoconus (KC) is the most common corneal ectasia. We aimed to determine the differential expression of coding and long noncoding RNAs (lncRNAs) in human corneas affected with KC.

Methods: From the corneas of 10 KC patients and 8 non-KC healthy controls, 200 ng total RNA was used to prepare sequencing libraries with the SMARTer Stranded RNA-Seq kit after ribosomal RNA depletion, followed by paired-end 50-bp sequencing with Illumina Sequencer. Differential analysis was done using TopHat/Cufflinks with a gene file from Ensembl and a lncRNA file from NONCODE. Pathway analysis was performed using WebGestalt. Using the expression level of differentially expressed coding and noncoding RNAs in each sample, we correlated their expression levels in KC and controls separately and identified significantly different correlations in KC against controls followed by visualization using Cytoscape.

Results: Using |fold change| ≥ 2 and a false discovery rate ≤ 0.05, we identified 436 coding RNAs and 584 lncRNAs with differential expression in the KC-affected corneas. Pathway analysis indicated the enrichment of genes involved in extracellular matrix, protein binding, glycosaminoglycan binding, and cell migration. Our correlation analysis identified 296 pairs of significant KC-specific correlations containing 117 coding genes enriched in functions related to cell migration/motility, extracellular space, cytokine response, and cell adhesion. Our study highlighted the potential roles of several genes (CTGF, SFRP1, AQP5, lnc-WNT4-2:1, and lnc-ALDH3A2-2:1) and pathways (TGF-β, WNT signaling, and PI3K/AKT pathways) in KC pathogenesis.

Conclusions: Our RNA-Seq–based differential expression and correlation analyses have identified many potential KC contributing coding and noncoding RNAs.

Keratoconus (KC) is the most common corneal ectatic disorder.^1,2 KC is characterized by bulging and thinning cone-shaped cornea, resulting in myopia, irregular astigmatism, and, eventually, visual impairment. KC usually starts at puberty and progresses until the third or fourth decade.^3,4 KC is a complex multifactorial disorder with the contribution from both genetic and environmental factors as well as biomechanical, metabolic, and hormonal factors.^5–10 KC patients report significantly impaired vision-related quality of life, which worsens with time due to disease progression.¹¹ 

RNA sequencing (RNA-Seq) allows the entire transcriptome of coding and/or noncoding RNAs to be surveyed in a high-throughput and quantitative manner.^12–15 Long noncoding RNAs (lncRNAs) are class of noncoding RNAs that are at least 200 nucleotides long,^16,17 and have been found to regulate gene expression transcriptionally and post-transcriptionally in both physiologic and pathologic conditions.^17,18 

Previous studies have attempted to understand processes that underlie KC using different approaches. Wentz-Hunter et al.¹⁹ have performed targeted expression analysis of the stromal layer of the cornea to identify differentially expressed transcripts in KC versus a healthy cornea. Soon afterward, microarray was used to assess RNAs isolated from either the KC-affected epithelium or keratocytes.^20,21 By 2005, differentially expressed proteins in KC patients have been investigated using the whole cornea, certain layer, or tears.^22–25 In 2017, Szcześniak et al.²⁶ and Kabza et al.²⁷ used RNA-Seq to determine differentially expressed RNAs between KC-affected and other disease affected corneas. 

For the first time, we used Total RNA-Seq to identify differentially expressed coding/noncoding RNAs in KC-affected corneas versus nondiseased healthy ones, and to determine the expression correlation. We used droplet digital PCR (ddPCR) to validate the differential expression of selected RNAs. These differentially expressed RNAs may provide a better understanding of the functions of related lncRNAs and the potential pathways involved in KC pathogenesis. 

Methods

Human Cornea Samples and RNA Extraction

We followed the Tenets of Declaration of Helsinki. Unaffected human cornea samples were collected from postmortem donors from North Carolina Eye Bank, and the KC-affected corneal samples were collected during corneal transplantation surgery at the Cedars-Sinai Medical Center (CSMC). The institutional research board offices at Duke University School of Medicine and CSMC approved both studies, respectively. Written informed consent was obtained from all CSMC patients. This study included 10 KC-affected corneas and eight healthy non-KC corneas (Table 1). The inclusion and exclusion criteria for KC diagnosis have been described previously.^28,29 The diagnosis of KC was done by a cornea fellowship–trained ophthalmologist based on clinical examination and confirmed with videokeratography. Clinical examination included slit-lamp biomicroscopy, cycloplegic retinoscopy, and fundus evaluations. The slit-lamp biomicroscope was used to identify stromal thinning, Vogts' striae, and Fleischer ring. The retinoscopy examination was performed with a fully dilated pupil (20 minutes after phenylephrine 2.5% and cyclopentolate 1% drops had been instilled in the eye) to determine the presence or absence of retro illumination signs of KC, such as the oil droplet sign and scissoring of the red reflex. Videokeratography evaluation was performed on each eye using the Topographic Modeling System (TOMEY, Nagoya, Japan). Subjects were assigned as having KC if they had at least one clinical sign and a confirmatory videokeratography map with an asymmetric bowtie with skewed radial axis above and below the horizontal meridian (AB/SRAX) pattern.³⁰ KC grade was assigned based on the average central keratometry readings as follows: mild: <50 diopters (D); moderate: >50 D but < 55D; moderate-severe: 55 to 60 D; severe: >60 D. 

Table 1

 

Clinical Phenotypes of Patients With KC and Healthy Controls

 

Only KC patients at the advanced stage of KC from whom all other therapeutic interventions have failed were included in this study. They were diagnosed with KC grade of above 55 D (based on the average central keratometry readings). Half of the corneal button removed from surgery was used for RNA extraction while the other half was sent for pathologic examination. The postmortem donors had no other ocular disorders affecting their cornea. These cornea samples were excluded from usage in corneal transplant mostly due to the lack or incomplete medical/health information or occasionally low number of corneal endothelial cells. The central button of the donor corneas was acquired using a 6-mm biopsy punch (Cat#33-38, Miltex, Inc., York, PA, USA) and the whole button was used for RNA isolation. The cornea was stored in RNALater solution (Qiagen, Valencia, CA, USA) to preserve RNA quality. Total RNA was extracted using mirVana miRNA isolation kit from ThermoFisher Scientific, and with NucleoSpin RNA/Protein isolation kit (MACHEREY-NAGEL, Inc., Bethlehem, PA, USA), according to the manufacturer's recommendations. RNA quality was evaluated using Agilent 2100 Bioanalyzer with Agilent's RNA 6000 Pico Kit (Agilent, Santa Clara, CA, USA). We selected RNA samples with a minimum RIN (RNA Integrity Number) score of 6 for RNA-Seq. 

Stranded Total RNA-Seq

We used 200 ng total RNA from all subjects to prepare sequencing libraries with the TaKaRa SMARTer Stranded RNA-Seq kit after ribosomal RNA depletion using the RiboGone - Mammalian kit from TaKaRa (TaKaRa Bio USA, Inc., Mountain View, CA, USA). Briefly, the RiboGone-Mammalian kit removes ribosomal RNA (rRNA) and mitochondrial RNA (mtRNA) sequences from human total RNA samples based on hybridization and RNase H digestion that specifically depletes 5S, 5.8S, 18S, and 28S nuclear rRNA sequences, as well as 12S mtRNA sequences. The rRNA-depleted total RNA was reverse-transcribed to synthesize the first-strand complementary (c) DNA followed by PCR amplification using universal forward PCR primer and reverse PCR indexing primer set. Purified RNA-Seq library was validated using the Agilent 2100 Bioanalyzer with Agilent's High Sensitivity DNA Kit (Agilent, Santa Clara, CA, USA). Pooled RNA-Seq libraries were sequenced with paired-end 50-bp reads using an Illumina HiSeq 2500 sequencer (Illumina, Inc., San Diego, CA, USA) at the Georgia Cancer Center Integrated Genomics Core of Augusta University. 

After quality check and control with all sequencing reads, demultiplexed reads were aligned by TopHat³¹ using paired-end reading with the approximation of the median library size. Counts of sequencing reads were normalized using Cufflinks in fragments per kilo bases and millions reads (FPKM).³² After normalization, we annotated the transcripts with a gene transcript file from Ensembl database at the gene isoform level. We performed differential expression analysis using the Cuffdiff package.³³ We used a fold change of ≥ 2 and false discovery rate (FDR) ≤ 0.05 to identify genes with significant expression changes in KC corneas versus healthy ones. The list of differentially expressed genes was uploaded to the Web-based Gene Set Analysis Toolkit (WebGestalt; available in the public domain, www.webgestalt.org) 2013 version for gene ontology and KEGG pathway analysis.³⁴ For lncRNAs, we annotated them using a file from the NONCODE version 5 database.³⁵ After normalization using Cufflinks, we performed differential expression analysis using Cuffdiff. For lncRNA data, missing expression data in specific samples were replaced with a background value of 0.001 to enable the differential analysis between cases and controls. Without such replacement of missing data, the related lncRNAs would have been excluded from the differential expression. For lncRNA differential expression analysis, to avoid the possible impact of potential DNA contamination and RNA isolation bias, we required a minimal size of 300 bp for all the lncRNAs. We used |fold change| ≥ 2 and FDR ≤ 0.05 to identify lncRNAs with significant expression changes in KC corneas versus healthy ones. 

Validation Using Droplet Digital PCR Technique

Twenty-one differentially expressed coding RNAs were chosen for validation with ddPCR as previously described.^36,37 We selected these coding RNAs based on the following criteria: (1) reasonably high corneal expression with ≥ 10 normalized read counts; (2) relevance to corneal integrity, transparency, and cellular adhesion; (3) potential involvement in oxidative stress; and (4) related with corneal wound healing. Specific EvaGreen-based ddPCR assays were obtained from Bio-Rad Laboratories, Inc. (Hercules, CA, USA). To validate the differentially expressed lncRNAs, we selected seven lncRNAs with reasonable corneal expression (≥10 normalized read counts) and strong correlations with various other coding and lncRNAs. Specific ddPCR expression assays with specific probes were designed and ordered using the online ddPCR assay design tool with the specific RNA sequence of each selected lncRNA for copy number determination (available in the public domain, https://www.bio-rad.com/digital-assays/#/). We used RNA samples from the Total RNA-Seq experiment for validation: eight KC-affected corneas and eight control corneas. The ddPCR assays were performed using the QX200 Droplet Digital PCR system (Bio-Rad, Hercules, CA, USA) and its associated QuantaSoft Software, as previously described.^28,29 We used four internal reference genes (GAPDH, PPIA, HPRT1, and RPLP0)³⁸ to be quantified in the samples for normalization purpose. All four reference genes were stably expressed across the KC and control cornea samples and had similar performance in expression normalization, leading to almost identical successful validation of selected target genes. Thus, we decided to normalize with GAPDH for the figure purpose. Mann-Whitney U test was used to analyze the differential expression data from ddPCR for the statistical analysis. 

Expression Correlation and Network Analysis

The expression levels of differentially expressed coding RNAs and lncRNAs in all 18 corneal samples were imported into RStudio Version 0.98.1062 (RStudio, Boston, MA, USA)³⁹ from CSV raw data files. The combined expression data was structured into a table where the rows were the mapped RNA names and the columns were the case or control individual de-identified with a code name. Repeated gene names were removed and genes with no expression levels in all controls or cases were also removed. The correlation between all RNA expression levels was determined for the controls data and the cases data separately using R's built-in Pearson correlation function (cor). Using R psych package's test of significance for correlations (r.test), the P value comparing the respective mapped gene's control and case correlations were calculated. The FDR was calculated from the P values using the R's built-in Benjamini–Hochberg procedure function (p.adjust with “BH” parameter). We removed redundant correlations between each gene and itself as well as duplicated gene correlations (A-B is the same as B-A). 

The list of expression correlations was filtered with the criteria of an FDR ≤ 0.05 and having at least a mean expression level ≥ 20 in the cases or the controls with one gene. The entire filtered list of correlated RNAs including lncRNAs was visualized using Cytoscape version 3.5.1 (Cytoscape Consortium, New York, NY, USA).⁴⁰ Expression network hubs were defined as network nodes with multiple correlations to other RNAs. It was expected that this correlation expression analysis among coding RNAs and lncRNAs would help identify biologically important lncRNAs in human cornea. The list of correlated coding RNAs was loaded to WebGestalt 2013 version for gene ontology and KEGG pathway analyses.³⁴ 

Results

Identification of Differentially Expressed Coding RNAs

Each sample generated 20 to 44 million paired reads with 100% sample index match, > 96% bases with Q30, and mean quality score > 38. Overall, 62% to 88% of sequencing reads for all RNA samples were successfully mapped to the transcriptome file from Ensembl database, containing a total of 212,368 transcript isoforms. There were a total 95,954 and 103,440 expressed transcripts in all the healthy controls and KC-affected cornea samples, respectively. The number of expressed transcripts in each sample ranged from 120,000 to 139,000. Using a cutoff of at least one count of normalized transcripts, the number of expressed transcripts in each sample ranged from 35,000 to 45,000. There is a total of 14,111 and 19,608 expressed transcripts (≥1 count) in all the healthy controls and KC-affected corneal samples, respectively. The mean expression levels of these shared transcripts (≥1 count) were 142 reads in control corneas and 63 read counts in the keratoconic corneas. 

Using the cutoff with |fold change| ≥ 2 and FDR ≤ 0.05, we identified 436 differentially expressed coding RNAs (159 upregulated and 277 downregulated; Supplementary Table S1). Our analysis was consistent with previous expression studies by identifying differential expression of many known KC-related genes, including but not limited to TGFBR3, TIMP1, FBLN1, TIMP3, AQP5, and SFRP1. The top 20 up- and downregulated coding genes, showing the highest fold changes, were listed in Table 2. 

Table 2

 

Top Coding Genes With Significant Differential Expression in KC-Affected Corneas

 

Using WebGestalt with the differentially expressed coding RNAs, our gene ontology analysis revealed disruption in different biological processes (cellular proliferation, differentiation, locomotion, migration, cellular stress, and wound response), molecular functions (glycosaminoglycan binding, iron ion binding, antioxidant activity, and insulin-like growth factor binding), and cellular components (adherens junctions, extracellular matrix region, and basement membrane) in the keratoconic corneas. 

Validation of Differentially Expressed Coding RNAs

Out of 436 differentially expressed coding genes, we selected 21 genes for validation by ddPCR based on the criteria described earlier (Table 3). Of these 21 genes, 19 were confirmed to have significant and similar expression trend to the RNA-Seq results (Figs. 1A, 1B). The other two genes (APOD and IFI6) demonstrated similar expression trends but did not reach statistical significance. 

Table 3

 

Twenty-One Differentially Expressed Coding Genes to be Validated Using ddPCR

 

Figure 1

 

Differentially expressed coding genes in KC-affected human corneas identified by RNA sequencing (A) and ddPCR (B). For panel A, **FDR < 0.01, *0.01< FDR < 0.05. For panel B, ***P value < 0.001, **P value < 0.01, * 0.01< P value < 0.05. Error bars: panels A and B represent standard error.

Identification of Differentially Expressed lncRNAs

A total of 15,966 lncRNAs were expressed in the control and keratoconic cornea samples. Using the cutoff with |fold change| ≥ 2 and FDR ≤ 0.05, we identified 586 differentially expressed lncRNAs (343 upregulated and 243 downregulated; Supplementary Table S2). Next, differentially expressed coding RNAs and lncRNAs were used to perform a correlation analysis. 

Validation of Differentially Expressed lncRNAs

Based on their high cornea expression (≥10 normalized read counts), functional annotation in relation to KC-related pathways, and expression correlation with coding and lncRNAs, we selected seven differentially expressed lncRNAs for ddPCR validation (Table 4). Specific probe-based ddPCR assays with fluorescein dye (FAM) labeling were designed to validate the differential expression of these lncRNAs in KC-affected corneas versus the controls (Fig. 2A). The four selected internal reference genes were also used to normalize the expression of these lncRNAs. Five lncRNA showed a similar trend to RNA-Seq results (lnc-WNT4-2:1, MEG3, lnc-BLID-5:1, lnc-SCP2-2:1, and lnc-ALDH3A2-2:1). However, NONHSAT174656 and lnc-DLG-1:1 did not show any significant differential expression in KC-affected corneas versus controls (Figs. 2A, 2B). 

Table 4

 

Seven Differentially Expressed lncRNAs to be Validated Using ddPCR

 

Figure 2

 

Differentially expressed lncRNAs in KC-affected human corneas identified using by RNA-Seq (A) and ddPCR (B). For panel A, **FDR < 0.01, *0.01 < FDR < 0.05. For panel B, ***P value < 0.001, **P value < 0.01, *0.01< P value < 0.05. Error bars: panels A and B represent standard error.

Correlation Analysis of Coding and Noncoding RNAs

The expression correlation analysis identified a total of 492,135 pairs of RNAs with expression correlations in cases or controls. After comparing the correlations between cases and controls, a total of 352 significant differential correlations were identified with a FDR ≤ 0.05. After applying a minimum mean expression of 20 for at least one of the two paired transcripts in cases or controls, we identified 296 pairs of significant differential correlations between cases and controls (Supplementary Table S3). Many coding and noncoding RNAs were identified in multiple correlation pairs. These included A2M, APOD, MUC16, lnc-SCP2-2:1, lnc-DLG5-1:1, lnc-BLID-5:1, lnc-ALDH3A2-2:1, NONHSAT175650_1, and NONHSAT217671 (Fig. 3). We used Cytoscape to visualize the whole expression network of these 296 pairs (Supplementary Fig. S1). This network analysis revealed a number of genes with significant correlations in the keratoconic corneas versus the unaffected healthy controls (Table 5). We identified 117 coding RNAs in the significantly differential expression correlations. These 117 genes were enriched for their functions in negative regulation of cell migration and cell motility, extracellular space, cellular response to cytokine stimulus, cell surface receptor signaling pathway, and regulation of cell adhesion. 

Figure 3

 

Correlation networks of the differentially expressed coding RNAs and lncRNAs displayed by Cytoscape. The nodes represent mapped RNAs colored and sized by the degree of correlations with other RNAs. The edges are sized and colored based on the significant differences (i.e., FDR) between controls and cases. Correlation networks were shown for the Lnc-ALDH3A2-2:1 (A), APOD (B), CLIC4 (C), SLC2A1 (D), and A2M (E) genes.

Table 5

 

Coding and lncRNAs With Significant Differences in Expression Correlation (With a Reasonable Corneal Expressions) in KC-Affected Corneas Versus the Controls

 

Discussion

KC is a complex disorder with both genetic and environmental factors. Epithelial and stromal layers of the cornea show the prominent histopathologic changes during KC progression. A better understanding of the key molecular players involved in KC-associated corneal pathogenesis may provide potential therapeutic targets for disease treatment. Researchers used various techniques including RNA-Seq, microarray, RT-PCR, mass spectrometry, and Western blots.^{20,21,25,27,41,42} Recently, Kabza et al.²⁷ have used RNA-Seq for KC-affected corneas; however, this study was complicated by the use of other diseased corneas as controls. Compared with this study, our data have shown overlapping genes with differential expression (35 up- and 135 downregulated genes²⁷; Supplementary Table S4). 

We identified decreased expression in the cytokeratin (KRT6A, KRT13, KRT15, and KRT19) and collagen (COL8A1 and COL8A2) genes, indicating compromised integrity of KC cornea. Type VIII collagens are normally expressed in the Descemet's membrane⁴³ and targeted inactivation of Col8a1/Col8a2 leads to corneal abnormalities in mice.⁴⁴ 

Our study confirms the potential role of oxidative stress in KC-affected corneas by identifying differential expression of many genes involved in response to oxidative stress including GSTO1, SGK1, APOD, KLF9, GPX3, STK39, TXNIP, and lnc-ALDH3A2-2:1.^45,46 The majority of these genes in defense mechanism for oxidative stress were downregulated. However, lnc-ALDH3A2-2:1 is found to have more than 3-fold increased expression and correlated to other lncRNAs (Fig. 3A). Lnc-ALDH3A2-2:1, according to Lncipedia database,⁴⁷ is classified as an antisense transcript that overlaps with the protein coding gene ALDH3A2 on the opposite strand. However, using BLASTN, we found that lnc-ALDH3A2-2:1 sequence overlaps with the last exon of ALDH3A1 gene. The ALDH3A1 transcript did not show significant differential expression in our study; however, lncRNAs may alter protein level without affecting mRNA.⁴⁸ Another factor impacting oxidative stress is disruption of iron homeostasis, which can lead to various corneal diseases including KC.^23,24,49,50 Fleisher ring, a ring of iron deposition in the corneal cone, is a typical clinical sign for KC.⁵¹ Importantly, transferrin receptor (TFRC), a cellular receptor for iron-ferritin complex uptake,⁵² has been identified to have elevated expression in KC cornea in our study. This elevated expression in KC-affected corneas may result from iron deficiency.⁵³ We also found downregulation of the Ferritin Light Chain transcript (FTL) in KC corneas. Significantly, mutations in FTL have been linked to iron accumulation in the brain, causing neurodegenerative disease.⁵⁴ 

We observed differential expression of many genes involved in corneal wound healing, such as EPPK1, AQP3, AQP5, SLC2A1, and CTGF. Corneal wound healing is a complex process controlled by the integrated actions of multiple growth factors, cytokines, and proteases produced by epithelial cells, stromal keratocytes, and inflammatory cells.⁵⁵ However, unregulated healing processes, due to genetic or epigenetic factors, can lead to impairment in cellular function.^56,57 Macé et al.⁴¹ have shown antiproliferative and hyperapoptotic phenotype in KC-affected cornea, which also reveals a defect in wound healing process. Aquaporins (AQPs) are a ubiquitous family of transmembrane water channels to maintain tear film osmolarity, stromal layer thickness, and corneal transparency.⁵⁸ AQP3 and AQP5, important for corneal wound healing, are highly expressed in the cornea.⁵⁸ Both AQPs show reduced expressions in the KC-affected corneas (Table 2). Consistent with our results, Rabinowitz et al.⁵⁹ have reported the lack of AQP5 expression in KC-affected corneas. Although Garfias et al.⁶⁰ did not report this altered expression, this discrepancy could be due to their heterogeneous controls (normal corneal epithelium, sclerocorneal rims, and healthy corneal buttons). EPPK1 encodes a cytoskeletal linker protein called epiplakin, which maintains the integrity of intermediate filaments networks in epithelial cells.^61,62 Migration-dependent wound healing improvement has been found in epiplakin-null epithelium.⁶³ Upregulation of EPPK1 in our study suggests the potential involvement of EPPK1 in the impaired corneal wound healing in KC. Interestingly, solute carrier family 2 member 1 (SLC2A1) is upregulated during corneal epithelium wound healing.⁶⁴ Reduced expression of SLC2A1 in KC-affected corneas may reflect a lesser glucose uptake. This is intriguing because there is a reduced risk of developing KC in diabetic patients.^65–68 Moreover, it has been reported that diabetic patients have thicker cornea than healthy individuals.⁶⁴ 

Connective tissue growth factor (CTGF) is expressed in multiple ocular tissues as a downstream signaling molecule in the TGF-β pathway.⁶⁹ Importantly, loss of CTGF impairs efficient corneal wound healing.⁷⁰ The decreased expression of CTGF in our study is consistent with previous studies.^27,28 Dysregulation of extracellular matrix (ECM) pathway has been reported with KC progression.^28,51 Previous studies have shown an upregulation of ECM degradative enzymes and/or downregulation of their inhibitors.^27,45,71 Our data show significant downregulation in proteinases inhibitors (TIMP-1, TIMP-3, and A2M, and SLPI), but interestingly no significant change with ECM proteinases. 

WNT signaling is essential for normal corneal development. Sequence variants in WNT7B and WNT10A genes have been associated with CCT and KC risk.^72,73 A recent RNA-Seq–based study showed dysregulation in WNT signaling in KC-affected corneal epithelium.⁴² We found some evidence of WNT signaling disruption represented by downregulation of ANGPTL7, TGFBR3, SFRP1 and upregulation of frizzled class receptor 8 (FZD8) and early growth response 1 (EGR1). ANGPTL7 encodes angiopoietin-like 7 and is abundantly expressed in keratocytes to maintain corneal avascularity and transparency.⁷⁴ SFRP1 inhibits WNT signaling by preventing the formation of the WNT-FZD complex. Iqbal et al.⁷⁵ have shown an increased expression of SFRP1 in the KC-affected corneal epithelium while SFRP1 expression is reduced in tears from KC patients.⁷⁶ TGFBR3 suppresses WNT signaling through binding to Wnt3a and β-catenin.⁷⁷ We observed increased expression of Lnc-WNT4-2:1 in KC versus controls. Our BLASTN analysis indicated that Lnc-WNT4-2:1 is a sense transcript for its overlapping with exon 5 of WNT4 gene. Interestingly, Coding-Potential Assessment Tool (CPAT) using an alignment-free logistic regression model predicts a coding potential for this lncRNA (in the public domain, https://hg19.lncipedia.org/).⁴⁷ Although no significant change in WNT4 mRNA in KC-affected corneas, this does not rule out the potential role of lnc-WNT4-2:1 in regulating WNT signaling pathway. 

Another pathway in KC pathogenesis is the PI3K/AKT pathway. This pathway is involved in cell cycle regulation and cell proliferation in the cornea.⁷⁸ Our data show downregulation of negative regulators of PI3K/AKT pathway, including PIK3IP1, DDIT4, IGFBP-3, -4, and -5. PIK3IP1 may suppress AKT activity by inhibiting the expression of PI3K.⁷⁹ PIK3IP1 has been reported to be downregulated in KC-affected corneas.²⁷ Insulin-like growth factor biding proteins (IGFBPs) belong to a family of six structurally similar proteins, which antagonize the binding of insulin growth factors to their receptors, thus, controlling cell survival, mitogenesis, and differentiation.⁸⁰ Our results indicate that IGFBP-3, -4, and -5 have 10-, 3-, and 8-fold reduced expression in KC corneas, respectively. Our findings are consistent with previous studies in keratocytes isolated from KC-affected corneas compared to healthy corneas.²⁰ 

For the first time, we have identified pairs of RNAs showing significantly different correlations in KC versus controls (Table 5). Functions of most of the identified genes are still not clear, especially in the cornea. APOD and SFRP1 showed a significant positive correlation in KC, with both being downregulated (Fig. 3B). Both genes were reported to have tumor-suppressive functions, including DNA repair, induction of apoptosis, detoxification, differentiation, and transcriptional regulation.⁸¹ Apolipoprotein D (APOD) is a stroma specific protein as lipophilic molecules carrier.⁸² SFRP1 is a negative regulator of WNT signaling, but it is unclear whether WNT signaling disruption is a cause for the inverse correlation between the two genes (in KC versus controls). 

Despite the significant findings, our study has several limitations. First is the relatively small sample size of cases and controls. Second, age, sex, and ethnicity were not matched completely between cases and controls. It will be appropriate to study how age affects the corneal expression in controls so that we can identify KC-specific but not age-related expression alterations. The source of corneal samples (i.e., surgical versus postmortem donors) might potentially contribute to the expression profile difference between KC and controls. Third, this study used whole cornea for expression. It will be important to examine tissue-specific differential expression in the epithelium and stroma separately if possible. Fourth, our expression study is limited by transcript quantification only. Additional proteomics or protein-based studies are necessary to confirm the differential expression in KC-affected corneas. 

We have not only confirmed many previous reported KC-related differential expressed genes, but also identified a number of novel coding and noncoding RNAs involved in KC pathogenesis. We have identified overabundance of genes with reduced expression, consistent with the destructive nature of corneal degradation associated with KC-related corneal thinning and with our previous RNA microarray study and another study from Macé et al.^28,41 

In summary, we have successfully used Total RNA-Seq to identify differentially expressed coding and noncoding RNAs in KC-affected whole-corneal tissues. These RNAs may disturb various corneal functions including cellular adhesion, oxidative stress response, proliferation, migration, apoptosis, and wound healing. Our expression correlation analysis between differentially expressed coding RNAs and lncRNAs has identified a number of potential expression regulators and provided potential functional annotation for lncRNAs with unknown functions. 

Acknowledgments

The authors thank all the donors for the human corneal tissue samples. Without such samples, this study would not be possible. 

Supported by National Institutes of Health (Bethesda, MD, USA) grants R01EY023242, R01EY023646, R01EY009052, and P30 EY005722. 

Disclosure: M.L. Khaled, None; Y. Bykhovskaya, None; S.E.R. Yablonski, None; H. Li, None; M.D. Drewry, None; I.F. Aboobakar, None; A. Estes, None; X.R. Gao, None; W.D. Stamer, None; H. Xu, None; R.R. Allingham, None; M.A. Hauser, None; Y.S. Rabinowitz, None; Y. Liu, None