Comprehensive Transcriptome Analysis of Patients With Keratoconus Highlights the Regulation of Immune Responses and Inflammatory Processes

Keratoconus (KTCN), characterized by the steeper curvature of the cornea and the thinner central corneal thickness, was a degenerative corneal ectasia with ambiguous etiology and mechanism. We aim to reveal underlying pathological mechanisms of KTCN by multi-level transcriptomic, integrative omics analyses. We performed RNA-sequencing on corneal epithelial samples from seven patients and seven control donors, as well as peripheral matched blood samples from three of the patients and three control donors. After RNA extraction, library construction, and sequencing, differentially expressed genes and splicing events were identified, followed by over-representation analysis. In total, 547 differential expressed genes were identified in KTCN corneal epithelial samples, which were mainly enriched in immune responses and inflammatory processes. WGCNA module analysis, the transcriptomic analysis of peripheral blood samples, multiple omics data, and the meta-analysis of GEO samples confirmed the involvement of immune and inflammatory factors. Besides, 190 and 1,163 aberrant splicing events were identified by rMATS combined with CASH methods in corneal epithelial and blood samples with KCTN. In conclusion, this comprehensive transcriptome analysis of KTCN patients based on patients’ tissue and blood samples uncovered a significant association between immune-inflammatory genes and pathways with KCTN, highlighting the contribution of the perturbed immune signaling to the pathogenesis of KCTN. Our study suggested the importance of measures to control inflammation in the treatment of KCTN.

Introduction

Keratoconus (KTCN) is a corneal ectatic disorder characterized by a steeper curvature of the cornea and a thinner central corneal thickness. The global average prevalence of KTCN is 1.38/1,000, with great regional variations, ranging from 0.17/1,000 in Russia to 47.99/1,000 in Saudi Arabia (Hashemi et al., 2020). Current interventions for KTCN, such as contact lens use, corneal ring implantation, and collagen cross-linking, can only slow the progression of the disease. Currently, KTCN cannot be cured by any treatment method. Thus, it is important to understand the fundamental pathogenesis which may be helpful for the development of effective therapeutic interventions.

The exact cause of KTCN is not completely understood. Proteomics studies revealed decreased expression of structural proteins (Chaerkady et al., 2013) and increased expression of degradative enzymes. Pathway analyses indicated the involvement of genes related to extracellular matrix organization (Kabza et al., 2017), transforming growth factor β signaling, WNT pathway (You et al., 2018), and nuclear factor erythroid 2–related factor 2-regulated network (Shinde et al., 2020). In recent years, omics studies improved our understanding of the molecular mechanisms of KTCN (McKay et al., 2016; Kabza et al., 2019; Karolak et al., 2020).

Although KTCN was previously considered a non-inflammatory disease, there is increasing evidence supporting the involvement of inflammation and immune responses in KTCN. Various studies found inflammatory mediators in the corneal tissue and tears of patients with KTCN, including interleukin (IL)-1β (Sorkhabi et al., 2015), IL-6, tumor necrosis factor (TNF)-α (Lema et al., 2009), matrix metallopeptidase-9 (Balasubramanian et al., 2012a), and nuclear factor kappa B. Molecules involved in the immune response, such as lactoferrin (LTF), zinc-α2-glycoprotein, and immunoglobulin kappa chain, are downregulated in the tears of patients with KTCN. The serum expression levels of Toll-like receptors 2 and 4 in monocytes and neutrophils were significantly higher in patients with KTCN compared to control subjects (Sobrino et al., 2017). The latest literature shows that KTCN is positively correlated with a variety of immune-mediated diseases, such as Hashimoto’s thyroiditis and inflammatory skin conditions (Claessens et al., 2021).

Our study aimed to determine the underlying pathological mechanisms of KTCN by multi-level transcriptomic study and integrative omics analyses. To the best of our knowledge, this is the first comparative transcriptomic study of corneal epithelial and blood samples from patients with KTCN and control donors. Through bioinformatics analyses, real-time polymerase chain reaction (RT-PCR) validation, and co-analysis with published omics data and GEO data, we revealed the involvement of immune and inflammatory factors in KTCN pathogenesis. In addition, we analyzed alternative splicing events which may contribute to the pathophysiology of KTCN and made a prediction for potential drugs associated with KTCN.

Materials and Methods

Samples

This study conformed to the Declaration of Helsinki and was approved by the Research Ethics Committee of Tianjin Eye Hospital (KY202107). All participants provided written informed consent. Patients who were diagnosed with KTCN and underwent de-epithelized collagen cross-linking at Tianjin Eye Hospital from May 2020 to May 2021 were enrolled. Controls were sex-matched patients undergoing laser-assisted subepithelial keratomileusis surgery for mild myopia without ocular and other systemic diseases. All participants underwent a complete medical history assessment, general examination, and specialist examination. A total of seven KTCN corneal epithelial tissues, three KTCN blood samples, seven control corneal epithelial tissues, and three control blood samples were collected. The removed epithelium was immediately submerged in RNALater (Qiagen, Germany) or immediately subjected to the workflow of RNA extraction. Peripheral blood was collected in ethylenediaminetetraacetic acid-containing tubes, and peripheral blood mononuclear cells (PBMCs) were isolated using red blood cell lysis buffer (TianGen, Beijing).

RNA Isolation, cDNA Library Construction, and Sequencing

RNA simple Total RNA Kit (TianGen, Beijing) was used to extract RNA from tissues and PBMCs. Briefly, RNALater was carefully removed, and 1 ml of TRIzol was added to each tube. After phase separation by adding 200 μl chloroform and centrifugation, the aqueous phase was removed to a new tube. Moreover, 0.5-fold volume of anhydrous ethanol was added, mixed well, and transferred to the adsorption column. The column was deproteinized and washed with buffer PD, followed by two more washes with buffer PW. The column was dried at 12,000 × g for 2 min at room temperature, and 30 μl of RNase-free water was added to elute the RNA. RNA quantity and integrity were detected by NanoDrop2000 and Agilent2100 BioAnalyzer, respectively.

The first strand of cDNA was synthesized using M-MuLV reverse transcriptase, random oligonucleotide, and cDNA fragment. After degrading the RNA strand with RNase H, the DNA polymerase I system was used to synthesize the second strand. The purified double-strand cDNA underwent cDNA terminal repair, poly-A tail addition, connector connection, length screening, and PCR amplification and purification and finally obtained the library. The obtained cDNA fragments were evaluated using the Agilent BioAnalyzer 2100 system, and sequencing was performed on the Illumina platform.

Transcriptomics Data Analysis

Clean reads were obtained by removing low-quality reads and ligand from the original sequence. The Q20/Q30 values and GC content were calculated to ensure quality. Clean reads were accurately aligned with the human reference genome (hg19) using HISAT2 software (Mortazavi et al., 2008), and the FeatureCounts in SubRead software (Liao et al., 2014) was used to calculate the read number. Read numbers lower than 10, unpaired reads, and reads in multiple regions of the genome were filtered out. Sequencing depth and gene length were corrected (Bray et al., 2016), and the expression values (FPKM) of all genes in each sample were calculated. To determine the associations between samples, we calculated the square of the Pearson correlation coefficient (R²) of the samples within groups and performed principal component analysis (PCA) using the calculation square of linear algebra based on FPKM.

Expressed genes were defined as significant with the DESeq2 cutoff of p value < 0.05 and |log 2fold change |> 1. Meanwhile, clustering analysis and visualization were performed for differentially expressed genes (DEGs) using the R package “vocano.” Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the R package “clusterprofile” and web-based gene set analysis toolkit WebGestalt (Liao et al., 2019) based on DEGs.

To better interpret the expression data, the weighted gene co-expression network analysis (WGCNA) algorithm in integrated Differential Expression and Pathway analysis (Ge et al., 2018) was applied to perform the cluster analysis of samples at the same time. WGCNA constructed a clustering tree based on the correlation of gene expression and then divided it into modules. The target genes and gene networks could be identified by exploring the correlation between modules and phenotypes. The most variable genes to include were set as 1,000, and the minimum module size was set to 20. In addition, the corresponding module information was extracted to perform GO enrichment analysis and visualized with a bubble chart.

Quantitative Polymerase Chain Reaction

7 pairs of corneal epithelial samples from patients and control donors were used for validation in RT-PCR. SYBR Green qPCR response system was 20 µl, including 2 × PerfectStart Green qPCR SuperMix (10 µl), cDNA (1 µl), and forward and reverse primer (0.4 µl), filled with water to total 20 µl. Furthermore, 45 cycles were set as follows: denaturation at 95°C for 5 s, annealing at 60°C for 30 s, and extension at 72°C for 10 s. The qPCR primer sequences for the detected genes are listed in Supplementary Table S10. 3 technical replicates were performed for each biological reaction and Graphpad Prism 6.0 was used for statistical analysis.

Comparison to Keratoconus Candidate Genes in the Literature

Genomic studies of KTCN by searching PubMed and GWAS repositories were included if they identified candidate gene loci through linkage study, genome-wide association study, and whole-exome sequencing study in patients with KTCN. Candidate genes associated with KTCN published by the original author were collected, and this process identified 46 candidate genes from 36 articles combined with the database. Proteomics studies of KTCN were included if they identified protein expression changes in patients with KTCN. An extensive search of PubMed collected 12 proteomics studies that met the requirements (Supplementary Table S9). Using the original author’s statistical criteria, 1844 differentially expressed proteins were identified. We then compared our transcriptome data with candidate gene loci and a proteomics gene list and visualized them using the GO plot package.

Identification of Differentially Splicing Events

Alternative splicing (AS) events were detected using two parallel approaches, rMATS (Shen et al., 2014) and CASH (Wu et al., 2018). The rMATS method quantifies the expression of AS events using a hierarchical model to simultaneously account for sampling uncertainty and variability. It calculates the p-value using a likelihood-ratio test to represent the difference in the inclusion level between groups of samples. The events with false discovery rate (FDR) < 0.05 were considered as significantly different between samples. Five types of splicing events were examined, including skipped exon (SE), alternative 3′ splice site (A3SS), alternative 5′ splice site (A5SS), mutually exclusive exons (MXE), and retained intron (RI). The CASH method directly self-constructs AS sites using the SpliceCons module from RNA-Seq data to explore novel AS events. Accordingly, SpliceDiff module merges the expression of AS inclusion/exclusion events and the expression of AS exons to infer differential levels of AS events between samples. The significant splice events (FDR < 0.05) detected in the patients were classified as cassette exon (i.e., skipped exon), multi-cassette exons, A3SS, A5SS, alternative start exon (AltStart), and alternative end exon (AltEnd), MXE and intron retention (IR i.e., retained intron). Differentially spliced genes were detected from differentially expressed alternative splicing events and were further applied for GO annotation.

RNAseq Meta-Analysis

Two datasets (GSE77938, GSE112155) from the Gene Expression Omnibus (GEO) database were included in the meta-analysis. GSE77938 comprises 132 RNA-Seq profiles (25 cases with a later time-point RNA-seq data file being removed due to conversion failure, 25 controls); GSE112155 contains 20 RNA-Seq profiles (10 cases, 10 controls). RNA-seq data of the three databases were included to meta-analysis. Hisat2 (Kim et al., 2015) was used to assign RNA-seq data and FeatureCounts (Liao Y et al., 2014) was utilized to assemble KCTN RNA-seq data with control RNA-seq data to gain raw counts. Raw counts were then processed by Bioconductor packages in R v4.0.5 (www.r-project.org). We removed batch effect and other unwanted variation for RNA sequencing data by sva package. Normalization, DEGs analysis and annotation were conducted by preprocessCore package, DESeq2 package and org.Hs.eg.db respectively. Benjamini and Hochberg method control the False Discovery Rate (FDR) for multiple-testing adjustment. DEGs were defined as FDR < 0.05, and |log (fold-change) | >2. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted by R package “clusterprofile” and Webgestalt based on DEGs.

Protein-Protein Interaction Analysis and Protein-Drug Prediction for Hub Genes

DEGs were imported into Network Analyst (Zhou et al., 2019) (https://www.networkanalyst.ca/) to construct a PPI network, PPI data were derived from STRING interactome (Szklarczyk et al., 2019) with a high confidence score of ≥ 900 and requirement for experimental evidence. Besides, the PPI network was visualized by Cytoscape (Cline et al., 2007), and hub genes were predicted via CytoHubba as previously described (Sun et al., 2021). Protein-chemical interaction analysis was performed based on data from the Comparative Toxicogenomics Database (CTD) (Mattingly et al., 2006) (downloaded on November 2016), a database identifies interactions between chemicals and genes, and promotes understanding about the effects of environmental chemicals on human health. Drugs and chemicals were ranked according to degree and betweenness.

Results

Assessments of Study Patients

We recruited patients with KTCN of East Asian ancestry from Tianjin, China (aged 12–19 years old). As controls for the KTCN samples, we used samples of mild myopia patients (aged 18–24 years old) from the same geographical and ethnic populations as cases. The mean age of patients with KTCN (16.00 ± 2.27 years) was comparable to that of the normal patients (20.14 ± 1.81 years). The corneal topography images showed significant differences in parameters between the two studied groups in anterior and posterior keratometry at flat axis (K1), steep axis (K2), and thinnest corneal thickness (Kmean) (Supplementary Table S1, Supplementary Figure S1).

The average RNA concentration and RNA integrity number value among tissue and blood samples were 265 ng/μl and 6.75 and 360.5 ng/μl and 6.47, respectively. Clean reads were obtained after raw data filtering, sequencing error rate assessment, and GC content distribution evaluation. The mean values of clean reads Q20 and Q30 were 97.96 and 94.26%, respectively. The sequencing error rates of the single bases were all less than 3% (Supplementary Table S2).

Gene Expression Distribution and the Correlation Between Samples

After calculating the expression values FPKM of all the genes in each sample, the distribution of gene expression levels in different samples is shown in the box diagram (Supplementary Figures S2A, S2B). Across all samples, 12,187–12,665 unique protein-coding genes per sample with FPKM ≥ 1 were detected, and more than 351 genes were highly expressed (FPKM ≥ 60). The correlation analysis showed a high degree of Pearson correlation coefficient (R² = 0.795–0.969) between samples from patients with KTCN and control patients (Supplementary Figures S2C, S2D). PCA showed close clustering by tissue type (Figure 1).

FIGURE 1

FIGURE 1. Principal component analysis (PCA) for corneal epithelial samples and blood samples. Red dots represent KTCN corneal epithelial samples, green triangles represent control corneal epithelial samples, blue squares represent KTCN blood samples, and purple crosses represent control blood samples.

Dysregulated Pathways Between Keratoconus and Control Corneal Epithelium

Based on the cutoff mentioned in the method, 547 genes were identified as DEGs (308 upregulated and 239 downregulated) between KTCN and control corneal epithelial samples (Figures 2A,B; Supplementary Table S3). The expression levels of secreted frizzled-related protein 1 (SFRP1), TPPP3, KCNB2, RASSSF2, CD97, BMP3, SLC16A14, HLA-DPA1, LIX1, FGF12, GRIK5, CRIP2, FGF9, PRR11, and CDCA2 ranked top 15 significantly upregulated genes in the KTCN samples, whereas LAMP3, SEPP1, CTSH, CYP4F12, TGM5, PILRB, KCNJ2, HPGD, ADAM23, ZNF469, SLC1A3, FMO1, PTN, and FAM21B were the top 15 downregulated genes (Table 1).

FIGURE 2

FIGURE 2. Cluster analysis of differentially expressed genes in corneal epithelial samples and RNA level expression validation. (A) Hierarchical clustering analysis of differentially expressed genes by a heatmap plot. Higher enrichment is colored in red, the lower ones are colored in blue. (B) Differentially expressed genes are shown by a volcano plot in corneal epithelial samples. Red dots indicate upregulated genes, and blue ones represent downregulated genes. (C) Eight genes were validated by qPCR. CON_T, control corneal epithelial samples, KC_T, KTCN corneal epithelial samples. (*p < 0.01).

TABLE 1

TABLE 1. The top dysregulated differentially expressed genes in corneal epithelial samples of patients with keratoconus.

The DEGs were validated by randomly selecting eight genes from the top 30 identified DEGs. Among the tested DEGs, consistent with our RNA-Seq findings, SFRP1, FGF9, and BMP3were significantly upregulated in the KTCN epithelium. In contrast, SLC1A3, HPDG, DCN, CTSH, and CYP4F1 RNA expressions were downregulated in the KTCN epithelium (Figure 2C).

The GO annotation of the 547 DEGs revealed that the biological processes of many DEGs were related to leukocyte cell-cell adhesion (such as HLA-DPA1 and S100A8), regulation of lymphocyte activation (such as SFRP1and CD247), and adaptive immune responses (such as CTSH and LAMP3) (Figure 3A; Table 2; Supplementary Table S4). Meanwhile, immune- and inflammatory-related KEGG pathways, including Th17 cell differentiation, Th1 and Th2 cell differentiation, and IL-17 signaling pathway, were observed in the DEGs of the samples (Figure 3B; Table 3; Supplementary Table S5). The enrichment of inflammatory immune diseases, such as rheumatoid arthritis, bowel disease, and autoimmune thyroid disease, was also observed, some of which have been reported to be a comorbidity accompanied by KTCN previously (Nemet et al., 2010; Tréchot et al., 2015).

FIGURE 3

FIGURE 3. Enrichment analysis of differentially expressed genes in corneal epithelial samples. (A) Gene Ontology analysis of differentially expressed genes. Biological process (BP) is colored in red. Cellular component (CC)is colored in green. Molecular function (MF) is colored in blue. (B) Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis based on differentially expressed genes. CON, control samples, KC, KTCN samples.

TABLE 2

TABLE 2. The significantly biological process of Gene Ontology analysis for 547 differentially expressed genes.

TABLE 3

TABLE 3. The significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways for 547 differentially expressed genes.

Using WGCNA, we identified seven co-expression modules that were constructed from the 896 genes. These co-expression modules were constructed and are shown in different colors (Figure 4A). The number of genes in these modules ranged from 23 to 424 (Supplementary Table S6). We found that the largest turquoise module was significantly correlated with the immune system processes. GO enrichment analysis was performed on the genes in the largest turquoise module. Genes in the turquoise module were mainly enriched in GO:0002283 (neutrophil activation involved in immune response), GO:0043312 (neutrophil degranulation), and GO:0002429 (immune response-activating cell surface receptor signaling pathway) (Figure 4B).

FIGURE 4

FIGURE 4. Identification of modules and functional annotation analysis for the module genes in corneal epithelial samples. (A) Module detection for the expressed genes of corneal epithelial samples. (B) Gene Ontology analysis for the turquoise module genes.

Expression Changes in the Blood of Patients With Keratoconus Emphasize Immune Regulation

We next determined whether networks of KTCN whole blood have biological processes similar to those of the corneal epithelium. RNA-Seq and enrichment annotations were performed for peripheral blood using the same method as that used for the tissue. Of the 28,080 genes analyzed, 509 genes were upregulated and 622 were downregulated in blood samples from patients with KTCN versus controls (Figures 5A,B). We found a significant excess of biological processes to be significantly enriched in the humoral immune response (padj < 0.05). The changes were prominent in genes related to HLA-DQB1, LTF, and PAX5. (Figures 5C,D; Supplementary Tables S3–S5).

FIGURE 5

FIGURE 5. Hierarchical cluster analysis of differentially expressed genes in blood samples. (A) Hierarchical clustering analysis of differentially expressed genes by a heatmap plot. Higher enrichment is colored in red, the lower ones are colored in blue. (B) Differentially expressed genes were showed by a volcano plot in blood samples. Red dots represent upregulated genes, and blue ones indicate downregulated genes. (C) Gene Ontology analysis of differentially expressed genes. (D) Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis based on differentially expressed genes.

We then explored molecules that act in the same direction in the whole blood and corneal epithelium. The overlap in dysregulated genes was not significantly larger than expected. However, it still has important implications, allowing us to determine common key factors between the two systems. Genes, such as TUBB1 and BANK1, that significantly (p < 0.05) enriched for dysregulation in both the corneal epithelium and whole blood are listed in Supplementary Table S7.

Co-Analysis of Previous Genomics and Proteomics Study of Keratoconus Reinforced the Conjecture

5 and 10 DEGs overlapped with the genetic candidate loci and the proteomics gene list, respectively (Figures 6A,B). We noted the inflammatory factor IL1B in the limited overlapping genes between genetic candidate loci and proteomics gene list. In addition, the GO clusters enriched in the list of overlapping genes from proteomics studies and our transcriptome study displayed evident immune and inflammatory involvement, which included acute inflammatory response (S100A8, C1QA, SERPINA1), humoral immune response (LYZ, S100A8, C1QA), and neutrophil-mediated immunity (LYZ, S100A8, SERPINA1).

FIGURE 6

FIGURE 6. Overlap gene and Gene Ontology (GO) cluster analysis. (A) The list of candidate gene loci that were also differentially expressed in our transcriptome study. (B) The enriched GO clusters for the dysregulated genes in our transcriptome and previous proteomics studies.

Alternative Splicing Events in Patients With Keratoconus

We also sought to identify transcripts that were differentially spliced between KTCN and control patients 67 candidate aberrant splicing events were identified with a threshold of FDR <0.05 by rMATS in cornea samples with KCTN. CASH detects more aberrant splicing events up to 132 compared to rMATS. On the other hand, CASH and rMATS method showed 1,076 and 218 differential splicing events in blood samples, among which 35 events were identified by both methods (Figure 7A).

FIGURE 7

FIGURE 7. An overview of alternative splicing events in corneal epithelial and blood samples. (A). The significantly differentially spliced events detected by two methods. (B) Pie charts showing the types and relative proportion of alternative splicing events in corneal epithelial and blood samples. (C) Venn diagrams depicting overlaps and differences between the events in corneal and blood samples. (D) The top 10 enriched GO biological processes among differentially spliced genes in blood samples.

The current study revealed 1,163 and 190 alternative splicing events in blood and cornea. The alternative splicing events in corneal samples included SE (38％), IR (22％), AltEnd (12％), AltStart (10％), MXE (9％), A5SS (5％) and A3SS (4％). In blood sample, the corresponding sequence was: SE (34%) > IR > AltEnd> AltStart > MXE> A3SS > A5SS (Figure 7B). 3 /71 SE, 4 / 42 IR, and 1 / 8 A3SS events overlapped in corneal and blood samples. The changes of A5SS, AltStart, AltEnd, and MXE were totally different between corneal epithelium and blood (Figure 7C).

The overrepresentation analysis for affected splicing events in corneal epithelial samples revealed several processes, but with high p-value. A few processes for aberrant splicing events in blood samples were captured in our functional annotation, including neutrophil activation, neutrophil activation involved in immune response, granulocyte activation, and neutrophil mediated immunity (Figure 7D).

Mata-Analysis of GEO RNA-Seq Data of Keratoconus Confirms Our Suspicion

Through meta-analysis and enrichment analysis performed on the two previously reported datasets, again, Th17 cell differentiation, T-cell receptor signaling pathway, and IL-17 signaling were found to be enriched in DEGs of KCTN. This confirmed that the change of immune gene expression and the disorder of immune-related signaling pathways may mediate the molecular mechanism of KCTN. We also observed the enrichment of the “ECM receptor interaction” pathway, which is consistent with previous reports that extracellular matrix degradation may be involved in the pathogenesis of KCTN (Supplementary Tables S3–S5; Supplementary Figure S3).

Potential Therapeutic Targets for Keratoconus

On analyzing the DEGs using the GDA interface of NetworkAnalyst to acquire subnetworks (at least 3 nodes), there were 29 subnetworks created of which only 1 (subnetwork 1) had a seed number above 3. The subnetwork1 was identified with 806 nodes and 1,092 connection edges. We used Cytoscape software to visualize the subnetwork1 (Figure 8A). The top 15 DEGs (JUN, RAC2, FOS, GNA14, JAK3, RET, EGR1, HSPA6, TUBB1, ETS1, GLI2, LCP2, FZD10, SFRP1, and AXIN2) were regarded as hub genes by node degrees. Among them, JUN was the most highly ranked gene (degree = 129; betweenness = 80856.21).

FIGURE 8

FIGURE 8. Protein-protein and protein-chemical interaction analysis for differentially expressed genes. (A) The protein-protein interaction (PPI) network and the top 15 hub genes in the network ranked by the Degree method. Nodes indicate the number of predicted gene interactions. Red squares are the hub genes. (B) The top 10 drugs associated with KTCN DEGs in the protein-chemical interactions network. (C) The top 10 drugs associated with hub genes. Overlapped drugs in two figures were marked in red.

The protein–chemical interaction analysis of DEGs help to explore potential therapeutic drugs in treating KTCN. A single network was identified with 1887 nodes and 8,116 edges. 1,470 chemicals were identified to interact with DEGs. The top ten compounds associated with KTCN DEGs were predicted in Figure 8B. We also found that valproic acid, SB431542, and dorsomorphin were drugs that could interact with most hub genes. The top ten compounds which interact with hub genes were listed in Figure 8C.

Discussion

Initial support for immune evidence in KTCN was obtained in an earlier proteomics study of KTCN tear samples in which immunoglobulin kappa chain (Lema et al., 2010) and polymeric immunoglobulin receptor (Balasubramanian et al., 2013) were decreased. Additionally, Nielsen et al., detected an increase in major histocompatibility complex class II, DR alpha (HLA-DRA1), using an expression microarray (Nielsen et al., 2003). In a recent study of two KTCN patient groups, using the canonical pathways method, researchers detected a large percentage of the RA pathway in the Middle Eastern ancestry group, emphasizing a dysregulation in inflammation and immune signals in the cornea (Shinde et al., 2020).

In the current study, we found evidence of the involvement of immune-inflammatory substances in the corneal epithelial transcriptome. Except for HLA-DPA1 and PILRB, SFRP1 was upregulated. SFRP1, which is an antagonist of Wnt signaling, is dysregulated in many studies concerning KTCN (Sutton et al., 2010; You et al., 2013a; You et al., 2013b; Iqbal et al., 2013; Khaled et al., 2018). SFRP1 participates in cell proliferation, migration, and differentiation. It might also serve as an anti-inflammatory factor by modulating the balance between pro-inflammatory and anti-inflammatory cytokines (Barandon et al., 2011). SFRP1 was reported to function in Th17 cell differentiation and was significantly correlated with IL-17 levels in patients with rheumatoid arthritis (Lee et al., 2012).

The internal consistency of our data between the blood and corneal epithelial transcriptomes allowed us to verify immune responses as core pathogenic changes in KTCN. The changes were prominent in humoral immune genes, such as HLA-DQB1 and LTF. Taking the LTF as an example in our study, downregulation of LTF was identified in the blood of patients with KTCN (Lema et al., 2010; Balasubramanian et al., 2012b; Priyadarsini et al., 2014). LTF, an iron-binding protein, is considered to play an important role in modifying innate and adaptive immune responses (Moreno-Expósito et al., 2018). LTF can act on antigen-presenting cells (Puddu et al., 2009) to mature dendritic cells, activate macrophages, accelerate antigen processing, and increase NK cell cytotoxicity. By acting on the maturation process of T lymphocytes, LTF can induce CD4 expression and direct differentiation of immature T lymphocytes toward CD4 + T lymphocytes, thus regulating the functional ability of T lymphocytes (Siqueiros-Cendón et al., 2014). LTF can promote TH1 response and inhibit the TH2 response, thereby activating cellular response and reducing the release of inflammatory factors. In addition, various clinical trials have found that LTF may exert anti-inflammatory effects by controlling TNF-α (Drago-Serrano et al., 2017). Knockout of LTF mice (LTF-KO) has been shown to be deficient in hematopoietic and immune systems. Neutrophil maturation, migration, phagocytosis, granule release, and antimicrobial response to bacterial challenge are affected in LTF-KO mice (Ward et al., 2008).

Although many studies have examined the changes in KTCN corneas, few have analyzed the blood samples of these patients, especially by transcriptome analysis. Our study, for the first time, quantified the gene expression profile in the whole blood of patients with KTCN, which allowed us to reexamine the mechanism of KTCN from the perspective of system theory. In addition, dysregulated genes in the blood may yield early biomarkers for KTCN. Most importantly, we have now found evidence to propose an “immune-inflammatory hypothesis” in KTCN. Using multilevel transcriptome data, we hypothesize that KTCN is a systemic disease with immune component and inflammatory basis.

A few limitations of this study must be considered. First, our experiment was restricted by the number of samples, and more samples from patients and controls are needed to increase the statistical power. Second, the corneal epithelium was used rather than the whole cornea button because of the limitation of the surgical procedure. As the most prominent layer of lesion, the biological data for the stromal layer of keratoconus are undoubtedly important. Further studies to perform on the corneal stroma lever if the surgical materials allowed and validate the protein levels of the altered gene expression identified in our study should be conducted. Based on the current results, it is also necessary to perform the next targeted molecular function experiments and discover potential biomarkers for KTCN detection and treatment.

In summary, we performed a comprehensive analysis of the expression profile in the corneal epithelium and peripheral blood of patients with KTCN and control donors to understand the pathogenesis of KTCN. A total of 547 DEGs were identified in the tissue expression profiles. Functional enrichment analysis showed that these DEGs were mainly related to leukocyte cell-cell adhesion, regulation of lymphocyte activation, and adaptive immune response. Concordantly, DEGs were significantly enriched in Th17 cell differentiation, Th1 and Th2 cell differentiation, and IL-17 signaling pathways. Meanwhile, through transcriptomics analysis of peripheral blood samples and co-analysis of previously reported omics studies, enrichment of immune-related biological processes was also observed. For the first time, RNA-Seq was used to compare gene expression in both the corneal epithelial and blood samples. Our results indicated the potential involvement of immune pathways in KTCN pathogenesis.

Data Availability Statement

RNA seq data were deposited in the NCBI Sequence Read Archive (SRA) under accession number: PRJNA799648.

Ethics Statement

The studies involving human participants were reviewed and approved by the Research Ethics Committee of Tianjin Eye Hospital. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

XS and YW conceived and designed the study. YW, LC, and LZ. performed the surgery. XS, YD and MS collected samples and clinical file information, and XS and HZ analyzed the data and computer figuring. XS and YW wrote the manuscript. All authors read and approved the final manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (Nos. 81873684 and 11701294), Tianjin Science and Technology Research Project (16JCYBJC25800), Key Projects of Tianjin Health Industry (15KG120), and Science and Technology Foundation Projects of Tianjin Eye Hospital (YKZD 2002).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Professor LC for supporting the keratoconus samples.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.782709/full#supplementary-material

References

Balasubramanian, S. A., Wasinger, V. C., Pye, D. C., and Willcox, M. D. (2013). Preliminary Identification of Differentially Expressed Tear Proteins in Keratoconus. Mol. Vis. 19, 2124–2134.

Google Scholar

Balasubramanian, S. A., Mohan, S., Pye, D. C., and Willcox, M. D. P. (2012a). Proteases, Proteolysis and Inflammatory Molecules in the Tears of People with Keratoconus. Acta Ophthalmol. 90, e303–e309. doi:10.1111/j.1755-3768.2011.02369.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Balasubramanian, S. A., Pye, D. C., and Willcox, M. D. P. (2012b). Levels of Lactoferrin, Secretory IgA and Serum Albumin in the Tear Film of People with Keratoconus. Exp. Eye Res. 96, 132–137. doi:10.1016/j.exer.2011.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Barandon, L., Casassus, F., Leroux, L., Moreau, C., Allières, C., Lamazière, J.-M. D., et al. (2011). Secreted Frizzled-Related Protein-1 Improves Postinfarction Scar Formation through a Modulation of Inflammatory Response. Atvb 31, e80. doi:10.1161/ATVBAHA.111.232280

CrossRef Full Text | Google Scholar

Bray, N. L., Pimentel, H., Melsted, P., and Pachter, L. (2016). Near-optimal Probabilistic RNA-Seq Quantification. Nat. Biotechnol. 34, 525–527. doi:10.1038/nbt.3519

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaerkady, R., Shao, H., Scott, S.-G., Pandey, A., Jun, A. S., and Chakravarti, S. (2013). The Keratoconus Corneal Proteome: Loss of Epithelial Integrity and Stromal Degeneration. J. Proteomics 87, 122–131. doi:10.1016/j.jprot.2013.05.023

CrossRef Full Text | Google Scholar

Claessens, J. L. J., Godefrooij, D. A., Vink, G., Frank, L. E., and Wisse, R. P. L. (2021). Nationwide Epidemiological Approach to Identify Associations between Keratoconus and Immune-Mediated Diseases. Br. J. Ophthalmol., doi:10.1136/bjophthalmol-2021-318804

CrossRef Full Text | Google Scholar

Cline, M. S., Smoot, M., Cerami, E., Kuchinsky, A., Landys, N., Workman, C., et al. (2007). Integration of Biological Networks and Gene Expression Data Using Cytoscape. Nat. Protoc. 2, 2366–2382. doi:10.1038/nprot.2007.324

PubMed Abstract | CrossRef Full Text | Google Scholar

Drago-Serrano, M., Campos-Rodríguez, R., Carrero, J., and de la Garza, M. (2017). Lactoferrin: Balancing Ups and downs of Inflammation Due to Microbial Infections. Ijms 18, 501. doi:10.3390/ijms18030501

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, S. X., Son, E. W., and Yao, R. (2018). iDEP: an Integrated Web Application for Differential Expression and Pathway Analysis of RNA-Seq Data. BMC Bioinformatics 19, 534. doi:10.1186/s12859-018-2486-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hashemi, H., Heydarian, S., Hooshmand, E., Saatchi, M., Yekta, A., Aghamirsalim, M., et al. (2020). The Prevalence and Risk Factors for Keratoconus: A Systematic Review and Meta-Analysis. Cornea 39, 263–270. doi:10.1097/ICO.0000000000002150

PubMed Abstract | CrossRef Full Text | Google Scholar

Iqbal, O., Fisher, G., Vira, S., Syed, D., Sadeghi, N., Freeman, D., et al. (2013). Increased Expression of Secreted Frizzled-Related Protein-1 and Microtubule-Associated Protein Light Chain 3 in Keratoconus. Cornea 32, 702–707. doi:10.1097/ICO.0b013e318282987a

PubMed Abstract | CrossRef Full Text | Google Scholar

Kabza, M., Karolak, J. A., Rydzanicz, M., Szcześniak, M. W., Nowak, D. M., Ginter-Matuszewska, B., et al. (2017). Collagen Synthesis Disruption and Downregulation of Core Elements of TGF-β, Hippo, and Wnt Pathways in Keratoconus Corneas. Eur. J. Hum. Genet. 25, 582–590. doi:10.1038/ejhg.2017.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kabza, M., Karolak, J. A., Rydzanicz, M., Udziela, M., Gasperowicz, P., Ploski, R., et al. (2019). Multiple Differentially Methylated Regions Specific to Keratoconus Explain Known Keratoconus Linkage Loci. Invest. Ophthalmol. Vis. Sci. 60, 1501–1509. doi:10.1167/iovs.18-25916

PubMed Abstract | CrossRef Full Text | Google Scholar

Karolak, J. A., Gambin, T., Rydzanicz, M., Polakowski, P., Ploski, R., Szaflik, J. P., et al. (2020). Accumulation of Sequence Variants in Genes of Wnt Signaling and Focal Adhesion Pathways in Human Corneas Further Explains Their Involvement in Keratoconus. PeerJ 8, e8982. doi:10.7717/peerj.8982

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaled, M. L., Bykhovskaya, Y., Yablonski, S. E. R., Li, H., Drewry, M. D., Aboobakar, I. F., et al. (2018). Differential Expression of Coding and Long Noncoding RNAs in Keratoconus-Affected Corneas. Invest. Ophthalmol. Vis. Sci. 59, 2717–2728. doi:10.1167/iovs.18-24267

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y.-S., Lee, K.-A., Yoon, H.-B., Yoo, S.-A., Park, Y. W., Chung, Y., et al. (2012). The Wnt Inhibitor Secreted Frizzled-Related Protein 1 (sFRP1) Promotes Human Th17 Differentiation. Eur. J. Immunol. 42, 2564–2573. doi:10.1002/eji.201242445

PubMed Abstract | CrossRef Full Text | Google Scholar

Lema, I., Brea, D., Rodríguez-González, R., Díez-Feijoo, E., and Sobrino, T. (2010). Proteomic Analysis of the Tear Film in Patients with Keratoconus. Mol. Vis. 16, 2055–2061.

PubMed Abstract | Google Scholar

Lema, I., Sobrino, T., Duran, J. A., Brea, D., and Diez-Feijoo, E. (2009). Subclinical Keratoconus and Inflammatory Molecules from Tears. Br. J. Ophthalmol. 93, 820–824. doi:10.1136/bjo.2008.144253

CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an Efficient General Purpose Program for Assigning Sequence Reads to Genomic Features. Bioinformatics 30, 923–930. doi:10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Wang, J., Jaehnig, E. J., Shi, Z., and Zhang, B. (2019). WebGestalt 2019: Gene Set Analysis Toolkit with Revamped Uis and APIs. Nucleic Acids Res. 47, W199–W205. doi:10.1093/nar/gkz401

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattingly, C. J., Rosenstein, M. C., Colby, G. T., Forrest Jr, J. N., and Boyer, J. L. (2006). The Comparative Toxicogenomics Database (CTD): a Resource for Comparative Toxicological Studies. J. Exp. Zool. 305A, 689–692. doi:10.1002/jez.a.307

CrossRef Full Text | Google Scholar

McKay, T. B., Hjortdal, J., Sejersen, H., Asara, J. M., Wu, J., and Karamichos, D. (2016). Endocrine and Metabolic Pathways Linked to Keratoconus: Implications for the Role of Hormones in the Stromal Microenvironment. Sci. Rep. 6, 25534. doi:10.1038/srep25534

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreno-Expósito, L., Illescas-Montes, R., Melguizo-Rodríguez, L., Ruiz, C., Ramos-Torrecillas, J., and de Luna-Bertos, E. (2018). Multifunctional Capacity and Therapeutic Potential of Lactoferrin. Life Sci. 195, 61–64. doi:10.1016/j.lfs.2018.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L., and Wold, B. (2008). Mapping and Quantifying Mammalian Transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. doi:10.1038/nmeth.1226

PubMed Abstract | CrossRef Full Text | Google Scholar

Nemet, A. Y., Vinker, S., Bahar, I., and Kaiserman, I. (2010). The Association of Keratoconus with Immune Disorders. Cornea 29, 1261–1264. doi:10.1097/ICO.0b013e3181cb410b

PubMed Abstract | CrossRef Full Text | Google Scholar

Nielsen, K., Birkenkamp-Demtro¨der, K., Ehlers, N., and Orntoft, T. F. (2003). Identification of Differentially Expressed Genes in Keratoconus Epithelium Analyzed on Microarrays. Invest. Ophthalmol. Vis. Sci. 44, 2466–2476. doi:10.1167/iovs.02-0671

PubMed Abstract | CrossRef Full Text | Google Scholar

Priyadarsini, S., Hjortdal, J., Sarker-Nag, A., Sejersen, H., Asara, J. M., and Karamichos, D. (2014). Gross Cystic Disease Fluid Protein-15/prolactin-Inducible Protein as a Biomarker for Keratoconus Disease. PLOS ONE 9, e113310. doi:10.1371/journal.pone.0113310

PubMed Abstract | CrossRef Full Text | Google Scholar

Puddu, P., Valenti, P., and Gessani, S. (2009). Immunomodulatory Effects of Lactoferrin on Antigen Presenting Cells. Biochimie 91, 11–18. doi:10.1016/j.biochi.2008.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, S., Park, J. W., Lu, Z.-x., Lin, L., Henry, M. D., Wu, Y. N., et al. (2014). rMATS: Robust and Flexible Detection of Differential Alternative Splicing from Replicate RNA-Seq Data. Proc. Natl. Acad. Sci. USA 111, E5593–E5601. doi:10.1073/pnas.1419161111

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinde, V., Hu, N., Mahale, A., Maiti, G., Daoud, Y., Eberhart, C. G., et al. (2020). RNA Sequencing of Corneas from Two Keratoconus Patient Groups Identifies Potential Biomarkers and Decreased NRF2-Antioxidant Responses. Sci. Rep. 10, 9907. doi:10.1038/s41598-020-66735-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Siqueiros-Cendón, T., Arévalo-Gallegos, S., Iglesias-Figueroa, B. F., García-Montoya, I. A., Salazar-Martínez, J., and Rascón-Cruz, Q. (2014). Immunomodulatory Effects of Lactoferrin. Acta Pharmacol. Sin. 35, 557–566. doi:10.1038/aps.2013.200

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobrino, T., Regueiro, U., Malfeito, M., Vieites-Prado, A., Pérez-Mato, M., Campos, F., et al. (2017). Higher Expression of Toll-like Receptors 2 and 4 in Blood Cells of Keratoconus Patiens. Sci. Rep. 7, 12975. doi:10.1038/s41598-017-13525-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Sorkhabi, R., Ghorbanihaghjo, A., Taheri, N., and Ahoor, M. H. (2015). Tear Film Inflammatory Mediators in Patients with Keratoconus. Int. Ophthalmol. 35, 467–472. doi:10.1007/s10792-014-9971-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Gao, X., Mu, B.-k., and Wang, Y. (2021). Understanding the Role of Corneal Biomechanics-Associated Genetic Variants by Bioinformatic Analyses. Int. Ophthalmol. doi:10.1007/s10792-021-02081-9

CrossRef Full Text | Google Scholar

Sutton, G., Madigan, M., Roufas, A., and McAvoy, J. (2010). Secreted Frizzled-Related Protein 1 (SFRP1) Is Highly Upregulated in Keratoconus Epithelium: a Novel Finding Highlighting a New Potential Focus for Keratoconus Research and Treatment. Clin. Exp. Ophthalmol. 38, 43–48. doi:10.1111/j.1442-9071.2009.02216.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Lyon, D., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING V11: Protein-Protein Association Networks with Increased Coverage, Supporting Functional Discovery in Genome-wide Experimental Datasets. Nucleic Acids Res. 47, D607–D613. doi:10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Tréchot, F., Angioi, K., Latarche, C., Conroy, G., Beaujeux, P., Andrianjafy, C., et al. (2015). Keratoconus in Inflammatory Bowel Disease Patients: A Cross-Sectional Study. Eccojc 9, 1108–1112. doi:10.1093/ecco-jcc/jjv151

CrossRef Full Text | Google Scholar

Ward, P. P., Mendoza-Meneses, M., Park, P. W., and Conneely, O. M. (2008). Stimulus-dependent Impairment of the Neutrophil Oxidative Burst Response in Lactoferrin-Deficient Mice. Am. J. Pathol. 172, 1019–1029. doi:10.2353/ajpath.2008.061145

CrossRef Full Text | Google Scholar

Wu, W., Zong, J., Wei, N., Cheng, J., Zhou, X., Cheng, Y., et al. (2018). CASH: a Constructing Comprehensive Splice Site Method for Detecting Alternative Splicing Events. Brief. Bioinform. 19, 905–917. doi:10.1093/bib/bbx034

PubMed Abstract | CrossRef Full Text | Google Scholar

You, J., Hodge, C., Wen, L., McAvoy, J. W., Madigan, M. C., and Sutton, G. (2013a). Tear Levels of SFRP1 Are Significantly Reduced in Keratoconus Patients. Mol. Vis. 19, 509.

Google Scholar

You, J., Corley, S. M., Wen, L., Hodge, C., Höllhumer, R., Madigan, M. C., et al. (2018). RNA-seq Analysis and Comparison of Corneal Epithelium in Keratoconus and Myopia Patients. Sci. Rep. 8, 389. doi:10.1038/s41598-017-18480-x

PubMed Abstract | CrossRef Full Text | Google Scholar

You, J., Wen, L., Roufas, A., Madigan, M. C., and Sutton, G. (2013b). Expression of SFRP Family Proteins in Human Keratoconus Corneas. PLOS ONE 8, e66770. doi:10.1371/journal.pone.0066770

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, G., Soufan, O., Ewald, J., Hancock, R. E. W., Basu, N., and Xia, J. (2019). NetworkAnalyst 3.0: a Visual Analytics Platform for Comprehensive Gene Expression Profiling and Meta-Analysis. Nucleic Acids Res. 47, W234–W241. doi:10.1093/nar/gkz240

PubMed Abstract | CrossRef Full Text | Google Scholar