Sequence variants contributing to dysregulated inflammatory responses across keratoconic cone surface in adolescent patients with keratoconus

Background Keratoconus (KTCN) is the most common corneal ectasia resulting in a conical shape of the cornea. Here, genomic variation in the corneal epithelium (CE) across the keratoconic cone surface in patients with KTCN and its relevance in the functioning of the immune system were assessed. Methods Samples from four unrelated adolescent patients with KTCN and two control individuals were obtained during the CXL and PRK procedures, respectively. Three topographic regions, central, middle, and peripheral, were separated towards the whole-genome sequencing (WGS) study embracing a total of 18 experimental samples. The coding and non-coding sequence variation, including structural variation, was assessed and then evaluated together with the previously reported transcriptomic outcomes for the same CE samples and full-thickness corneas. Results First, pathway enrichment analysis of genes with identified coding variants pointed to “Antigen presentation” and “Interferon alpha/beta signaling” as the most overrepresented pathways, indicating the involvement of inflammatory responses in KTCN. Both coding and non-coding sequence variants were found in genes (or in their close proximity) linked to the previously revealed KTCN-specific cellular components, namely, “Actin cytoskeleton”, “Extracellular matrix”, “Collagen-containing extracellular matrix”, “Focal adhesion”, “Hippo signaling pathway”, and “Wnt signaling” pathways. No genomic heterogeneity across the corneal surface was found comparing the assessed topographic regions. Thirty-five chromosomal regions enriched in both coding and non-coding KTCN-specific sequence variants were revealed, with a most representative 5q locus previously recognized as involved in KTCN. Conclusion The identified genomic features indicate the involvement of innate and adaptive immune system responses in KTCN pathogenesis.


Introduction
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)(18)(19)(20)(21)(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, TNFa, 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][25][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)(32)(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.
were involved in this study. The study protocol was approved by the Bioethics Committee at Poznan University of Medical Sciences, Poznan, Poland. The possible consequences of the study were explained, and informed consent was obtained from all participants, according to the Declaration of Helsinki.
Each individual underwent a complete ophthalmological examination, including the assessments of both uncorrected (UCVA) and best-corrected visual acuity (BCVA), intraocular pressure (IOP), corneal tomography with rotating Scheimpflug camera WaveLight ® Oculyzer ™ II (Alcon, TX, USA), epithelial thickness mapping [Spectral-Domain Optical Coherence Tomography (SD-OCT) device, Zeiss Cirrus 5000, Carl Zeiss Meditec, Dublin, California, USA], and slit-lamp and dilated funduscopic examination. A questionnaire comprising the behavioral, environmental, and socioeconomic aspects, including eye rubbing, use of contact lenses, atopy, UV exposure, smoking, reading habits, time spent in front of a screen, hormone intake, education level, and place of living, was completed by each participant.
The inclusion and exclusion criteria for adolescents with KTCN, and control individuals are described in detail in Supplementary Methods 1.1.

CXL and PRK procedures
CXL in patients with KTCN was performed in accordance with the standard Dresden Protocol (34), while the PRK was performed as a refractive error correction procedure (35) in control individuals as described in Supplementary Methods 1.2, 1.3 respectively.

Material collection and sample preparation
Stamps towards the nose and eyebrow were made on the CE before the excision in the CXL/PRK procedures to enable correct tissue orientation during cutting and separation of the topographic regions. The obtained tissues were submersed in an RNA stabilization solution (RNAlater; Qiagen, Hilden, Germany) immediately after excision and stored at −80°C until further proceeding.
The procedure of designation of the particular topographic region is shown in Figures 1, S1. The details of the sample preparation are described in Supplementary Methods 1.4.

DNA extraction
Separated CE samples were transferred from the microscope slides to the lysis solution (Norgen Biotek, Thorold, ON, Canada). Total RNA, DNA, and proteins were extracted and purified according to the instructions of the RNA/DNA/Protein Purification Plus Micro Kit (Norgen Biotek). The quantity of DNA samples was measured by Qubit dsDNA HS Assay Kit (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) and quality was assessed by 1% gel electrophoresis.

Library preparation, sequencing, and WGS data analyses
WGS of the 18 CE samples was performed with a TruSeq Nano DNA HT Library Prep Kit (Illumina, San Diego, CA, USA) and the HiSeqX platform (Illumina) with mean coverage depth 30× at CloudHealth Genomics (Shanghai, China) as previously described (36). WGS data processing is summarized in the Supplementary Methods 1.5.
The following subset of sequence variants were analyzed: high impact variants (nonsense mutations/stop gain, stop lost, frameshifts, splice site mutations, etc.), missense variants (missense mutations considered deleterious by either SIFT or PolyPhen), and variants in regulatory regions [overlapping promoters, promoter flanking regions, enhancers, CCCTC-binding factor binding sites, transcription factor (TF) binding sites, or open chromatin regions, based on the Ensembl Regulatory Build for fibroblasts], classified as single-nucleotide variant (SNV)/deletion/insertion/indel/sequence alteration (substitution).
Motif analysis of variants located in regulatory elements (REs) was executed using the web-based FABIAN-variant application (37), and SNVs were additionally verified using MotifbreakR (38).
A set of impacted genes (high impact/missense/regulatory regions variants) was assessed for each patient after removing variants identified in control patients and variants with a MAF_AF (Maximum observed allele frequency in 1000 Genomes, ESP, and gnomAD) ≥0.1 (1000 Genomes Project phase three, ESP6500SI-V2, gnomAD v2.1). The non-redundant set of impacted genes in KTCN patients was analyzed using ConsensusPathDB (39), ShinyGO (40), STRING tool (41), and Reactome database (42). In pathway enrichment and gene ontology analyses, the p-value ≤ 0.05 and/or FDR ≤ 0.05 were the cutoffs.
The identification of genomic regions significantly enriched in genes with KTCN-specific sequence variants, compared to the density of genes in the background, was performed using ShinyGO (40). A "sliding window" is applied to scan the genome and multiple hypergeometric tests to determine whether the evaluated genes are significantly overrepresented within the analyzed window with a size of 6 Mbp and the FDR ≤ 0.05 as a cutoff.
The Integrative Genomics Viewer (IGV) (43) was used for visual exploration of genomic data.
To determine if there is any difference between the genomic characteristics of the particular CE topographic regions and fullthickness corneas, we compared the data obtained for the CE samples with our previously generated molecular data for the same population (Polish Caucasians) (19,32,44,45) and with other published data on candidate genes and gene pathways (31,(46)(47)(48)(49)(50)(51).
All experimental samples were analyzed separately and at no step of the analysis (including the analysis of raw next generation sequencing reads) were data from the three regions of one individual combined. In the Results section, the variant(s) referred as recognized in the KTCN patient were present in all three topographic regions.

Data integration
Simultaneously, the transcriptomic and proteomic assessments as well as morphological evaluation in the same CE samples were performed, as described in detail elsewhere (13). All the data of the mentioned investigation were considered in final interpretation of the study outcomes.

Characteristics of patients and DNA samples
Four unrelated adolescent patients with KTCN (one female/ three male patients) and two control individuals (two male individuals) were involved in this study. There was no history of KTCN and genetic diseases (including autoimmune diseases) in families of ascertained patients and controls. Clinical characteristics of the examined individuals and the eyes subjected to the surgery are presented in Table 1, while the information collected for both eyes in the studied individuals is compiled in S1 Table. Selected behavioral, environmental, and socioeconomic aspects evaluated in the questionnaire are presented in S2 Table. A total of 18 experimental samples of CE (three topographic regions in each of the six ascertained individuals) were collected and used in the WGS experiments, in accordance with a study scheme (Figures 1, S1). The results of quantity and quality control of the 18 DNA samples and quality control of WGS reads are presented in S3 Table. 3

.2 Genes shaping the inflammatory responses
Analyzing the coding sequence, a total of 646 variants were identified in KTCN patients, including high impact and missense variants (classified as SNVs, deletions, and insertions). In aspects of molecular function, biological process, and their cellular location, genes with identified variants contributed to ECM structural constituent, collagen-containing ECM, and cytoskeleton organization (S2 Figure and S4 Table). What is important, the pathway enrichment analysis comprising all genes with recognized sequence variation revealed the "Antigen presentation" and "Interferon alpha/beta signaling" as the most overrepresented Reactome pathways (FDR < 0.001, S5 Table). In STRING analysis, embracing the same set of genes with the 646 variants, the significant enrichment of 63 clusters, including network of interferon alpha/beta signaling, was recognized (S3 Figure and S6 Table).
Furthermore, analyzing the high-impact and missense variants identified in at least two patients with KTCN, the "Graft-versushost disease" and "Antigen processing and presentation" were found to be the most enriched pathways (Figure 2A and S8 Table).

KCTN-specific non-coding sequence features
Excluding variants recognized in control individuals and variants with MAF>0.1, 15,213 sequence variants in REs in KTCN patients, classified as SNVs, deletions, insertions, indels, and substitutions, were identified. Next, using the FABIAN-variant application that predicts the effects of DNA variant located in RE on transcription factor binding sites (TFBSs), we narrowed down the number of variants and genes for further evaluation to 7,886 and 5,709, respectively. We assessed genes in which close proximity in at least one sequence variant in RE per KTCN patient was recognized as potentially forming or abolishing TFBS. Pathway enrichment analysis of these genes disclosed "Focal adhesion", "Rap1 signaling pathway", and also "Hippo signaling pathway", "Wnt signaling", and "TGF-beta signaling" pathways ( Figure 2B and S9 Table).

Coexistence of coding and non-coding sequence variants in genes influencing the KTCN-specific pathways
Both coding and non-coding sequence variants of KTCN patients were found in the genes/genes' close proximity linked to the revealed KTCN-specific cellular components, namely, "Actin cytoskeleton", "Extracellular matrix", and "Collagen-containing extracellular matrix", and KTCN-specific pathways such as "Focal adhesion", "Hippo signaling pathway", and "Wnt signaling" (S4 Figure and S10 Table). Genes with coding and non-coding sequence variants, attributed to biological pathways from Reactome database, are listed in Table 2, and in detail in S11 Table. To determine if the evaluated genome regions are significantly enriched, genes with both coding and non-coding sequence variants were analyzed and 35 enriched regions were revealed, as presented in Figure 3, with majority of these enriched regions at chromosomes 5q (13 regions), 9q (9 regions), and 16q (5 regions).

Genomic heterogeneity across corneal surface and structural variation
All identified coding-sequence variants were assessed for differences between topographic regions, separately for each adolescent patient with KTCN. No genomic heterogeneity across corneal surface was found. Although identified sequence variants in KIR3DL1 (in two KTCN patients), UMAD1 (in two KTCN patients), and GOLGA6A2 (in one KTCN patient) genes have met quality criteria of reads in the applied bioinformatic algorithms and differentiated the topographic CE regions, the BAM files visualized with the IGV showed a variation in one of 26-43 reads only in these listed genes (the total number of reads for a particular genome site varied), indicating uncertain results. Expanding the same analysis for control samples (but applying the criterion of MAF ≤ 0.1), the additional coding sequence variants in 15 genes (AC018682.2, AC092384.1, AC092490.1, AC117481.1, AC118281.1, ANKRD36BP2, CA15P2, CRACD, HLA-DRB5, KRTAP2-2, LINC00634, NBEAP1, OR2T35, SAMD1, UGT2B25P, and ZNF571) were recognized as potentially differentiating the corneal surface, but none of these variants was confirmed by the BAM files. Representative results of the pathway enrichment analysis. Representative results of the pathway enrichment analysis of (A) genes with the same coding sequence variants identified in more than one KTCN patient, pointing to the involvement of innate and adaptive immune system responses in KTCN pathogenesis and (B) genes located close to the variants in REs. On x axes, −log(FDR) values of enriched pathways are presented, and on y axes, the KEGG pathways are presented. The size of the bubbles indicates the number of genes attributed to the particular pathway and present in our set.  Bold text indicates genes with coding sequence variants identified in more than one KTCN patient; † stands for genes with variants in both and non-coding genome sequence; ‡ stands for genes with the same variant identified in more than one KTCN patient. Variants identified in control samples and variants with MAF≥0.1 in the gnomAD database were excluded (the completed list of genes with sequence variants is presented in S9 Table).
Assessing structural variation, deletions longer than 50 bp were characterized in detail. The length of the deletions and their density in genome are presented in S5 Figure. Most of the deletions were localized in non-coding regions of the genome and in gene introns. We did not observe a difference in the number of deletions comparing patients with KTCN and control individuals, assessing the previously published KTCN-specific loci (S12 Table)

Influence of the presence of RE variants on gene expression
We evaluated the influence of the presence of identified RE variants on expression of genes localized in close proximity, based on data of RNA-seq performed in the same CE samples, and full-thickness corneas (45). No substantial features were recognized in the assessment of the CE. However, 722 of the previously recognized differentially expressed genes (DEGs) in full-thickness KTCN corneas were re-recognized as the gene containing variants in REs. Then, in the pathway enrichment analysis of the rerecognized genes and transcripts, the ECM pathway as the most overrepresented was revealed (Figure 4 and S13 Table).

Discussion
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 The display of genes, in which both coding and non-coding variants unique for adolescent patients with KTCN were identified, across chromosomes. Genes are represented by red dots. The violet lines indicate 35 regions statistically enriched. The genome was scanned with a sliding window in of size 6 Mbp. Within each window, the hypergeometric tests were used to determine if genes of both coding and non-coding variants were significantly overrepresented. The used FDR cutoff was ≤0.05.
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 keratinrelated 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 noncoding sequence variation) suggest that KTCN largely occurs FIGURE 4 Representative results of the pathway enrichment analysis of 722 genes previously recognized as differentially expressed (DEGs) in full-thickness KTCN corneas and re-recognized as containing variants in the REs localized in close proximity. On the x axis, −log(FDR) values of enriched pathways are presented, and on the y axis, the KEGG pathways are presented. The size of bubbles indicated the number of genes attributed to the particular pathway and present in our set. The Reactome database was chosen for the definition of pathways.
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 corneablood 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 gdT cells have been reported in KTCN studies (86). Also in many studies of tear fluid collected in KTCN patients, the higher levels of interleukins (IL1b, IL6, and IL17A), interferons (IFNa/b/g), and other pro-inflammatory factors (TNFa, TGFb1, 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.

Data availability statement
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.

Ethics statement
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.

Author contributions
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.