Keratoconus (KTCN) is the most common corneal ectasia, characterized by progressive stromal thinning and the presence of the keratoconic cone (1). The structural changes, which often manifest asymmetrically, lead to a loss of visual acuity as an effect of irregular astigmatism (2). The first symptoms usually occur in adolescence, and the disease progresses through the third or fourth decade of life (1). Keratoconus affects both genders, but the reason of its higher occurrence in male patients remains unclear (3). The occurrence of KTCN may vary depending on the ethnicity of patients (1), with an estimated prevalence of 1.38 per 1,000 in a general population (4).

Advanced KTCN remains one of the most common indications for full-thickness (penetrating keratoplasty) and partial-thickness (deep anterior lamellar keratoplasty) corneal transplantation in developed countries (5). In the face of the organ shortage crisis in worldwide transplantation, there is a need for other therapies that inhibit the progression of KTCN to an advanced stage. A corneal collagen cross-linking (CXL) is an effective alternative used to stabilize the progression of KTCN and often improves the quality of patients’ vision (6).

Histopathological abnormalities are present in all corneal layers of the KTCN cornea (7), with diverse manifestations in central and peripheral zones being an effect of different biomechanical tension across corneal curvature in KTCN (8). Changes in collagen fibers and other elements of the extracellular matrix (ECM) form stromal irregularities (8) contributing to the KTCN cone formation and also affect the corneal epithelium (CE). The resulting rebuilt epithelial structure is characterized by a central thinning surrounded by a thickened annulus, named doughnut pattern (9). In the morphological observations of the CE (including confocal and immunofluorescence microscopy), the decreased cell density and their elongated shape, disruptions in basal integrity, increased number of apoptotic cells, and depositions of iron particles have been previously described (7, 10, 11). Moreover, some molecular proteomic and transcriptomic findings showed disruption of elements of Wnt signaling (e.g., WNT10A) (12), cell–cell communications (e.g., CDH13) (13), pro-inflammatory cytokines (e.g., IL6) (14), matrix metalloproteinases (e.g., MMP9) (15), and innate immune system (e.g., TLR2 and TLR4) (16).

The available reports regarding molecular and environmental findings in KTCN studies (17–22) indicate the multifactorial nature of the disease. The behavioral, environmental, and socioeconomic factors supposedly induce disease manifestation/progression in genetically predisposed individuals. In clinical studies, KTCN was found to be related to chronic eye rubbing, allergy, atopy, asthma, and UV exposure (18).

Historically, KTCN has been described as a non-inflammatory disease (2), but in the last decade, this aspect has been widely studied and questioned by ophthalmologists. However, the currently valid Global Consensus on Keratoconus and Ectatic Diseases (23) has not addressed this issue.

Previously published transcriptomic and proteomic findings derived a hint for the influence of inflammatory factors in KTCN as the increased levels of inflammatory cytokines (as IL6, TNFα, and MMP9) were found in patients’ tear fluid and serum samples and suggested to contribute to the thinning and weakening of the corneal tissue in KTCN (17, 24–26).

So far, the genetic aspects regarding the functioning of the immune system in the KTCN have not been demonstrated. Numerous genetic findings, including familial inheritance (27, 28), a concordance between monozygotic twins in contrast to dizygotic twins (29), and the occurrence of syndromic KTCN (30), imply the presence of a genetic component for the disease development. The previous results obtained by our research group confirmed the postulated genetic heterogeneity between patients and the involvement of numerous genetic factors in the pathogenesis of KTCN (19, 31–33). Besides the identification of the 13q32 (31), 5q31.1-q35.3 (32), 2q13-q14.3, and 20p13-p12.2 (33) KTCN loci, the numerous sequence variants in VSX1, TGFBI, DOCK9, STK24, IPO5, SKP1 and IL17B candidate genes (19, 32) have been demonstrated. So far, whole-genome sequencing (WGS) has not been applied in KTCN research. As only WGS provides an overview of the entire human genome, and it is a suitable method for causative variant discovery in genetically heterogeneous diseases, this approach is a natural next step in investigating the pathogenesis of KTCN. We hypothesize that variants in non-coding regions of the genome complement the previously identified features specific to the KTCN exome.

The multilevel structural and functional changes in cornea, including the varied morphology and various intensities of inflammatory-causing internal biomechanical tension and external environmental stimuli (e.g., eye rubbing) acting on particular zones of CE, suggest the diverse molecular features across corneal curvature. These premises support the designation and separation of three topographic regions of the CE, namely, central, middle, and peripheral. Here, we assessed the genomic diversity across keratoconic cone surface implementing a unique study design embracing the three topographic regions of the CE of adolescent patients with KTCN.

The complex background of KTCN, in addition to environmental factors, includes heterogeneous genetic characteristics. The multiple candidate genes and variants identified using different molecular techniques indicate that this heterogeneity requires the introduction of high-throughput investigation methods, which may be more effective in unraveling the KTCN-specific determinants.

Each individual genome differs from the reference human genome by an average of 3.5 million SNPs, but only a minority of those variations contributes to the phenotypic differences and disease predisposition (53). Only WGS can provide an overview of the entire human genome, allowing the global analysis and the full set of variants to be assessed (54–56). The majority of SNVs and indels is identifiable by both WGS and exome sequencing (ES), but still approximately 3% of coding variants is missed in the ES approach (54). Therefore, WGS is more efficient than ES for detecting coding sequence variants, besides the undoubted benefit of identification of non-coding sequence variation. To date, WGS has not been applied in KTCN research. Since WGS is an effective method of detecting mutations in genetically heterogeneous diseases, including ophthalmic diseases (57), and is more efficient than ES in the identification of sequence variation, this technique is crucial in deciphering the genetic aspects of KTCN.

Here, we found that the genetic features of patients with KTCN correspond with those previously reported regarding sequence variants in components of ECM (19, 50). Also, we identified other sequence variants in genes previously reported in KTCN research in the same population (19), including COL17A1 [playing a role in the integrity of hemidesmosome and the attachment of basal keratinocytes (58), whose expression was decreased in corneal epithelial cells in single-cell RNA-seq analysis of KTCN corneas (59)]; FREM2 [protein required for maintenance of the integrity of the skin epithelium (60) and eye morphogenesis (60)]; SSX2IP [which may connect the nectin–afadin and E-cadherin–catenin system through alpha-actinin and may be involved in the organization of the actin cytoskeleton (61)]; MINK1 [mediating stimulation of the stress-activated protein kinase MAPK14/p38 MAPK downstream of the Raf/ERK pathway (62)]; CAMSAP1 [Calmodulin-regulated spectrin-associated protein 1; key microtubule-organizing protein (63)]; and RGS12 [Regulator of G-protein signaling 12 (64)]. Moreover, we identified variants previously reported in EML6 (50), which is an element of cytoskeletal network and microtubules in cultured keratocytes (50, 65); MAGI1 (50), which regulates cell–cell and cell–matrix adhesion (66); and COL5A1 (48, 49), whose knockout in corneal stroma in mouse showed severe dysfunctional regulation of fibrillogenesis resulting in decreased fibril number and a disorganized lamellae structure with decreased stroma thickness (67).

Also, we have revealed variants in genes involved in keratin-related pathways (namely, the “keratinization” pathway from the Reactome database), such as EVPL, PCSK6, KRTAP10-7, KRT35, FLG, and KRT33B. As keratins participate in forming intermediate filaments (IF), which provide mechanical support, and form desmosomes between cells and hemi-desmosomes with basement membranes for epithelium integrity (68), this result may be of special interest regarding changes in biomechanics of keratoconic cornea.

Importantly, non-coding variants are now considered to be potentially pathogenic (69, 70). Sequence variation in the form of SNVs or indels in the regulatory regions can either generate or abolish TF binding sites, modulating the transcriptional activity of genes and regulating chromatin state (71). The disruption of the TF binding site and its effect on the disease phenotype have already been reported in many cancer types such as T‐cell acute lymphoblastic leukemia (72) and esophageal and gastric cancers (73). However, while great advances in predicting the effects of coding variants have been made, the assessment of non-coding variants remains challenging. In particular, the immense number of variants revealed per patient sample impedes the analysis and biological interpretation. In our investigation, all variants located in REs, including promoters, promoter flanking regions, enhancers, CCCTC-binding factor binding sites, TFBS, or open chromatin regions of fibroblasts, were evaluated in silico towards functional effect (TFBS gain/lost) based on the motif analysis. The obtained list of RE variants showed unity on the pathway level with previously reported variants, namely, focal adhesion (19), Wnt signaling (19), and TGF-beta signaling (74). Moreover, the LAMC2, COL17A1, SDK1, PECAM1, TJP2, MINK1, NEIL3, EVPL, and PCSK6 genes were revealed to have unique variants for patients with KTCN in both coding and non-coding genome sequence, which again were found as involved in the ECM, cell–cell communications, cytoskeleton organization, and keratinization pathways. To further evaluate the effect of variants located in REs, the functional approaches such as gene reporter assay are needed, but firstly transcriptome data from adequate biological material is necessary to be integrated with the genomic data. Therefore, we evaluated the data overlap between DEGs in full-thickness KTCN corneas (45) and variants recognized in the REs. We found the involvement of RE variants in numerous elements of ECM and that fact supports our assumption concerning the role of RE sequence variation in KTCN pathogenesis.

Although single variants could influence the KTCN phenotype in some families (75), our genomic data (both coding and non-coding sequence variation) suggest that KTCN largely occurs through the interplay of multiple variants that interfere with major physiological processes in the cornea. The presented connection between coding and non-coding sequence variants in the gene ontology (“Actin cytoskeleton”, “Extracellular matrix”, and “Collagen-containing extracellular matrix”) and in the pathway assessments (“Focal adhesion”, “Hippo signaling pathway”, and “Wnt signaling”) corroborates this assumption, as well as the identified 13 regions enriched in variants in the 5q locus that had previously been found to be specific to KTCN in distinct populations (32, 46, 47).

It has been known that DNA structural variation modifies the interactions between REs and their target genes (71). Here, the majority of the deletions were found to be located in non-coding and intron regions. Overall, assessing the previously published KTCN-specific loci, we did not observe a difference in the number/density of deletions comparing patients with KTCN and control individuals. However, numerous genes involved in the innate immune system; adaptive immune system; ECM including collagen biosynthesis, formation, and trimerization; and others were found to contain deletions in size of 51–5,098 bp in patients with KTCN in contrast to control individuals. The extent to which this finding is important should be clarified in further research on this aspect. To date, deletions >50 bp in size themselves and their effect on the KTCN phenotype were only once evaluated, in an Ecuadorian family with KTCN (31), but that study did not reveal any copy number variation on evaluated 13q32 in the tested individuals.

We have anticipated that the various intensities of both internal biomechanical tension and external environmental stimuli (e.g., eye rubbing) acting on various topographic regions of CE could result in diverse molecular features across corneal curvature. Our study performed in the designated three topographic regions of the CE, central, middle, and peripheral, enabled the conclusion that there is no genomic heterogeneity across the corneal surface in KTCN. However, in a comprehensive morphological, transcriptomic, and proteomic evaluation performed in the same corneal samples, we identified the variability shaping the distinct features of CE zones (13). Continuing the topographic aspects of the CE, we have considered features of mosaicism to be involved in the process of remodeling of the KTCN CE. To date, based on the results of exome and/or Sanger sequencing performed in matched KTCN cornea–blood pairs, we found no evidence of inter-tissue variants in human KTCN corneas, since all variants possibly related to KTCN were identified in both corneal tissue and peripheral blood of patients (19). Again, here we did not identify features of mosaicism.

In our opinion, the most important findings of this study are those concerning the inflammation in KTCN. Nowadays, the inflammatory aspects in KTCN are widely discussed. The multiple analyses of the tear cytokine and protease profiles of patients with KTCN have shown a persistent inflammation (17, 76). Previously, we identified numerous, randomly distributed variants in genes that encoded inflammatory molecules (IL1A and IL1B) and one variant (c.214+242C>T in the IL1RN) presented in all KTCN individuals (33). Here, the enrichment in variants of genes involved in the innate (neutrophil degranulation) and adaptive (Class I MHC antigen processing and presentation) immune system was revealed. In our study, TRAV19 and TRVB12-5 genes were recognized with three recurring variants (one as a stop gain). Because these genes encode variable regions of T-cell receptors, it could indirectly influence the process of antigen processing (77). Importantly, in a single-cell study of full corneas, the upregulation of the T-cell receptor signaling pathway has been recognized (59). Identified recurring missense variants in KIR2DL1 and KIR3DL1 could influence the activity of the immune system, since they are members of the highly conserved killer cell immunoglobulin-like receptors (KIRs) family (78), which are transmembrane glycoproteins expressed by natural killer (NK) cells and subsets of T cells (79). It was previously reported that the variant in KIR3DL1 has affected the NK cell function (80). NK cells activate neutrophils and induce both their infiltration to the inflamed sites and degranulation (81). Also, NK cells were reported to participate in acute immune response in dry eye syndrome and trigger Th-17 response, manifested by an increase of IL17 (82). Moreover, the involvement of NK cells in the corneal wound healing has been confirmed in a murine model (83). The evidence for migration of cells into corneal limbus and central stroma after a  2-mm epithelial abrasion of the mouse cornea and altered inflammatory reaction after antibody-induced depletion of NK cells were presented (83). Previously, as elements of innate immune responses, the expression of TLR2 and TLR4 has been studied by flow cytometry in corneal epithelial cells and blood samples of KTCN patients and their relatives (16, 84, 85), and significant upregulation has been shown, indicating a potential of the TLR2 and TLR4 proteins to be the KTCN biomarkers. However, no sequence variants in the TLRs were recognized here as well as in our current gene and protein expression assessment (13). Instead, we recognized sequence variants of other innate immune system elements (e.g., MUC5AC, as mucins consist the first-line defense). Moreover, the increased proportions of activated NK cells as well as neutrophils and γδT cells have been reported in KTCN studies (86). Also in many studies of tear fluid collected in KTCN patients, the higher levels of interleukins (IL1β, IL6, and IL17A), interferons (IFNα/β/γ), and other pro-inflammatory factors (TNFα, TGFβ1, LIF, granzyme-B, perforin, and MMP2) were recognized (14, 24, 86–88). Besides published transcriptomic and proteomic findings, the genomic aspects regarding disruption of immune responses have not been demonstrated in KTCN research. Also, the here-identified single high-impact sequence variants (CEP290, MMRN1, and FN1) and missense variants (PECAM1, VWF, CD109, HRG, TLN1, HPSE, EPO, SIGLEC10, DOCK8, ERBB3, and AP3B1) could contribute to molecular dysregulation of wound healing, followed by inflammation, although we have not experimentally determined their impact.

Therefore, we conclude here that the disturbance of the immune responses in KTCN affects various elements of corneal hemostasis, including the antigen presentation and neutrophil degranulation processes, and possibly results in thinning and weakening of the KTCN corneas.

Summarizing, applied WGS allowed for identification and confirmation of numerous genomic features in KTCN. In future studies, WGS data should be integrated with other genome-wide analyses data, including chromatin accessibility (Assay for Transposase Accessible Chromatin with sequencing, ATAC-Seq), transcriptomics, and proteomic data, in order to obtain the most reliable research results. In this study, owing to a limited number of included patients, we intently did not focus on particular sequence variations but rather on the sets of impacted genes and the pathways in which those genes are involved. Additionally, we decided to evaluate adolescent patients with KTCN only, assuming the more pronounced genetic background of the disease. The recruitment of the control group was especially challenging, as nowadays the photorefractive keratectomy (PRK) procedure, during which the CE is removed, is rarely performed and is replaced with the newer techniques of vision correction without CE removal. The youngest possible patients with minimal vision impairment undergoing the PRK procedure were ascertained into the control group to minimize potential bias. The PRK procedure is performed in adults, and this fact also results in the difference in age between the compared groups of patients. Although only four patients and two controls were involved, 18 samples were tested in the WGS experiments and data obtained from those 18 samples were carefully curated. In addition, the sequence variant interpretation was performed using additional allele frequency data derived from the reference databases (1000 Genomes, ESP and gnomAD). Owing to the small size of the compared study groups, we were unable to subdivide these groups into subgroups of patients, taking into account the presence of an allergic disease or dry eye syndrome.

The multiple coding and non-coding sequence variants were found in genes (or in their close proximity) contributing to the previously discussed KTCN-specific cellular components, and newly presented aspects of innate and adaptive immune system responses, pointing to the involvement of inflammatory aspects in KTCN. Therefore, we conclude that variants in non-coding regions of the genome have complemented the previously identified features specific to the KTCN exome.

Additionally, the 35 chromosomal regions shown enriched in both coding and non-coding KTCN-specific sequence variants confirmed the heterogenous background of KTCN with a limited number of common elements.

We did not find the diversity in genetic features of the assessed topographic regions of CE. This lack of genomic heterogeneity across the corneal surface and no evidence of inter-tissue variants in human KTCN corneas suggest other molecular mechanisms of KTCN-specific changes in CE.

The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/clinvar/ under the accession number SUB13057607.

The studies involving human participants were reviewed and approved by Bioethics Committee at Poznan University of Medical Sciences, Poznan, Poland. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

KJ: participated in research design, sample collection, and patient surveys; performed sample preparations towards WGS experiments; participated in bioinformatics analyses; figures and tables preparation; participated in the genomic and clinical data interpretation; and wrote the manuscript. MM-K: performed recruitment of the patients and clinical evaluation; facilitated sample collection; performed the clinical data interpretation; participated in the genomic data interpretation; and manuscript writing. MK: performed setup of bioinformatics pipeline; bioinformatics analyses of the genomic data; and preparation of figures and tables. JAK: participated in the genomic data interpretation. MG: research design; secured research funding; participated in the genomic data interpretation; and edited the manuscript. All authors contributed to the article and approved the submitted version.