Transcriptional profiling of corneal stromal cells derived from patients with keratoconus

Keratoconus (KC) is a multi-factorial corneal ectasia with unknown etiology affecting approximately 1:2000 people worldwide. Dysregulated gene expression, using RNA-Seq technology, have been reported in KC corneal tissue. However, the differential expression of genes, in KC corneal stromal cells have been widely ignored. We utilized mRNA-Seq to analyze gene expression in primary human corneal stromal cells derived from five non-Keratoconus healthy (HCF) and four Keratoconus (HKC) donors. Selected genes were further validated using real time PCR (RT-PCR). We have identified 423 differentially expressed genes with 187 down- and 236 up-regulated in KC-affected corneal stromal cells. Gene ontology analysis using WebGestalt indicates the enrichment of genes involved in cell migration, extracellular matrix, adherens junction, and MAPK signaling. Our protein-protein interaction network analysis identified several network seeds, such as EGFR, NEDD4, SNTA1, LGALS3BP, HSPB1, SDC2, MME, and HIF1A. Our work provides an otherwise unknown information on the transcriptional changes in HKCs, and reveals critical mechanisms of the cellular compartment. It also highlights the importance of human-based in vitro studies on a disease that currently lacks strong biomarkers and animal models.

We have previously shown that human corneal stromal cells isolated from KC patients, termed HKCs, maintain characteristics of the disease phenotype in vitro and secrete and assemble a thinner extracellular matrix (ECM) high in the pro-fibrotic marker, collagen type III, compared to their healthy counterparts 26,27 . Furthermore, consistent with tear analysis of KC patients in vivo [28][29][30] , we have found that HKCs maintain altered cellular metabolism 31 correlating to increased oxidative stress 32 . In this study, for the first time, we have investigated the transcriptome profiles of primary human corneal stromal cells from Healthy and KC donors using RNA-Seq to identify global changes in transcript levels that may contribute to the disease phenotype. Our bioinformatics analysis and RT-PCR validation implicated the potential involvement of genes in several KC-related pathways. If in vitro studies are to lead the way towards future KC therapeutic targets, it is critical to examine the primary corneal cells used. Our study provides, otherwise unknown, information to our understanding of the disease mechanisms.

Methods
Ethics and inclusion criteria. Primary human corneal stromal cells were isolated from five healthy/normal and four KC individuals. The Institutional Review Board at the University of Oklahoma Health Science Center -Dean McGee Eye Institute approved the protocol and studies, prior to initiation. This study adhered to the tenets of the Declaration of Helsinki. Prior to collection of corneas, written permission was obtained from all subjects following corneal transplantation. None of the KC donors was from a vulnerable population and all donors or next of kin provided written informed consent that was freely available. Patients who had any other ocular or systemic disease, or previously underwent collagen crosslinking, were excluded from this study. All donors were examined comprehensively using Pentacam HR, refraction, and slit lamp to confirm KC diagnosis. Stabilization of corneal thinning or KC progression was not reported by the clinician prior to tissue isolation and remains unknown.
Healthy corneal tissues with no history of ocular or systemic diseases was provided by, the National Development and Research Institutes (NDRI).
Cell isolation and expansion. Primary corneal stromal fibroblasts were isolated as previously described 33,34 .
KC corneal tissue was provided from the Dean McGee Eye Institute in Oklahoma City, OK. Corneal stromal cells were isolated, by removing both the epithelial and endothelial layers with a sterile surgical scalpel. Tissues were then cut into small pieces (~2 × 2 × 2 mm), and incubated in sterile flasks to promote adhesion. Explants were supplemented with EMEM containing 10% fetal bovine serum (FBS, Atlanta biologicals, Flowery Branch, GA) and antibiotic/antimycotic (anti/anti, Life Technologies, Grand Island, NY), and incubated at 37 °C/5% CO 2 for 2-3 weeks until cells migrated from the explant (~70-80% confluence).
Both Healthy Corneal Fibroblasts (HCFs) and Human Keratoconus cells (HKCs) were cultured on T75 tissue culture flasks. Cells were seeded at 5,000 cells/cm 2 and maintained at 37 °C with 5% CO 2 in a humidified incubator. Cells used in all experiments were between 2 nd and 4 th passages.
RNA extraction. Total RNA was extracted using the Ambion RNA mini extraction kit (Ambion TRIzol ® Plus RNA Purification Kit: Life technologies, Carlsbad, CA). Briefly, TriReagent ® (Life Technologies Corporation, Carlsbad, CA), was added to the cell layer after aspiration of medium from the culture and brief washing with Phosphate Buffered Saline (PBS). The cellular layer was mechanically disrupted from tissue culture plate using gentle pipetting. Phase separation was conducted with chloroform and the total RNA contained in the aqueous phase was purified using RNeasy ® mini kit column (QIAGEN, Hilden, Germany), according to the manufacturer's protocol. Three extractions were carried out for each sample and pooled at the end of the RNeasy protocol. The purity and quantity of total RNA were evaluated using an ultraviolet spectrometer (CLARIOstar, BMG LABTECH, Cary, NC).

RNA-Seq and bioinformatics analysis.
RNA sequencing library was prepared using TruSeq RNA Library Prep kit, from Illumina with sample-specific indexes at the Oklahoma Medical Research Foundation (OMRF) Genomics facility. Briefly, 1 µg of total RNA was used to purify and fragment mRNA, followed by first and second strand cDNA synthesis. After purification and ligation of Illumina adapters and sample-specific indexes, sequencing libraries were validated using Agilent DNA 100 kit with Agilent Bio analyzer 2100. The libraries were normalized, and pooled for sequencing, with Illumina NextSeq 550 system using the paired end 75 bp with the high-output kit.
The data analysis for RNA-Seq was done using our established pipeline as described previously 13 . Briefly, after quality check and quality control with all the sequencing reads, demultiplexed reads were aligned by TopHat in paired-end reading with the approximation of the median library size 13 . Counts of sequencing reads were normalized using Cufflinks in fragments per kilo bases and millions reads (FPKM) 13 . Normalized sequencing read counts were analyzed by Cuffdiff 13 , with a transcript file from Ensembl database for the annotations at the gene as well as the isoform levels in a group comparison manner such as 4 cases vs 5 controls in this experimental design.
Missing expression data in specific samples was replaced with a value of 0.001 to enable the differential analysis between cases and controls. Without such replacement of missing data, the related mRNAs would have been excluded from the differential expression analysis. In the final gene list with the positions at the corresponding chromosomes, it showed comparative fold changes and adjusted P values (false discovery rate, FDR) through the resulting analysis of four cases vs. five controls. Differentially expressed genes were defined to have an FDR value ≤ 0.05 and a |fold change| ≥ 2.
The differentially expressed genes were loaded into WEB-based Gene Set Analysis Toolkit (WebGestalt) 35,36 , for bioinformatics analysis to identify enriched gene ontologies, pathways, network modules, and associations with phenotypes/disease/drug. Potential targeted miRNAs were identified for the differentially expressed genes using WebGestalt.
Validation using Real-time PCR. To validate the differential expression of selected genes, we used RT-PCR with all the donors/samples, as previously described 32,37,38 . The cDNA synthesis was followed using a SuperScript ™ III First-Strand Synthesis SuperMix (Invitrogen, Carlsbad, CA) according to the manufacturer's protocol. The TaqMan gene expression assays (Applied Biosystems, Foster City) for GAPDH (Hs99999905_m1) and 18S (Hs99999901_s1) were used as the reference assays to normalize target gene expression. We selected the following genes for validation: KRT7 (Hs00559840_m1), MME (Hs00153510_m1), ANKRD1 (Hs00173317_m1), ERG1 (Hs00152928_m1), IL1B (Hs01555410_m1) CXCL1 (Hs00236937_m1), and GDF15 (Hs00171132_m1). Furthermore, 10 ng of cDNA was used for initiating the PCR reaction for a 20-µl reaction mixture containing our desired probes and the TaqMan Fast Advanced MasterMix (Applied Biosystems, Life technologies, Foster City, CA). Amplification of samples was performed using the StepOnePlus TM real-time PCR system (Life Technologies) in accordance with the manufacture's protocol. Graph Pad Prism 7 and MS-Excel were used for data analysis.
Statistical analysis. GraphPad Prism 7.02 was used to determine statistical significance using ANOVA or t-test, where appropriate. A p ≤ 0.05 was considered statistically significant.

Results
Clinical phenotypes. Our study included HCF cells from 5 postmortem donors (3 males/2 females) with postmortem delay less than 24 hours. All four HKC cells were derived from surgically removed keratoconic corneas immediately after surgery. It included 1 male and 3 females. The average age was 47.2 years for controls and 60.8 years for cases.

HCFs and HKCs -transcriptional profiles.
We performed RNA-Seq with primary human corneal stromal fibroblast cells derived from five unaffected controls (HCFN21, HCFN20, HCFN23, HCFN22, and HCFN4) and four KC patients (HKCWV1, HKCWV2, HKCDM1, and HKCDM2). Individual information for all the samples is provided in Table 1.
All nine RNA samples were sequenced with 33.7-47.9 million paired-end reads of 75 nucleotides. All the samples had at least 30 million paired sequence reads aligned. Sequence reads with multiple alignments were removed during quality control process. Out of the total 15,159 genes, 11,540, 11,637, 11,568, 11,577, and 11,740 genes were expressed in the individual HCF cells with FPKM ≥ 1.0 and 9,870 genes were expressed in all five primary HCF cells with FPKM ≥ 1.0. 11,791, 11,768, 11,751, 12,122, and 11,544 genes were expressed in the individual HKC cells with FPKM ≥ 1.0 and 9,923 genes were expressed in all five primary HKC cells with FPKM ≥ 1.0. With all corneal stromal fibroblast cells from 10 donors, 9,302 genes were expressed with FPKM ≥ 1.0.
Our differential analysis using Cuffdiff identified 423 differentially expressed genes in HKC cells with at least 2 fold change and false discovery rate q value ≤ 0.05 (supplemental Table 1). There were 187 down-regulated and 236 up-regulated genes in KC-affected corneal stromal cells. After applying a minimum expression of 10 normalized reads in either cases or controls, a total of 208 genes (121 down and 87 up-regulated) were differentially expressed in HKC fibroblast cells (supplemental Table 2). The top 20 up-and down-regulated genes with largest fold changes are listed in Table 2. www.nature.com/scientificreports www.nature.com/scientificreports/ Diverse and dominant pathways in HKCs. Gene ontology analysis using all 423 genes indicated the significant enrichment of genes coding for proteins involved in or related with cell migration, collagen-containing ECM, adherens junction, intrinsic to plasma membrane, cytokine and growth factor activity, growth factor binding, and kinase activity. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis identified the enriched pathways including MAPK signaling and Rap1 signaling. Reactome analysis indicated the involvement of PTK6 promoting HIF1A stabilization, non-integrin membrane-ECM interactions, ECM organization, and syndecan interactions. Wikipathway analysis suggested the enrichment of genes in MAPK signaling, differentiation, and focal adhesion-PI3K-Akt-mTOR-signaling pathway. Seven genes (DAPK1, ELK4, HIF1A, MEF2C, RPS6KA3, and SGK1) are targets of MAPK7 kinase activity. Many of the differentially expressed genes are enriched targets of miR-26a/b, miR-203, miR-380-3p, and miR-96 with 28, 23, 12, and 22 target genes respectively. Our miRNA profiling in seven normal human corneal tissues 39 indicated the expression of these miRNAs, suggesting their potential role in corneal function. Protein network analysis identified a subnetwork with several connection hubs or seeds, including EGFR, NEDD4, SNTA1, LGALS3BP, HSPB1, SDC2, MME, HIF1A, CBL, and ERRFI1 (Fig. 1).  www.nature.com/scientificreports www.nature.com/scientificreports/ With the relatively highly expressed 208 differential genes, gene ontology analysis indicated the enrichment of genes involved with cell migration/motility, focal adhesion, endoplasmic reticulum lumen, and ECM, similar to the bioinformatics findings with all 423 genes.

RT-PCR validation.
We further validated the differential gene expression by selecting seven genes from the top 20 genes identified by transcriptomics. We focused on novel pathways associated with proliferation, wound healing, and pro-inflammatory pathways. We compared our RNA-Seq results with RNA expression using RT-PCR. Of the genes tested, ANKRD1 and KRT7 were significantly down-regulated in HKCs when compared to controls (Fig. 2), consistent with our RNA-Seq findings. In contrast, EGR1, IL1B, GDF15, MME, and CXCL1 RNA expression was not significantly different between HKCs and Controls, where RNA-Seq analysis showed significant up-regulation in HKCs (Fig. 2).

Discussion
Previous reports applying transcript and gene analysis approaches in KC have identified altered expression ranging from a diverse group of pathways, including expression of keratocan, ECM proteins, Wnt signaling, lysyl oxidase, and inflammatory genes [12][13][14][15][16][17][18][19][20][21][22][23][24] . These studies have highlighted the heterogeneity of the KC condition and the need to identify common biomarkers of KC. In this study, we investigated whether gene expression in corneal stromal cells is altered by KC. For the first time, RNA-Seq was used to compare the global gene expression of the www.nature.com/scientificreports www.nature.com/scientificreports/ corneal stroma derived cells from KC and healthy controls. Bioinformatics were used to delineate critical pathways and highly modulated genes were validated by RT-PCR. Our analysis indicated the potential involvement of several pathways in KC pathogenesis via stromal cells, including cell migration, cellular adherence, and ECM.
The transcriptional factor, ANKRD1, is known to function in wound healing responses in cardiac tissue 40 via regulation of the ERK and TGF-β pathways, among others 41,42 , though its expression and role in the cornea has yet to be reported. Furthermore, knockout of ANKRD1 in the mouse has shown delayed wound healing in skin tissue, due to defects in interactions between dermal fibroblasts and the surrounding ECM 43 . The interplay of ANKRD1 expression and the TGF-β pathway has also been reported 44,45 , further establishing a functional role in wound healing 46 . Interestingly, in our study, we have identified a downregulation of ANKRD1 in HKCs under homeostatic conditions validated in both the RNA-Seq and RT-PCR-data, which supports the hypothesis of an abnormal wound healing phenotype in KC 47,48 . TGF-β is known to also regulate Keratin7 49 . Keratin 7 is a protein coding gene, which plays a role in DNA synthesis. Krenzer and Freddo reported Keratin 7 expression in the conjunctiva 50 . Elder et al. described Keratin 7 in the basal and suprabasal epithelial cells of the central cornea 51 . However, it was not found in the central corneal epithelium 52 . While Keratin 7 is not normally reported in the corneal stromal layer or in KC studies, our data shows downregulation of Keratin 7 in HKCs. This finding could indicate a role in the context of wound healing and loss of corneal transparency. Keratin 7 isoform polymerizes to form the intermediate cytoskeletal filaments that provide structure and stability to corneal epithelial cells. Thus, it is associated with cytoskeletal signaling and remodeling pathways 49 . Through many factors and specifically TGF-β regulation, Keratin 7 has also been implicated in cellular stress responses to wound healing and tissue repair 53,54 . Acting through the recruitment and activation of SMADs to regulate gene expression. We, and others, have reported the role of SMADs in KC [55][56][57][58][59] . Given that KC has previously been associated with an altered response to TGF-β1 and TGF-β3 ligands 57,60 , the downregulation of both ANKRD1 and Keratin 7 in HKCs may contribute to this complex interplay favoring a pro-fibrotic phenotype in HKCs compared to healthy controls. The functional relevance of our findings to the KC condition may be related to these downstream pathways regulated by TGF-β signaling, a key pathway important in wound healing and ECM deposition in corneal biology 61 . Further studies to determine the functional effects of decreased gene expression in the context of KC development or progression are warranted. Moreover, validation of protein level changes are needed to determine that altered gene expression identified in our study contributes to protein level changes.
Our study also has a few limitations. Firstly, the sample size is relatively small. Large sample size of both cases and controls with similar age range will definitely increase the statistical power to detect additional signaling pathways and networks. Secondly, we have used primary corneal stromal fibroblast cells instead of the isolated stromal layer tissue. The culture process may influence the expression profile. However, studies from our lab 27,57 Figure 2. Real-time PCR validation of selected genes for their differential expression from RNA-Seq. The expression of (a) ANKRD1, and (b) Keratin7 was successfully validated using RT-PCR with significant differential expression with p values < 0.0001, and 0.0365 respectively. The expression of (c) ERG1, (d) GDF15, (e) MME1, (f) IL1B, and (g) CXCL1 was not validated for their differential expression from RNA-Seq (p value > 0.05). (2019) 9:12567 | https://doi.org/10.1038/s41598-019-48983-8 www.nature.com/scientificreports www.nature.com/scientificreports/ and others 62,63 have consistently supported maintenance of an altered disease phenotype by KC-derived stromal cells cultured in vitro with characteristics parallel to the in vivo condition, including altered ECM expression and regulation 9 . Thirdly, the identified KC-associated signaling pathways and networks could either be the disease outcome or may actively contribute to KC pathogenesis, which is difficult to discern in a snapshot of transcript levels. A systematic approach to determine the role of specific pathways in KC development from onset to further progression to the late-stage disease is required to validate disease causation. However, this is difficult with the lack of tissue/cell isolation from early-stage KC patients. Lastly, the age differences between control and KC patients was ~13.6 years with an average age of 47.2 ± 18.3 years for controls and 60.8 ± 11.6 years for KC. The effects of age on gene expression has previously been reported with downregulation associated with collagen 64 and elastin expression by fibroblasts of the skin 65 . Likewise, aging influences corneal structure with increased keratocyte senescence and altered ECM rigidity 66 . Thus, the age gap in our study between groups may have contributed to differential gene expression independent of the KC disease.
In summary, using RNA-Seq technology, we have successfully profiled the differential gene expression in primary human corneal stromal cells derived from patients with KC. This data will advance our molecular understanding of the pathogenesis of KC.

Data Availability
The datasets generated during the current study are available from the corresponding author on reasonable request.