Altered Regulation of mRNA and miRNA Expression in Epithelial and Stromal Tissue of Keratoconus Corneas.

Purpose: Evaluation of mRNA and microRNA (miRNA) expression in epithelium and stroma of patients with keratoconus.
Methods: The epithelium and stroma of eight corneas of eight patients with keratoconus and eight corneas of eight non-keratoconus healthy controls were studied separately. RNA was extracted, and mRNA and miRNA analyses were performed using microarrays. Differentially expressed mRNAs and miRNAs in epithelial and stromal keratoconus samples compared to healthy controls were identified. Selected genes and miRNAs were further validated using RT-qPCR.
Results: We discovered 170 epithelial and 1498 stromal deregulated protein-coding mRNAs in KC samples. In addition, in epithelial samples 180 miRNAs and in stromal samples 379 miRNAs were significantly deregulated more than twofold compared to controls. Pathway analysis revealed enrichment of metabolic and axon guidance pathways for epithelial cells and enrichment of metabolic, mitogen-activated protein kinase (MAPK), and focal adhesion pathways for stromal cells. Conclusions: This study demonstrates significant differences in the expression and regulation of mRNAs and miRNAs in the epithelium and stroma of Patients with KC. Also, in addition to the well-known target candidates, we were able to identify further genes and miRNAs that may be associated with keratoconus. Signaling pathways influencing metabolic changes and cell contacts are affected in epithelial and stromal cells of patients with keratoconus.

Keratoconus (KC) is a disease of the cornea associated with thinning and deformation. The disease usually begins during puberty with an incidence of about 1 in 400 individuals. The disease affects the stroma, as well as the epithelial layer of the cornea; however, the causes of KC are still largely unexplained. Several studies have been conducted to better understand the pathogenesis of KC. For these studies, among others, tear fluids were examined indicating altered inflammatory cytokines and proteinases. Cell culture experiments and studies from aqueous humor of patients with KC have also shown altered metabolic pathways, changes in hormonal components, and inflammatory cytokines.– In genome-wide association studies and genome-wide linkage studies, numerous candidate genes have already been evaluated. A familial association in addition to environmental factors or mechanical factors such as eye rubbing has also been discussed., In order to obtain information about functional genes and pathways, RNA subanalyses using mRNA arrays or RNA sequencing have been performed.– When comparing the regulated protein-coding RNAs of these studies, several KC-related genes were identified, including SFRP1, AQP5, FBLN1, KRT16, and S100A9. These studies have been performed on corneal epithelium, on total corneal tissue, or on cultivated corneal fibroblasts, which means that the epithelial cells and stromal cells (i.e., the keratocytes) were examined in the same sample.– However, there appear to be major expression differences between epithelium and stroma, as shown in the study by Yam et al. regarding protein profiles of KC corneas. In both cell types, different metabolic pathways are affected. The epithelial proteome showed more changes in cell metabolism and mitochondrial involvement, whereas the stromal proteome showed more changes in cellular assembly and tissue organization. To the best of our knowledge, no studies have been performed on the epithelium and stroma of patients with KC separately for microRNA (miRNA) and mRNA expression using microarray. In contrast to coding mRNAs, the miRNAs are short non-coding RNAs of approximately 22 nucleotides. The major role of miRNAs is to regulate protein translation of mRNAs on a post-transcriptional level. Mature miRNAs are incorporated into the RNA-induced silencing complex (RISC) and bind to their target mRNAs in a sequence specific manner, which leads to either inhibition of protein translation or degradation of target mRNAs. Numerous miRNA–target gene interactions have already been identified. Through their regulation of multiple target genes, miRNAs are able to affect a multitude of cellular pathways, including proliferation, migration, and angiogenesis in carcinomas. It is conceivable that deregulated miRNAs could also have a significant impact on the disease process. The role of miRNAs and their functionality have not been thoroughly explored in diseases of the eye and KC in particular. Therefore, it is of great importance not only to study the deregulations at mRNA and protein level but also to analyze the influence of miRNAs. As of yet, to the best of our knowledge, there has been no systematic study on the expression of miRNAs in epithelium and stroma of patients with KC.

To systematize the pathogenesis of KC, it could be important to study gene and miRNA expression of epithelial and stromal cells separately, to explore different starting points and possible regulatory components here as well. The purpose of this study was to investigate the differences in mRNA and miRNA expression separately in the epithelium and stroma of patients with KC to gain deeper insight into the cellular mechanisms deregulated during KC development.

Materials and Methods

All subjects had KC confirmed by clinical records and by tomography using a Pentacam Oculus (Arlington, WA, USA). Demographic characteristics and KC grading of the included corneas according to Belin and Duncan are given in Supplementary Table S1.

Patient and Donor Samples

Eight corneas of eight patients (50% male; mean age, 49.1 ± 13.3 years) with KC from elective penetrating keratoplasties were recruited from the Department of Ophthalmology, Saarland University Medical Centre, Homburg (Saar), Germany. Prior to corneal excision, the corneal epithelium was removed so that the epithelium and the stroma of each KC sample could be analyzed separately.

As controls for the epithelium, eight samples (eight eyes of eight patients; 50% male; mean age, 58.8 ± 15.0 years) with epithelial basal membrane dystrophy (EBMD) and one sample from recurrent erosion syndrome were used. The control stromal tissues for mRNA expression profiling were obtained from the LIONS Cornea Bank Saar-Lor-Lux (Trier/Westpfalz, Germany). Donor tissues were used for Descemet membrane endothelial keratoplasty (DMEK; 63% male; mean age, 75.0 ± 11.2 years).

The stromal samples for miRNA microarray profiling were obtained from patients that underwent small incision lenticule extraction (SMILE) at the Eye Department of the Center for Refractive Surgery (St. Francis Hospital, Münster, Germany) and at the Visual Center for Eye Laser and Eye Surgery (Graz, Austria). These samples were pooled by age and sex to obtain the same amount of RNA as from the KC corneas. In total, we used the lenticules of 48 patients for this study (50% male; mean age, 32.1 ± 6.9 years). Detailed demographic data of all patients and donors are displayed in Supplementary Table S1. Written informed consent was obtained from all subjects. This research conformed to all tenets of the Declaration of Helsinki. The Ethics Committees of the Medical Association of Saarland (No. 74/19), of Westfalen-Lippe (No. 2019- 548-b-S) and of Graz (No. 31- 473 ex 18/19) approved the study.

Tissue Preparation

The epithelial samples of patients with KC and patients with EBMD that were used as normal controls were centrifuged, and the cell pellet was taken up in lysis buffer for storage in −80°C. After penetrating keratoplasty of patients with KC or use of donors for DMEK surgery, the stromal samples were cut into several small pieces of approximately 3 mm3, were directly taken up in lysis buffer (Buffer RL of the Total RNA Purification Plus Micro Kit; Norgen Biotek, Thorold, Ontario, Canada), and were frozen at −80°C until RNA isolation. The lenticules of the patients who underwent SMILE were snap frozen at −18°C and were shipped with dry ice for RNA extraction.

RNA Isolation and Quality Control

No additional preparation of epithelial samples was required for the isolation of RNA. To obtain RNA from KC, normal, and lenticule stroma samples, the samples were first homogenized using a tissue homogenizer (Precellys Evolution Super Tissue Homogenizer; Bertin Technologies, Saint Quentin, France), for four cycles, each at 4500 rpm for 30 seconds, followed by a 2-minute cooling time. Stromal samples were than incubated with 0.4 mg/mL proteinase K (Norgen Biotek) for 15 minutes at 55°C. The RNA isolation was performed according to the manufacturer's protocol, and the RNA was stored at −80°C. The RNA including the miRNA was eluted in 30 µL RNase-free water. The RNA concentration was measured using a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA), followed by measurement of the RNA integrity using an Agilent Bioanalyzer and RNA Nano Kit (Agilent Technologies, Santa Clara, CA, USA).

mRNA Expression Profile

Gene expression profiles of epithelial and stromal samples of KC and normal samples were measured using SurePrint G3 Human Gene expression v3 8×60K microarrays and the Low Input Quick Amp Labeling Kit (Agilent Technologies). Total RNA (100 ng) of each sample was used for complementary DNA (cDNA) synthesis using a T7 primer. The cDNA was transcribed into cRNA using a Cy3-labeled pCp and T7 polymerase master mix. Following determination of RNA concentration, samples were hybridized to microarrays and scanned using an Agilent microarray scanner with a resolution of 3 µm, followed by analysis using Agilent Feature Extraction software (version 10.10.1.1).

miRNA Microarray Profiling

The Agilent SurePrint G3 miRNA microarray and the Complete Labeling and Hybridization Kit (Agilent Technologies) were used to measure the expression of 2549 human mature miRNAs of miRBase 21 (http://www.mirbase.org). The RNA was labeled using Cy3-pCp following dephosphorylation. Following hybridization for 20 hours at 55°C and 20 rpm, microarrays were scanned using an Agilent G2505C DNA Microarray Scanner.

Statistical Analysis

MiRNA and gene expression data were analyzed using Agilent GeneSpring software. Data was quantile normalized, log transformed, and filtered for protein-coding genes and miRNAs expressed in at least 25% of samples, leaving 734 miRNAs and 16,800 mRNAs for analysis. Differentially expressed miRNAs and mRNAs in KC samples were identified using unpaired two-sided t-tests. Benjamini–Hochberg false discovery rate adjustment was used to correct P values for multiple testing. Genes and miRNAs with adjusted P < 0.05 were considered significantly altered. Microarray data for mRNA and miRNA are available at Gene Expression Omnibus (GEO) database accession numbers GSE204791 and GSE204838, respectively. Heat maps to demonstrate biological variances between samples were created using Morpheus (https://software.broadinstitute.org/morpheus).

Validation of Deregulated mRNAs and miRNAs Using RT-qPCR

From the top 20 down- and upregulated protein-coding genes of the epithelial and stromal samples, seven candidates were selected for validation by RT-qPCR. For miRNA validation, eight were randomly selected from the top 40 dysregulated miRNAs. The selected transcripts are listed in Supplementary Table S1.

RT-qPCR Validation of Deregulated miRNAs and Transcripts

To validate the mRNA and miRNA expression data obtained from the microarray analyses, we performed RT-qPCR for selected candidates in the same samples that were used in the microarray experiment. The cDNA was generated from 500 ng RNA (mRNAs) and 50 ng RNA (miRNAs) using the QuantiTect RT Kit and the miScript II RT Kit (both from QIAGEN, Hilden, Germany), respectively, according to the manufacturer's recommendations. qPCR was performed in duplicate using the QuantiTect SYBR Green PCR Kit (QIAGEN) with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and RNU6B as respective endogenous controls for mRNA and miRNA analysis. Primers used for mRNA and miRNA are displayed in Supplementary Table S1. RT-qPCR was run on a StepOne Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The expression of selected mRNAs or miRNAs was normalized to their respective endogenous controls, and deregulation of expression between KC and normal tissues was calculated using the ΔΔCq method (PMID: 11846609).

Pathway Analysis

To identify cellular pathways enriched for deregulated genes, the GeneTrail platform (https://genetrail2.bioinf.uni-sb.de/) was used to perform an overrepresentation in biological processes and pathways (PMID: 26787660) utilizing a list of deregulated mRNAs with more than twofold expression change and P < 0.05 as input list.

Results

Sample Characteristics

A detailed overview of the samples is shown in Table 1. The epithelium and stroma of corneas from eight patients with an advanced stage of KC who underwent keratoplasty were used for microarray analyses. Because it was not possible to obtain epithelium and stroma from the same subjects as control tissue, the epithelium of patients with EBMD or recurrent erosion syndrome and the stroma of donor corneas (previously used for DMEK) and stroma of subjects who underwent SMILE were used as controls.

Table 1.

Demographic Data for Subjects Used in This Study

Demographic Data for KC Samples Used for Epithelial and Stromal mRNA and miRNA Profiling
Epithelium	Stroma	Sex	Age (y)	ABCD Grading
1 E KC	1 S KC	F	59	A4 B4 C3 D4 +
2 E KC	2 S KC	M	28	A4 B4 C3 D1 −
3 E KC	3 S KC	F	58	A4 B4 C3 D4 +
4 E KC	4 S KC	M	56	A4 B4 C3 D4 ++
5 E KC	5 S KC	F	50	A2 B4 C1 D3 −
6 E KC	6 S KC	M	67	A2 B4 C1 D1 +
7 E KC	7 S KC	F	29	A4 B4 C3 D4 +
8 E KC	8 S KC	M	46	A4 B4 C2 D3 −
		M, 50%	Mean, 49.1 ± 13.3
Demographic Data for Subjects Used for Epithelial mRNA and miRNA Profiling
Epithelium		Sex	Age (y)	Disease
1 E N		F	56	EBMD
2 E N		M	65	EBMD
3 E N		F	58	EBMD
4 E N		M	81	EBMD
5 E N		F	50	EBMD
6 E N		M	30	Erosio
7 E N		F	53	EBMD
8 E N		M	77	EBMD
		M, 50%	Mean, 58.8 ± 15.0
Demographic Data for Healthy Donors Used for Stromal mRNA Profiling (Corneas Used for DMEK)
Stroma		Sex	Age (y)
1 S N		M	87
2 S N		F	87
3 S N		M	87
4 S N		M	70
5 S N		M	79
6 S N		F	70
7 S N		F	55
8 S N		M	65
		M, 63%	Mean, 75.0 ± 11.2
Demographic Data for Subjects Used for Stromal miRNA Profiling (SMILE)
Stroma		Sex	Mean Age (y)
1 S N L		F	25
2 S N L		M	27
3 S N L		F	29
4 S N L		M	29
5 S N L		F	35
6 S N L		M	26
7 S N L		F	42
8 S N L		M	44
		M, 50%	Mean, 32.1 ± 6.9

Each sample was pooled from two lenticules of six patients (age- and sex-matched). F, female; M, male; Erosio, recurrent erosion syndrome.

Expression Profiles

The top 20 down- and upregulated mRNAs and miRNAs are displayed in Tables 2 and 3, respectively.

Table 2.

The 20 Most Down- and Upregulated mRNAs in KC Epithelial Samples and KC Stromal Samples, Ordered by Increasing FC

Deregulated mRNAs in KC Epithelium			Deregulated mRNAs in KC Stroma
Gene Symbol	FC	Adjusted P	Gene Symbol	FC	Adjusted P
PDE8B	−8.777	0.0300	DHRS9	−311.730	0.0018
AKR1B10	−8.560	0.0295	S100A9	−276.773	0.0006
CCDC3	−8.554	0.0177	S100A8	−260.582	0.0009
MGST1	−8.257	0.0493	SAA2	−252.721	0.0019
SEMA3C	−7.066	0.0136	PI3	−231.211	0.0002
NELL2	−6.929	0.0325	SPRR3	−209.532	0.0043
AKR1B15	−6.509	0.0192	KRT16	−127.746	0.0053
HPGD	−5.380	0.0439	KRT6C	−119.859	0.0103
AKT3	−5.267	0.0407	S100P	−106.877	0.0073
GALNT14	−5.158	0.0236	KRT6A	−98.903	0.0204
FMO1	−5.142	0.0321	SPRR1B	−98.501	0.0074
ALDH2	−4.940	0.0463	SPRR2A	−95.719	0.0152
SYNPR	−4.910	0.0355	LCN2	−89.694	0.0029
FST	−4.754	0.0443	FMO1	−89.386	0.0005
CRYM	−4.383	0.0463	SPRR2E	−83.699	0.0114
FAM110B	−4.335	0.0377	CXCL1	−81.210	0.0038
KHDRBS3	−4.210	0.0284	SPRR2F	−76.949	0.0113
TMEM64	−4.154	0.0306	KRT13	−67.664	0.0059
DHRS9	−4.098	0.0323	SPRR2D	−67.215	0.0113
ARNTL	−4.097	0.0148	DEFB4A	−65.555	0.0067
NANOG	2.697	0.0288	CNTNAP2	9.847	0.0280
ARL17B	2.730	0.0245	CD36	9.940	0.0435
ZBTB1	2.758	0.0233	DAPL1	10.221	0.0283
GPR135	2.782	0.0190	CD86	10.424	0.0272
AVIL	2.888	0.0333	LRRC2	10.435	0.0279
CDRT1	2.891	0.0284	PVALB	11.676	0.0006
SERPINB13	2.948	0.0454	LEFTY1	11.849	0.0332
SMG1	2.977	0.0319	SPARCL1	11.899	0.0039
GPR75	3.014	0.0199	GPR179	11.941	0.0235
SEMA6C	3.162	0.0494	TLDC2	13.112	0.0161
NYNRIN	3.228	0.0132	SLC16A12	13.210	0.0009
SUSD2	3.501	0.0418	RASGRF1	13.256	0.0465
PLEKHN1	3.692	0.0204	ANGPTL7	14.293	0.0080
TNNT2	3.980	0.0203	PCK1	14.854	0.0056
CHST2	4.063	0.0170	FILIP1L	14.892	0.0029
SLC7A4	4.120	0.0101	HLF	15.132	0.0004
METTL21A	4.301	0.0056	CA6	26.418	0.0232
ABO	4.563	0.0264	CYP26A1	54.922	0.0045
ANKRD36B	5.678	0.0133	MTUS2	85.062	0.0013
SFRP1	7.159	0.0063	RPRML	123.205	0.0002

Table 3.

The 20 Most Down- and Upregulated miRNAs in KC Epithelial Samples and KC Stromal Samples, Ordered by Increasing FC

Deregulated miRNAs in KC Epithelium			Deregulated miRNAs in KC Stroma
miRNA	FC	Adjusted P	miRNA	FC	Adjusted P
miR-634	−11.404	0.0168	miR-936	−37.140	3.36E-08
miR-1229-3p	−10.618	0.0169	miR-208a-5p	−31.275	2.05E-09
miR-98-3p	−10.103	0.0340	miR-519e-5p	−24.006	1.33E-06
miR-4664-3p	−10.102	0.0188	miR-3945	−17.742	2.02E-04
miR-6757-3p	−8.737	0.0365	miR-3610	−16.919	1.33E-07
miR-3180-5p	−8.555	0.0266	miR-198	−14.830	6.14E-06
miR-6507-3p	−8.176	0.0365	miR-6739-5p	−14.193	2.18E-05
miR-3190-5p	−7.963	0.0349	miR-125b-1-3p	−13.913	8.75E-04
miR-23c	−6.229	0.0344	miR-3131	−13.378	2.16E-08
miR-2116-3p	−5.082	0.0277	miR-4476	−13.167	3.07E-04
miR-4652-3p	−4.386	0.0365	miR-3180-3p	−12.672	5.78E-10
miR-139-3p	−4.293	0.0253	miR-548q	−12.094	1.24E-05
miR-423-5p	−4.119	0.0136	miR-1247-3p	−11.765	4.40E-07
miR-1908-3p	−4.056	0.0311	miR-1343-5p	−11.721	3.59E-05
miR-204-3p	−3.993	0.0215	miR-5699-5p	−11.321	1.40E-04
miR-4310	−3.868	0.0212	miR-4299	−10.882	5.10E-06
miR-1539	−3.783	0.0110	miR-601	−10.800	4.29E-06
let-7b-3p	−3.531	0.0252	miR-4651	−10.668	1.17E-06
miR-6752-3p	−3.528	0.0215	miR-659-3p	−10.573	9.98E-09
miR-874-3p	−3.465	0.0134	miR-3158-5p	−10.337	1.96E-05
miR-584-5p	11.289	0.0253	miR-210-3p	8.787	2.16E-08
miR-7108-5p	11.371	0.0311	miR-7975	8.966	1.82E-06
miR-371b-5p	11.419	0.0365	miR-34b-5p	8.979	3.64E-05
miR-5195-3p	11.693	0.0253	miR-194-5p	9.487	5.46E-08
miR-3141	13.115	0.0188	miR-197-3p	9.764	1.25E-06
miR-3911	14.852	0.0136	miR-204-5p	10.051	2.77E-03
miR-374a-5p	15.493	0.0310	miR-1260b	10.139	8.17E-06
miR-134-5p	16.897	0.0212	miR-1260a	12.552	5.88E-06
miR-4695-5p	18.809	0.0169	miR-183-5p	12.615	5.55E-03
miR-6757-5p	19.445	0.0056	miR-126-3p	12.874	1.08E-06
miR-6768-5p	19.653	0.0040	miR-149-5p	13.025	2.79E-06
miR-1972	21.387	0.0017	miR-4286	13.761	4.07E-06
miR-6786-5p	22.730	0.0120	miR-96-5p	15.481	2.57E-03
miR-1233-5p	22.789	0.0096	miR-429	17.828	3.56E-03
miR-6778-5p	25.555	0.0029	miR-455-3p	18.912	1.09E-05
miR-6076	26.707	0.0056	miR-205-5p	19.581	5.08E-04
miR-6723-5p	27.429	0.0065	miR-7977	19.822	3.77E-06
miR-6829-5p	29.869	0.0034	miR-15a-5p	22.154	4.29E-06
miR-7114-5p	52.053	0.0018	miR-141-3p	33.241	1.91E-04
miR-211-3p	59.684	0.0001	miR-135b-5p	33.961	9.55E-05

mRNA

A total of 16,800 protein-coding transcripts were expressed in at least 25% of the epithelial and stromal samples. After setting the cutoff (fold change [FC] ≥ 2 between KC and normal samples; t-test adjusted P ≤ 0.05), we identified 170 differentially expressed transcripts in the epithelial samples, from which 82 transcripts were downregulated and 88 transcripts were upregulated in KC samples.

In stromal samples, 1498 differentially expressed transcripts from KC in comparison to healthy controls were identified, of which 915 transcripts were downregulated and 583 transcripts were upregulated in KC samples. The most downregulated mRNAs in KC epithelium or KC stroma were PDE8B, with a FC of –8.7, and DHRS9, with a FC of –311.7, respectively. The most upregulated RNAs in epithelium or stroma were SFRP1 (FC = 7.2) and RPRML (FC = 123.9), respectively. A complete list of deregulated mRNAs in epithelium and stroma is shown in Supplementary Table S5.

The overlap between deregulated transcripts in epithelium and stroma is displayed in a Venn diagram in Figure 1. The Venn diagram clearly shows only limited overlap of deregulated genes in the epithelium and stroma (maximum 1%), whereas the percentage of overlapping deregulated miRNAs is much higher. Notably, there is a 5.8% overlap of miRNAs that are upregulated in the stroma and downregulated in the epithelium. In our analysis, we found several genes known to play a role in KC from previous expression studies, including AQP5, S100A8, EGLN3, EGFR1, and SFRP1.,

Figure 1.

Overlap of deregulated mRNAs and miRNAs in epithelial and stromal cells of KC corneas. The number and identity of the overlapping genes (A) and miRNAs (B) more than twofold deregulated with adjusted P < 0.05 are displayed. A complete list of deregulated genes and miRNAs can be found in Supplementary Table S3. Compared to mRNAs, a higher percentage of overlap is observed for miRNAs.

Through the use of GeneTrail 3.0, eight significantly enriched signaling pathways were identified in epithelial cells, and 27 significantly enriched signaling pathways were identified in stromal cells of patients with KC. Among them, there was an enrichment of deregulated genes in the metabolic signaling pathways for both cell types, as well as an enrichment in the axon guidance signaling pathway in epithelial cells and the focal adhesion signaling pathway in stromal cells of patients with KC (Tables 4 and 5).

Table 4.

Significantly Enriched Pathways and Associated Differentially Regulated Genes in Corneal Epithelial Cells of Patients With KC, Relative to Healthy Corneal Epithelial Cells

Name	Expected Number of Genes	Observed Number of Genes	Adjusted P	Upregulated Genes	Downregulated Genes
Metabolic pathways	10.21	20	0.045	ABO, ACP5, LDHAL6A, PIKFYVE, UGT1A6, UGT1A8	AKR1B10, AKR1C3, ALDH2, AZIN2, DHRS9, GALNT14, GMPR, HGD, HSD17B2, MGST1, NDST1, P4HA2, PDE8B, RIMKLB
Axon guidance	1.32	8	0.007	ABLIM2, NCK2, PLXNA4, RAC3, SEMA3C, SEMA4F	DPYSL2, SEMA6C
Regulation of actin cytoskeleton	1.50	7	0.024	BAIAP2, BDKRB2, PIKFYVE, PPP1R12B, RAC3	ITGB8, MRAS
Steroid hormone biosynthesis	0.4	4	0.024	UGT1A6, UGT1A8	AKR1C3, HSD17B2
Drug metabolism	0.47	4	0.028	UGT1A6, UGT1A8	FMO1, MGST1
Ascorbate and aldarate metabolism	0.18	3	0.024	UGT1A6, UGT1A8	ALDH2
Pentose and glucoronate interconversions	0.22	3	0.028	UGT1A6, UGT1A8	AKR1B10

Table 5.

Significantly Enriched Pathways and Associated Differentially Regulated Genes in Corneal Stromal Cells of Patients With KC, Relative to Healthy Corneal Stromal Cells

Name	Expected Number of Genes	Observed Number of Genes	Adjusted P	Upregulated Genes in Pathway	Downregulated Genes in Pathway
Metabolic pathways	90.15	121	0.023	ABO, ACSL3, ACSS1, ADI1, ALDH1A1, ALDH3A1, ALDH6A1, ALDH7A1, AMACR, ATP6V1G2, BCKDHB, CA14, CA6, CKMT2, COX4I2, COX6B2, CYP17A1, CYP26A1, CYP2B6, DBT, EPHX2, EXTL2, FECH, GART, GLS, GLUL, GSTT2, GSTT2B, HIBCH, HSD17B8, HYAL4, ITPKB, KMO, LRTOMT, MOCS1, NAPRT, NMNAT1, OGDHL, PANK1, PCK1, PIGV, PLCB4, PLCD1, PLCD4, PMVK, SC5D, SELENBP1, SQOR, ST6GAL1, ST6GALNAC5, TBXAS1, TK2	ACSL4, ACSL5, ADA, ADCY3, ADPGK, AK2, AKR1B1, ALOX5, AMPD2, ARG2, BDH1, CA13, CANT1, CERS4, CHDH, CHST8, CMPK2, DEGS1, DGKD, DHRS9, ENPP4, ETNK1, EXT1, FPGS, FUT2, FUT3, FUT6, GALNT1, GALNT14, GALNT5, GCNT3, GNE, GYS1, HK2, HPRT1, IDH2, IDO1, IDS, INPP5E, ITPKA, LAP3, MGAT4B, MGST1, MRI1, MTHFD1L, MTHFR, NAMPT, NDUFA4L2, NT5C3B, NT5E, PDE4A, PDE7B, PDE9A, PFKFB3, PGS1, PIGL, PLB1, PRPS1L1, PRPS2, PTGS1, PTGS2, QPRT, RIMKLA, SGSH, SPTLC1, SYNJ2, TYMS, UGT1A6, UGT1A8
MAPK signaling pathway	18.7	37	0.021	CACNA1H, CACNA2D3, FGF9, FGFR2, FGFR4, HSPA1A, HSPA1B, MAP2K6, MAP3K20, MAP4K2, MAPK10, NF1, NLK, NTRK2, RASGRF1, RPS6KA5	AKT2, ARAF, AREG, CDC25B, DUSP4, DUSP8, GNA12, IL1B, IL1R1, IL1RAP, MAP3K14, MAP3K8, MAPK9, MAPT, MAX, PAK1, PGF, PPP3CC, RAC2, RELB, VEGFA
Focal adhesion	12.5	28	0.010	COL4A4, COMP, LAMA4, MAPK10, RASGRF1, THBS4	ACTN1, AKT2, CCND2, COL4A2, COL9A3, ITGA5, ITGB6, LAMA1, LAMC3, MAPK9, MYLK, PAK1, PGF, PPP1R12B, PXN, RAC2, SHC1, SPP1, THBS1, VAV1, VEGFA, ZYX
TNF signaling pathway	7.14	18	0.021	AKT2, BCL3, CXCL1, CXCL2, CXCL5, IL15, IL1B, IRF1, MAP3K14, MAP3K8, MAPK9, MMP9, PTGS2, SOCS3	CFLAR, MAP2K6, MAPK10, RPS6KA5
IL-17 pathway	6.16	15	0.037	MAPK10	CXCL1, CXCL2, CXCL5, DEFB4A, IL17RA, IL1B, LCN2, MAPK9, MMP9, PTGS2, S100A8, S100A9, TRAF3IP2, USP25
Prolactin signaling pathway	4.52	13	0.025	CYP17A1, ESR1, MAPK10	AKT2, CCND2, IRF1, JAK2, MAPK9, SHC1, SOCS1, SOCS2, SOCS3, STAT1

The heat maps in Figures 2A and 2B show the biological variance of the samples included in this study. It can clearly be seen that there is a significantly higher variance between the individual samples in the KC group compared to the control group. These differences are particularly pronounced in the stroma samples. Several candidates were selected from the top 20 deregulated transcripts of epithelial and stromal samples for RT-qPCR validation. The FCs calculated from the results of qPCR analysis were comparable to the FCs calculated from microarray data (Fig. 3).

Figure 2.

Heat maps of the differentially expressed genes (DEGs), including the 20 most upregulated and 20 most downregulated genes of the epithelial and stromal group. (A, B) Heat maps of the 40 most differentially expressed genes of the epithelial group (EKC vs. EN, A) and the stromal group (SKC vs. SN, B). (C, D) Heat maps of the differentially expressed miRNAs of the epithelial samples (EKC vs. EN, C) and stromal samples (SKC vs. SN, D), including 20 upregulated (red) and 20 downregulated (blue) gene hubs. Rows are identified by Entrez gene names or names of miRNAs. Intensities are shown by a color range, from red (row max) to white (row average) to blue (row minimum). The graphic was assessed by the versatile matrix visualization and analysis software Morpheus (https://software.broadinstitute.org/morpheus). EKC, keratoconus epithelium; EN, normal epithelium; SKC, keratoconus stroma; SN, normal stroma.

Figure 3.

qPCR validation of deregulated mRNAs in epithelial and stromal samples of patients with KC. The genes for validation were randomly selected from the top 20 deregulated genes of the microarray analysis. Of the selected genes, NANOG and SFRP1 were upregulated (5.3-fold and 17.4-fold, respectively) in the epithelial samples, and AKT3 and DHRS9 showed downregulation (1.6-fold and 4.5-fold, respectively). The CYP26A1 and DAPL1 genes showed upregulation in the stromal samples (46-fold and 11.6-fold, respectively) and downregulation in the KRT16 and DHRS9 genes (10.5-fold and 55.5-fold, respectively). Statistical analysis was performed using t-tests; n.s., not significant. *P < 0.05; **P < 0.01. E, epithelial; S, stromal; N, normal.

miRNA

The most downregulated miRNAs in KC epithelium or KC stroma were miR-634, with a FC of –11.4, and miR-936, with a FC of –37.4. The most upregulated miRNAs in epithelium or stroma were miR-211-3p (FC = 59.8) and miR-135b-5p (FC = 33.9), respectively. A complete list of deregulated miRNAs in the epithelium and stroma is shown in Supplementary Table S6. The target genes of the top 20 deregulated miRNAs were assessed using miRTargetLink (https://ccb-web.cs.uni-saarland.de/mirtargetlink). Only the targets with strong evidence were used for further evaluation. A detailed table is displayed in Supplementary Table S2. A total of 58 genes were identified as targets of two or more deregulated miRNAs. The biological variance of miRNA expression in the samples is shown in Figures 2C and 2D using a heat map. Like the situation in the mRNA analysis, miRNA expression also showed a significantly higher interindividual variance in the KC group compared to the control group. Especially outstanding is the opposite regulation between normal and KC samples in the stroma. The two miRNAs with the most validated targets were miR-141-3p and miR-429, each with 24 or 23 assigned target genes. Several targets of the top 40 deregulated miRNAs were selected for RT-qPCR validation in epithelial and stromal samples. The FCs calculated from the results of the qPCR analysis were comparable to the FCs calculated from microarray data (Fig. 4). To obtain further information about miRNA–mRNA interactions with potential relevance in KC, we used miRTargetLink (https://ccb-web.cs.uni-saarland.de/mirtargetlink) to identify experimentally validated miRNA–mRNA pairs for the 20 most up- and downregulated miRNAs in the stroma and epithelium. Table 6 lists all targeted mRNAs with significant opposite deregulation based on our microarray data (i.e., downregulated mRNAs that are targets of upregulated miRNAs and vice versa). In epithelium, no interactions of either upregulated miRNAs with downregulated mRNAs or downregulated miRNAs with upregulated mRNAs were found. Although numerous interactions were found for upregulated miRNAs in the stroma, only one interaction could be verified for downregulated miRNAs—namely, the interaction of miR-198 with NTRK3 mRNA.

Figure 4.

qPCR validation of deregulated miRNAs in epithelial and stromal samples of patients with KC. The miRNAs for validation were randomly selected from the top 40 deregulated miRNAs of the microarray analysis. Of the selected miRNAs, miR-101-3p and miR-151-3p were upregulated (4-fold and 2.7-fold, respectively) in the epithelial samples, and miR-320b showed downregulation (1.6-fold). MiR-125-5p showed no change in the epithelial samples of patients with KC in contrast to downregulation by microarray analysis (onefold). In stromal samples of lenticules, miR-205-5p and miR-200c-3p showed upregulation (3.2-fold and 2.6-fold, respectively), and miR-320a and miR-134 showed downregulation (2.3-fold and 3.0-fold, respectively). Statistical analysis was performed using t-tests. *P < 0.05; **P < 0.01.

Table 6.

Analysis of Deregulated miRNAs and Their Interactions With mRNAs in Epithelium or Stroma of KC Subjects

Upregulated miRNAs in KC Stroma	Downregulated Target Genes in KC Stroma	Fold Change
miR-141-3p	MAPK9	−2.75
	E2F3	−4.11
	SHC1	−3.42
	TM4SF1	−3.70
	CDC25A	−3.45
miR-210-3p	P4HB	−2.09
	E2F3	−4.11
	IGFBP3	−50.03
	COL4A2	−2.54
	HIF1A	−2.76
	BNIP3	−4.19
miR-126-3p	PLK2	−5.00
	SLC45A3	−3.42
	VEGFA	−7.95
	TWF1	−2.06
	PITPNC1	−2.17
	MERTK	−3.61
	EGFL7	−2.23
	ADM	−4.00
miR-15a-5p	CDC25A	−3.45
	CCND2	−3.94
	VEGFA	−7.95
	HMGA1	−3.53
miR-205-5p	E2F5	−3.54
	E2F1	−6.74
	VEGFA	−7.95
	ACSL4	−3.29
	ITGA5	−6.02
	HMGB3	−3.11
miR-429	SHC1	−3.42
miR-96-5p	PRMT5	−2.80
	SCARB1	−2.54
	SLC6A6	−2.61
miR-183-5p	IDH2	−2.95
	NFIL3	−5.25
miR-194-5p	SOCS2	−4.71
miR-135b-5p	MID1	−2.26
Downregulated miRNAs in KC Stroma	Upregulated Target Genes in KC Stroma	Fold Change
hsa-miR-198	NTRK3	2.54

Discussion

In the present study, eight corneas, each from patients with KC, divided into epithelium and stroma, were compared with non-KC samples with respect to mRNA and miRNA expression using microarrays. To the best of our knowledge, this is the first study in which the expression profiles of mRNAs and miRNAs in corresponding epithelial and stromal cells of KC corneas have been observed.

There are several studies that have identified some candidate genes by RNA sequencing (RNA-Seq) from whole corneas or epithelium of KC corneas or corneal fibroblasts from cell cultures. Although different study materials were used and different bioinformatics analysis methods were applied, there is some overlap of deregulated genes. Among the top 20 up- and downregulated mRNAs from the current study, 12 mRNAs were previously found deregulated in KC in the studies of Kabza et al., You et al., Khaled et al., and Sharif et al., with each of these mRNAs identified in at least two different studies. From the present study, the mRNAs KRT16, CXCL1, and SA1009 from the top 20 deregulated stromal mRNAs overlap with the top 20 from the study by Khaled et al. (Supplementary Table S4).

The present study shows significant differences in the expression of deregulated mRNAs in the epithelium and stroma of patients with KC. In the epithelium, only 170 differentially expressed mRNAs were found with a maximum eightfold expression change, but the number of deregulated mRNAs and their maximum expression change are distinctively higher in the stroma, with 1498 deregulated genes and a maximum FC of >300.

mRNA Analysis

Among the top 20 deregulated mRNAs in the epithelial samples of Patients with KC, phospodiesterase B8 (PDEB8) was found with an FC of –8.7 and secreted frizzled-related protein 1 (SFRP1) with an FC of 7.2. PDEB8 induces cell proliferation via cAMP signaling pathway. Reduced expression of PDEB8 would therefore also be indicative of reduced cell proliferation as observed in KC. SFRP1, which is significantly elevated with an FC of 7.2, acts as an inhibitor of Wnt signaling with an antiproliferative effect. Again, it can be speculated that cell proliferation is disturbed in KC. On the one hand, expression of the cell proliferation–promoting gene (PDEB8) is decreased, but, at the same time, expression of an antiproliferative gene (SFRP1) is increased.

The top 20 deregulated mRNAs in stroma of Patients with KC include dehydrogenase/reductase (SDR family) member 9 (DHRS9) with a FC of –311.7 and CYP26A1 with a FC of 54.9. DHRS9 is an important molecule in the retinoic acid pathway, and it is involved in the biosynthesis of all-trans-retinoic acid (atRA), which in turn induces the tumor-suppressing gene XAF1, which is important due to its inhibition of cell proliferation, induction of cell differentiation, and apoptosis. On the other hand, CYP26A1 is highly upregulated in KC stroma. CYP26A1 is the key regulator in RA catabolism. Considering these two molecules, one could speculate that less RA is available due to decreased DHRS9, but at the same time more RA is degraded due to increased CYP26A1. In addition, RA also exhibit antioxidative activity. As already mentioned, the antioxidants in KC are reduced.

Also worth noting is the number of deregulated small-proline rich proteins (SPRRs) among the top 20 deregulated mRNAs in the stroma. The SPRRs belong to a protein family with four members (SPRR1–4), of which SPRR 1 to SPRR 3 show several isotypes. SPRRs are known to promote epidermal differentiation. Little is known about their function in other cell types, but they have been demonstrated in smooth muscle cells and fibroblasts, including their association in fibrosis progression. Among the top 20 downregulated mRNAs in the stroma, five belong to the SPRRs that are closely associated with keratins. Thus, studies indicate that an increase in SPRR may lead to increased formation of extracellular matrix (ECM), which is characteristic in fibrosis. Conversely, this could mean that a decrease in SPRR could lead to decreased collagen formation and reduced corneal cells proliferation, as observed in KC.

The SPRRs are also associated with inflammatory processes, and their expression is elevated in psoriasis and allergic lung inflammation. Similarly, elevated SPRR levels have been demonstrated in dry eye disease in response to proinflammatory cytokines. KC has long been classified as a non-inflammatory condition, but an inflammatory component is now beyond question. Interestingly, none of the genes associated with inflammatory response or regulation is upregulated in stromal cells of Patients with KC, but several mRNAs associated with inflammatory response or regulation are found among the 915 downregulated genes, including S100A8, S100A9, and KRT16 among the top 20 downregulated genes. Further downregulated genes were IL-1β, IL-15, and RELB. This change may be attributed to a counter-regulation, its analysis at protein level would be required for future studies. Another observation was the high interindividual variance of mRNA expression, especially in the data for stromal cells. A possible explanation for this phenomenon may be the individual acute, active changes in corneas of Patients with KC which cannot be captured by the clinical data, considered in this study.

miRNA Analysis

miRNAs are considered key factors in the regulation of gene expression. The group of Drewry et al. identified several miRNAs associated with KC in the ciliary body, cornea, and trabecular meshwork using small RNA-Seq. In their study, miR-143-3p, miR-182-5p, and miR-92a-3p were identified as highly expressed miRNAs. In the present study, we could not detect miRNA-143-3p in either the epithelium or the stroma. The miRNAs miR-182-5p and miR-92a-3p were also not expressed in the epithelium, but their expression was significantly increased in the stroma of Patients with KC, with FCs of 2.79 and 2.079, respectively. These discrepancies may reflect the use of different methods in miRNA detection, such as small RNA-Seq versus microarrays, or differences in normalization or filtering procedures.

Members of the miR-200 family were co-expressed and upregulated in KC stromal tissue in the current study. The miR-200 family includes miR-200a/-200b/-429 and miR-200c/-141. Members of the miR-200 family are known to regulate epithelial–mesenchymal transition and were found downregulated in, for example, pterygium samples. On the other hand, upregulation of members of the miR-200 family is associated with oxidative stress, induced by reactive oxygen species (ROS). It is well known that ROS are elevated in KC. Furthermore, overexpression of members of the miR-200 family represses proliferation and migration of primary epithelial cells. Of the top 20 upregulated miRNAs in the current study, miRNA-141-3p and miRNA-429 were highly upregulated, with FCs of 33.24 and 17.28, respectively. The most downregulated miRNA was miR-936 with a FC of −37.14. Downregulation of miR-936 leads to reduction of proliferation in laryngeal squamous cell carcinoma. This would suggest that the proliferation rate of KC stromal cells is reduced and thus fewer cells are available to produce ECM components. The miR-204 and miR-126-3p associated with corneal neovascularization do not appear to play a decisive role in KC, although miR-204-5p and miR-126-3p are both significantly overexpressed in stromal cells relative to the normal corneas, with FCs of 10.05 and 12.87, respectively.

Pathway analysis using GeneTrail 3.0 showed an enrichment of deregulated genes in metabolic signaling pathways in both epithelial and stromal cells of patients with KC. The fact that metabolic signaling pathways are affected has already been described several times. This can be explained by the fact that oxidative stress is increased in KC. Increased oxidative stress upsets the balance of the cells, which affects the metabolic functions. This result supports the data of Karamichos et al., who were also able to identify a significant change in different metabolic signaling pathways in their study, including urea cycle, citric acid cycle, and oxidative stress. However, the question remains as to whether the metabolic alteration is a consequence of the disease or the disease is a consequence of the metabolic alterations.

Another highly altered signaling pathway in the stroma is the focal adhesion signaling pathway. In focal adhesion, those molecules play a role, which are crucial as signal transmitters for the ECM interaction. The ECM plays a major role in stability and transparency of the cornea, and it is known that collagen synthesis in the ECM of Patients with KC is altered. Changes in focal adhesion have also been described by Nowak-Malczewska et al. The cornea of patients with KC shows thinning and an associated conical bulge in clinical expression. The progressive thinning is caused by, among other factors, the altered production of the ECM. A key molecule could be heat shock protein 47 (Hsp47), a chaperone responsible for the correct folding and ejection of collagens. Although our data indicate no change in Hsp47 (SERPIN1) transcript expression, the group led by Foster et al. demonstrated a significant decrease in Hsp47 protein levels. Thus, the collagens would not be properly shuttled out of the cells and would accumulate in them, which may lead to decreased ECM production in the stroma. It would be of particular interest to determine whether the data obtained here could be linked to the severity or current inflammatory activity of KC. Nevertheless, corneas from patients with KC were exclusively processed in the present study from patients in advanced stages of the disease leading to corneal transplantation.

Conclusions

The results of the current study show a significantly higher level of perturbation of the gene expression profile in the stroma than in the corneal epithelium of patients with KC. Predominantly affected are genes that show a connection to the regulation of inflammatory processes and affect metabolic pathways. The results of this study show that different mRNAs are deregulated in the epithelium and stroma of patients with KC. We could also show an altered miRNA expression pattern in epithelium and stroma of patients with KC. Based on these results, further studies could provide insight into the extent that different cell types influence each other.