Characteristics of Autophagy-Related Genes, Diagnostic Models, and Their Correlation with Immune Infiltration in Keratoconus

Introduction

Keratoconus (KTCN) is a corneal dystrophy marked by a central or paracentral corneal cone-like protrusion that causes uneven astigmatism and different degrees of visual impairment.^1–3 Acute corneal edema is a potential complication of serious cases, which usually leaves a corneal scar, and it is one of the primary justifications for corneal transplantation.^4–6 The incidence and prevalence of KTCN vary widely globally due to differences in race, geographic location, and diagnostic criteria, ranging from 0.2 to 4790 per 100,000 people and 1.5 to 25 per 100,000 people/per year.⁷ KTCN typically develops during adolescence and accelerates in the third or fourth decade of life.³ Early detection and diagnosis of preclinical KTCN are crucial for effective management, making it necessary to explore the pathogenesis of KTCN.

Numerous previous studies have shown that KTCN is a complicated disease combining both genetic and environmental components, with a prevalence of first-degree relatives amongst KTCN patients that is about 15 to 67 times greater than that of the general population.⁸ Given the critical role of genetic factors in KTCN pathogenesis, a large number of studies have sought to investigate the disease’s etiology at the DNA, RNA, and protein levels using bioinformatics analysis. So far, dozens of candidate genes have been proposed to increase disease susceptibility. Yet, the function of these genes, including the key candidate genes Visual System Homeobox 1 (VSX1) and Superoxide Dismutase 1 (SOD1) genes, is still controversial, and the common pathogenesis of KTCN is still uncertain.^9,10

Autophagy is a highly conserved cellular response to oxidative stress that involves removing damaged cellular components and organelles, providing protection against various pathologies.^11–13 It promotes cellular growth, differentiation, survival, and homeostasis.¹² Several ocular diseases, such as age-related macular degeneration (AMD), cataracts, diabetic retinopathy (DR), glaucoma, and keratoconus, have been linked to dysfunction failure in preliminary research.^14,15 Recent investigations have highlighted autophagy’s significance in the development of KTCN. Specifically, heightened expression of microtubule-associated protein light chain 3 (LC3) and secreted frizzled-related protein-1 (SFRP-1) indicates programmed cell death and autophagy activity, which might potentially contribute to corneal stromal thinning in KTCN.^16,17 Furthermore, increased expression of LC3, an autophagy marker, correlates with disease severity, indicating its potential as a marker for monitoring disease progression and staging.¹⁷ However, further research call to assess the usefulness of LC3 expression levels in disease management. In another study, reduced levels of LC3-II and lysosomal-associated membrane protein 1 (LAMP1) were observed in both the cone and peripheral regions of KTCN patients at different clinical grades.¹⁴ These results imply that dysregulation of autophagy might be crucial in KTCN pathogenesis and progression. Hence, exploring the potential autophagy-related genes (ARGs) of KTCN would present promising novel avenues for the diagnosis and therapy of KTCN.

Although KTCN has traditionally been thought of as a non-inflammatory disease, mounting data suggests that inflammatory and immunological responses might be involved. In patients with KTCN, there is an activation of T helper cell 17 (Th17), Th1, and Th2 cells in corneal tissue and peripheral blood, activation of IL-17 signaling pathways, and high expression of multiple inflammatory mediators, including interleukin 1 (IL-1), IL-6, tumor necrosis factor-alpha (TNF-α), transforming growth factor-beta (TGF-β), and matrix metalloproteinase 9 (MMP-9) in the tear fluid.^18–21 Autophagy is known to modulate immune cell function and cytokine production during immune responses,²² although some immune cells and cytokines could also inhibit autophagy signaling.²³ To investigate the immune mechanisms involved in disease development, it is crucial to assess the connection between ARGs and immune cells in KTCN.

Here, with a comprehensive assessment of the ARGs expression profile in KTCN, we aimed to identify diagnostic ARGs and construct a well-performing diagnostic model. Subsequently, we analyzed the immune-infiltrating microenvironment of KTCN and explored the correlation between these immune-infiltrating cells and ARGs for the first time. Overall, this research sheds light on the function of autophagy-related genes in the pathophysiology of KTCN and underlines their possibilities as therapeutic and diagnosis targets.

Materials and Methods

Data Collection and Processing

The “GEOquery” program fetched the relevant datasets of KTCN from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/),²⁴ including training (GSE112155 and GSE151631) and validation (GSE77938) datasets.^25–27 The dataset information is shown in Table 1. The batch effect of the GSE112155 and GSE151631 data sets was eliminated by sva to obtain the combined dataset, which contains a total of 29 KTCN corneal samples and 17 normal corneal samples.²⁸

Table 1 GEO Dataset Information List

In addition, a search of the GeneCards database (https://www.genecards.org/) for “Autophagy”, with only “Protein Coding” retained, yielded 7412 ARGs. A search for “Autophagy” in the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) yielded 574 ARGs. A total of 537 ARGs were obtained after merging and de-duplication, and the details are shown in folder Table S1.

Identification of ARGs with Differential Expression

Before and after the batch effect was eliminated, principal component analysis (PCA) was performed on the dataset.²⁹ PCA is a data dimensionality reduction method that extracts feature vectors (components) from high-dimensional data, converts them to low-dimensional data, and displays these features in a two- or three-dimensional graph. R packages of “limma” were applied to assess the combined dataset. Gene expression values of |log fold change|(|log FC|)>0.5 and adjusted p <0.05 were regarded as DEGs of KTCN. Genes with log FC > 0.5 and adjusted p < 0.05 are differentially expressed genes with upregulated expression, and genes with log FC < 0.5 and adjusted p < 0.05 are differentially expressed genes with downregulated expression. DEGs and ARGs were taken as intersections and plotted in Wayne Diagram to obtain autophagy-related differential expression genes (ARDEGs). R software packages “ggplot2” and “pheatmap” were used to draw volcano maps and heatmaps.

Construction of PPI Networks of ARDEGs

Protein-protein interaction (PPI) was formed based on the STRING database, with a minimum correlation factor greater than 0.400, and visualized using Cytoscape software.³⁰ The cytoHubba was conducted to score each node gene by five algorithms, named, Degree Correlation (Degree), Edge Percolated Component (EPC), Maximal Clique Centrality (MCC), Maximum Neighborhood Component (MNC), and Closeness Centrality (Closeness).³¹ The top 15 genes sorted by each algorithm were used as hub genes after taking the intersection.

Construction of mRNA-miRNA, mRNA-TF Interaction Network

The interaction between miRNAs and hub genes was predicted based on the miRDB database and ENCORI database (https://starbase.sysu.edu.cn/).³² Data having a Target Score > 85 in the miRDB database and mRNA-miRNA data in the ENCORI database were intersected to create the mRNA-miRNA interaction network. The hTFtarget database (http://bioinfo.life.hust.edu.cn/hTFtarget) and the CHIPBase database (version 3.0) (https://rna.sysu.edu.cn/chipbase/) are then used to identify transcription factors (TFs) that bind to hub genes. Using the RCircos program, we also investigated the chromosomal location of hub genes.^33–35

GO and KEGG Pathway Enrichment Analysis

We conducted enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms with the R package clusterProfiler,³⁶ to predict the functionality of ARDEGs in KTCN.

Diagnostic Model Based on Hub Genes

The program “glmnet” was conducted to perform least absolute shrinkage and selection operator (LASSO) regression analysis on hub genes, with parameter seed=2020, family = “binomial”.^37,38 Subsequently, logistic regression models of acquired genes were constructed by logistic regression analysis, and then the rms package was employed to generate a nomogram.³⁹ The nomogram is based on multifactor regression analysis, in which a scale is set up to characterize the variables in a multifactor regression model, and a total score is calculated to predict the probability of an event occurring. The prediction accuracy of the nomogram was assessed using the concordance index (C index) of the calibration curve. Finally, the accuracy and discriminatory power of the logistics regression model were assessed using decision curve analysis (DCA) plots by R package ggDCA.⁴⁰

Immune Infiltration Analysis (CIBERSORTx and ssGSEA)

The combined dataset’s immune cell infiltration was examined using the CIBERSORTx program.⁴¹ With the help of this analysis tool, 22 immune cells’ constitutions and population could be estimated. Besides, the accumulation of 28 immune cell types infiltrated in the combined dataset was measured using the single sample gene set enrichment analysis (ssGSEA) algorithm.⁴² To evaluate the relationship between immune cell infiltration and hub gene expression, the Pearson correlation was carried out.

Validation of Clinical Specimens

To confirm the selected genes, we enrolled 6 KTCN patients diagnosed with KTCN (case group) and 2 age sex-matched healthy controls (HC group) from the Sichuan Provincial People’s Hospital. All participants included in the study underwent a complete history assessment, general examination, and detailed ophthalmologic examination. And their corneal surface morphology, curvature variations, and central corneal thickness were obtained using pentacam (OCULUS Optikgeräte GmbH, Wetzlar, Germany). The diagnosis of KTCN was based on the 2015 keratoconus global consensus (essential clinical manifestations: abnormal posterior corneal surface heights, abnormal distribution of corneal thicknesses, and clinical noninflammatory corneal thinning).⁴³ A hospital-based health screening center gave two healthy volunteers, which were excluded from ocular and other systemic diseases. The study, which followed the Declaration of Helsinki, received approval by the hospital’s medical ethics committee. All participants signed written consent forms after being fully informed. All trial participants—cases and controls—had their venous blood drawn.

The hub genes were expressed quantitatively using the quantitative reverse transcription-polymerase chain reaction (qRT-PCR). Peripheral blood mononuclear cells (PBMC) were extracted from KTCN and HC blood samples using Ficoll solution (Cytiva, USA). Total RNA was isolated from PBMC using TRIZol Up reagent (TransGen Biotech, Beijing, China). Then the cDNA was generated according to TransScript^® All-in-One First-Strand cDNA Synthesis SuperMix for qPCR (TransGen Biotech, Beijing, China). ABI 7500 Real-Time PCR System (Applied Bioscience, Foster City, CA, USA) was utilized to analyze the amplified cDNA after using PerfectStart^® Green qPCR SuperMix(+DyeI /+Dye II). The relative expression levels were normalized to GAPDH. Table S2 includes primers (Sangon Biotech Co., Ltd, Shanghai, China). With GAPDH as the standard, the 2^−ΔΔCT method was used for calculating the relative expression of mRNA.

Statistical Analysis

R software (Version 4.1.2) serves as the foundation for all data analyses, and for comparisons of two consecutive groups, the statistical significance of normally distributed variables was estimated by an independent Student’s t-test. The Kruskal–Wallis test was employed for group comparisons involving three or more. To evaluate and examine the statistical significance between two groups of categorical variables, the chi-square test or Fisher’s exact test was used. The correlation coefficients between different molecules were computed for the findings using Spearman correlation analysis. If not specifically indicated, the statistical significance threshold of 0.05 was met by all two-sided statistical P-values.

Results

Figure 1 depicts the study’s methodology.

Figure 1 Flowchart of the Research Methodologies.

Identification and Functional Enrichment Analysis of Differentially Expressed ARGs

As shown in Figure 2A and B, the batch effect of the samples of GSE112155 and GSE151631 datasets was eliminated using the R package “sva”. The combined dataset detected 1454 DEGs in total, including 676 upregulated and 777 downregulated genes, and DEGs were visualized by the volcano map (Figure 3A). After taking the intersection of DEGs and ARGs, 42 ARDEGs were found, as shown in the Wayne map (Figure 3B). Heatmap was drawn to demonstrate the differences in expression of these identified ARDEGs (Figure 3C).

Figure 2 Data set integration. (A) Two-dimensional PCA of the merged dataset before correction. (B) Two-dimensional PCA of the merged dataset after correction.

Figure 3 Analysis of autophagy-related differentially expressed genes. (A) Volcano map of DEGs in the KTCN and Normal groups of the combined dataset. (B) Wayne diagram of DEGs and ARGs. (C) Heatmap of ARDEGs in the KTCN and Normal groups of the combined dataset. (D) A histogram representing the ARDEGs’ GO analysis findings. (E) A histogram representing the ARDEGs’ KEGG analysis findings.

The biological processes and pathways of ARDEGs were elucidated using GO and KEGG analyses. As presented in Figure 3D, in the category of biological process (BP), the ARDEGs were notably enriched in terms of regulation of autophagy, autophagy, regulation of binding, and macroautophagy. In the cellular component (CC) category, the ARDEGs were markedly concentrated in terms of endosome membrane, autophagosome, outer membrane, and microtubule. The ARDEGs were noticeably enriched in the molecular function (MF) category for GTP binding, the structural constituent of cytoskeleton, and protein serine/threonine kinase activity GTPase activity. To determine the connection between DEGs and signaling pathways, KEGG analysis was performed in Figure 3E. The pathways of amyotrophic lateral sclerosis, neurodegeneration - multiple diseases, and Alzheimer’s disease were largely enriched in the ARDEGs. Overall, ARDEG might contribute to the pathogenesis of KTCN by modulating autophagy and kinases and is mainly related to degenerative diseases of the nervous system.

The PPI Network and mRNA-miRNA, mRNA-TF Interaction Network of ARDEGs

To determine the interactions between the 42 ARDEGs, the Cytoscape program was conducted to construct and display a PPI network (Figure 4A). A total of 33 ARDEGs interacted with other genes. Then we examined the ARDEG values under the five algorithms of MCC, EPC, Degree, MNC, and Closeness using the cytoHubba insert-in, and the intersection of top 15 genes sorted by each of the five algorithms was taken as the hub genes (Figure 4B). In total, 12 hub genes were selected, including SQSTM1, UBC, IL6, GFAP, CFTR, HMOX1, PLIN2, TUBB4A, DDIT3, BAG3, BNIP3, and TGFB1 (Table S3).

Figure 4 PPI network and mRNA-miRNA, mRNA-TF interaction network. (A) PPI network of ARDEGs. (B) Venn diagram of shared genes for the top 15 ARDEGs selected under the five algorithms of MCC, MNC, EPC, Degree, and Closeness. (C and D) mRNA-miRNA (C), mRNA-TF (D) interactions network of hub genes. (E) The chromosome of hub genes Localization map. In C-D, green circular blocks are mRNAs; purple-pink circular blocks are miRNAs; Orange circular blocks are specific TFs.

Then, by developing an mRNA-miRNA, mRNA-TF interaction network, we explored the regulatory process of 12 hub genes in KTCN. As shown in Figure 4C, the mRNA-miRNA interaction network consisted of 9 targeting mRNAs and 31 miRNAs. In the mRNA-TF interaction network, 11 genes and 186 TFs interacted, among which the hub gene UBC had the most interaction relationships with TFs, targeting 160 TFs at the same time (Figure 4D). Finally, by mapping the localization of 12 hub genes on human chromosomes, we discovered the hub genes were primarily distributed on chromosomes 5, 7, 9, 10, 12, 17, 19, and 22 (Figure 4E).

Diagnostic Model Based on Hub Genes

To evaluate hub genes’ diagnostic utility, we used LASSO regression analysis (Figure 5A and B) and identified six hub genes (BAG3, CFTR, DDIT3, GFAP, PLIN2, TUBB4A) according to the optimal lambda value. Figure 5C visualizes the six hub genes’ expression levels between KTCN and the control group, and the average expression of BAG3 is the highest.

Figure 5 Construction of the diagnostic model of hub genes and diagnostic performance. (A) Value of the 12 genes’ log (Lambda) in the LASSO model. (B) The most proper log (Lambda) value in the LASSO model. (C) Forest plot of the six selected genes in the diagnostic model. (D) Nomogram of the LASSO risk score. (E) Nomogram of the expression of 6 hub genes. (F) Calibration plot of the LASSO risk score model. (G) DCA plot of the 6 hub genes. (H) DCA plot of the LASSO risk score. In the DCA plot, the x-axis denotes the threshold probability, and the y-axis, the net benefit.

Logistic regression analysis was performed separately for the LASSO risk score and the 6 hub genes’ expression in the combined dataset. The diagnostic nomogram of the LASSO risk score and expression of the six hub genes was then built by nomogram analysis (Figure 5D and E). For a visual evaluation of the diagnostic model’s performance based on the LASSO risk score, the calibration plot (Figure 5F) was developed, and the result showed that the nomogram’s projected results were reliable (C-index is 0.961). In addition, two diagnostic models were evaluated using DCA plots, as a result, the diagnostic model of the LASSO risk score was superior to the model with six hub genes, and this model might improve the outcomes of patients (Figure 5G and H).

Diagnostic Value of Hub Genes in KTCN

To assess whether individual genes could be used to distinguish KTCN from healthy controls, 12 hub genes were first analyzed for differential expression in different subgroups. In Figure 6A, the expression of BAG3, CFTR, DDIT3, HMOX1, IL6, PLIN2, TGFB1, and UBC was considerably decreased in KTCN patients, while BNIP3 expression in the KTCN group was significantly elevated. Subsequently, ROC analysis of these nine genes showed that their area under the curve (AUC) values were greater than 0.70 (Figure 6B-J), indicating that they all have certain accuracy for the diagnosis of KTCN. The AUC values of BAG3, DDIT3, HMOX1, IL6, and UBC are more than 0.80, which is considered as having excellent diagnostic specificity and sensitivity.

Figure 6 Expression of hub genes in the combined dataset. (A) Comparison graph of the grouping of hub genes in the KTCN and Normal groups in the combined dataset. (B-J) Correlation scatter plots of TGFB1 and BAG3 (C), SQSTM1 and HMOX1 (D), UBC and HMOX1 (E), UBC and SQSTM1 (F), CFTR and BAG3 (G), HMOX1 and CFTR (H), TGFB1 and CFTR (I), and TUBB4A and CFTR (J). *P < 0.05; ** P < 0.01; *** P < 0.001.

The Immune Microenvironment of Keratoconus

To reveal the potential mechanisms of immune pathways in KTCN, we first analyzed the infiltration abundance of 22 immune cell types using CIBERSORTx. The results revealed that the KTCN patients had a much higher fraction of M2 Macrophages than controls, while the proportion of M0 Macrophages and T cells regulatory (Tregs) was significantly lower than in controls (Figure 7A). Figure 7B-D shows the association analysis between these immune cells and hub genes. M0 Macrophages and HMOX1 were in the highest positive correlation; M2 Macrophages and HMOX1 had the highest negative correlation; Tregs and SQSTM1 had the highest positive correlation.

Figure 7 Immuno-infiltration analysis of the combined data set (CIBERSORTx). (A) Grouped comparison plots of immune cells in different subgroups in the combined dataset. (B-D) Lollipop charts of correlation of immune cells M0 Macrophages (B), M2 Macrophages (C), and Tregs (D) with hub genes. *P < 0.05, **P < 0.01, ***P < 0.001.

Subsequently, the differences between the infiltration of 28 immune cell types in the KTCN and Normal groups were examined using the ssGSEA approach (Figure 8A), which demonstrated that KTCN patients had a higher proportion of Gamma delta T cells, whereas the proportion of plasmacytoid dendritic cells (pDC), NK cells, central memory CD8 T cells, and effector memory CD8 T cells were observably lower than those in the controls. Correlation analysis between these immune cells was then evaluated (Figure 8B-F). Among them, NK cells and central memory CD8 T cells (r=0.608, p<0.001), and central memory CD8 T cells and plasmacytoid dendritic cells (r=0.674, p<0.001) have the highest positive correlation. The highest negative correlations were found between effector memory CD8 T cells and gamma delta T cells (r=−0.392, p=0.007), and between NK cells and Gamma delta T cells (r=−0.438, p=0.003). Besides, to investigate the link between immunity and autophagy, we further analyzed the correlation between hub genes and immune cells. As shown in Figure 8G-K, pDC has the highest positive correlations with TGFB1 (Figure 8H) and UBC (Figure 8I) at the same time; the highest negative correlations were found between CFTR and pDC (Figure 8J), and DDIT3 and Gamma delta T cells (Figure 8K).

Figure 8 Immuno-infiltration analysis of the combined data set (ssGSEA). (A) Grouping comparison plots of immune cells under different groupings of KTCN and Normal. B. Heatmap of correlation between immune cells. (C-F) Correlation scatter diagram of central memory CD8 T cells and Natural killer cells (C), central memory CD8 T cells and pDC (D), gamma delta T cell and effector memory CD8 T cell (E), and natural killer cell and gamma delta T cell (F). (G) Heatmap of correlation between hub genes and immune cells. (H-K) Correlation scatter diagram of TGFB1 and pDC (H), UBC and pDC (I), CFTR and pDC (J), DDIT3 and gamma delta T cell (K). ns, not significant, * P < 0.05, ** P < 0.01, *** P < 0.001.

Expression Verification of Hub Genes by GSE77938 Dataset and qRT-PCR in Clinical Samples

To verify whether hub genes have the same differential expression and expression trends in other datasets, we analyzed the 12 hub genes’ expressions in the GSE77938 dataset (Figure 9A). Among them, BNIP3 and CFTR had significant (p<0.05) expression differences under KTCN and normal groupings, and the expression trends were consistent with the combined dataset. And ROC curve (Figure 9B and C) showed that BNIP3 (AUC=0.875, CI: 0.679–1.000) and CFTR (AUC=0.859, CI:0.669–1.000) had certain diagnostic specificity and sensitivity on KTCN and Normal groupings.

Figure 9 Expression validation of hub genes in the GSE77938 dataset. (A) Grouping comparison plot of hub genes under KTCN grouping and normal grouping in the GSE77938 dataset. (B) ROC curve of BNIP3. (C) ROC curve of CFTR. * P < 0.05; ** P < 0.01.

In addition, from patients (n=6) and healthy controls (n=2), we obtained blood samples for qRT-PCR. The findings demonstrated that SQSTM1, DDIT3, and BINP3 mRNA expression levels were considerably downregulating in KTCN patients. While the level of HMOX1 was significantly elevated in KTCN patients. Only the trends of DDIT3 and BINP3 were consistent with our bioinformatics analysis (Figure 10A and B). These results suggest that they might be the diagnostic biomarker for KTCN.

Figure 10 Quantitative real-time PCR results for hub genes in healthy controls (HC) and KTCN patients. (A) qRT-PCR result for DDIT3 in HC and KTCN patients. (B) qRT-PCR result for BNIP3 in HC and KTCN patients. **P < 0.01, ****P < 0.0001.

Discussion

Despite more than 160 years of research exploring the pathophysiology of KTCN since it was accurately defined by Nottingham,⁴⁴ its pathogenesis remains an enigma. Although various treatment options, such as corneal cross-linking, contact lenses, and corneal transplant surgery, each has potential benefits and drawbacks, none of them can fundamentally prevent the progress of the disease. Therefore, early identification and intervention are still a great challenge to be addressed.

In recent years, the importance of autophagy has attracted more attention, as a double-edged sword, in the development of various ocular diseases.^45,46 Growing amounts of research are being done on the biomarkers of KTCN. For instance, SQSTM1, an autophagy receptor, was found to increase in the cone area of KTCN patients compared to healthy controls and non-disease areas.¹⁴ Conversely, Yildiz et al, in their study, did not observe any significant difference in SQSTM1’s level in the corneal epithelium between the KTCN group and healthy controls.⁴⁷ In contrast to the above findings, our study revealed lower expression levels of SQSTM1 in KTCN patients than in healthy controls. Moreover, several previous studies have reported significantly elevated IL-6 levels in the tears of KTCN patients, which has a significant correlation with disease severity and corneal structural damage.^48–50 A recent study identified potential genes associated with cone corneas and allergic diseases, IL6 was identified as one of 12 hub genes, suggesting that the inflammatory environment might be critical in the pathogenesis of KTCN. In our study’s validation dataset, the AUC value of IL6 was 0.829, demonstrating its excellent diagnostic effect. Furthermore, earlier studies have reported overexpressed TGF-β1 mRNA in KTCN corneal fibroblasts, which was three times more than that in non-KTCN corneal fibroblasts,⁵¹ while our analysis found no discernible difference. In integrating text mining, TGFB1 was also discovered to be one of the hub genes of KTCN. In our study, TGFB1 had certain diagnostic accuracy as a diagnostic biomarker of KTCN. To our knowledge, only Wang et al established a diagnostic model for KTCN based on immune-related DEGs and identified eight diagnostic biomarkers, including Smad7, CDKN1A, JUN, AREG, BBC3, DUSP2, map3k8, and LIF with an AUC value of 0.860 for the combination of eight genes).⁵² Yet, nothing is known regarding the expression characteristics of ARGs in KTCN.

This study identified 42 potential ARDEGs and 12 hub genes of KTCN (SQSTM1, UBC, IL6, GFAP, CFTR, HMOX1, PLIN2, TUBB4A, DDIT3, BAG3, BNIP3, and TGFB1) by bioinformatics analysis. These ARDEGs’ biological functions were assessed by biological function enrichment analysis, and the results indicated that autophagy and its regulation are significant in the pathophysiology of KTCN. In addition, the KEGG pathway was enriched in various neurodegenerative diseases, suggesting that the pathogenesis of KTCN might have similar characteristics to degenerative diseases. Subsequently, using the LASSO risk score, we constructed a nomogram, including 6 hub ARDEGs (BAG3, CFTR, DDIT3, GFAP, PLIN2, TUBB4A), and calibration plots (C-index of 0.961) and DCA validated this nomogram’s the effective performance and confirmed its ability to effectively predict patients. To further search for diagnostic biomarkers of KTCN, we compared the expression variations of 12 hub genes within the KTCN and normal groups. As a result, we identified 9 biomarkers: BAG3, BNIP3, CFTR, DDIT3, HMOX1, IL6, PLIN2, TGFB1, and UBC, all of which had great diagnostic specificity and sensitivity.

Notably, CFTR, PLIN2, DDIT3, and BAG3 showed exceptional diagnostic values in the LASSO risk score and ROC curve. Additionally, the diagnostic value of down-regulation of BINP3 expression in KTCN has been verified through the validation dataset and qRT-PCR. As a cystic fibrosis gene, cystic fibrosis transmembrane conductance regulator (CFTR) is a unique ATP-binding cassette transporter family member. It is not only an anion channel that is vital in the pathogenesis of cystic fibrosis, but also an important regulator of the development of several tumors.⁵³ CFTR is overexpressed in ocular surface epithelial cells and is vital in chloride ion-driven tear secretion.⁵⁴ CFTR activators are potential targets for current therapy for dry eye diseases (DED).⁵⁵ In our study, CFTR exhibited remarkable diagnostic performance in both test and validation datasets, which might be a valuable biomarker for the diagnosis of KTCN.

The Perilipin2 (PLIN2) gene encodes a lipid droplet protein marker whose absence enhances autophagy and affects lipid metabolism in the liver.⁵⁶ PLIN2 has been identified as a potential diagnostic and therapeutic biomarker for gastric cancer through gene regulation involved in ferroptosis and lipid metabolism.⁵⁷ Nevertheless, there is much little known about the role of PLIN2 in ocular diseases. We proposed that, for the first time, PLIN2 might be a potential biomarker of KTCN. Although our results demonstrate the value of PLIN2 in the nomogram and ROC curve analysis, it was neither verified through qRT-PCR nor in the validation dataset. Thus, further investigation is necessary to confirm its potential role as a reliable KTCN biomarker.

DNA Damage Inducible Transcript 3 (DDIT3) is a typical endoplasmic reticulum stress marker and a vital downstream gene for autophagy initiation.⁵⁸ It has been demonstrated in earlier research that the transcription factor DDIT3 regulates retinal ganglion cell damage in a variety of neurodegenerative disorders, including glaucoma.^59,60 In addition, DDIT3 is also known to be involved in mediating corneal epithelial death, thereby contributing to the pathogenesis of DED.⁶¹ In our study, DDIT3 had the highest diagnostic efficacy in the training dataset, and qRT-PCR verified its excellent ability to differentiate KTCN patients and healthy controls. It is also meaningful in the LASSO diagnostic model. We speculate that DDIT3, a pro-apoptosis gene, might control or indicate the balance between autophagy and apoptosis in KTCN; additional research is required to determine if this is cause or effect.

Bcl2-associated athanogene-3 (BAG3) is a multi-structural domain protein belonging to the BAG family that restricts apoptosis by coupling to B-cell lymphoma-2 (bcl2).⁶² Among the family, BAG3 is the only protein identified so far responsible for determining cell fate due to its important paracrine role in the autophagic process.⁶³ Extensive research has focused on BAG3’s role in various cancers, heart diseases, and degenerative diseases, given that targeting BAG3 is an emerging strategy to regulate autophagy processes.^64–67 In our established LASSO diagnostic model, BAG3 exhibited the highest average expression and was also accurate as an independent diagnostic biomarker. There are no reports about BAG3 in KTCN yet, and further exploring the autophagy regulation of BAG3 in KTCN might provide us with novel therapeutic strategies.

BCL2 interacting protein 3 (BNIP3), an atypical member of the BCL-2 family, encodes a direct receptor for mitochondrial autophagy. It regulates apoptosis through the BH3 structural domain and also induces mitochondrial autophagy through the LC3 interacting region (LIR).⁶⁸ BINP3-mediated mitochondrial autophagy has been demonstrated to be necessary for the preservation and acquisition of stem cell pluripotency.⁶⁹ BINP3 was significantly upregulated in human keratoconus cells after different subtypes of TGF-β treatment, suggesting that mitochondrial dysfunction is a key factor in KTCN.⁷⁰ However, a recent study found no discernible difference in the expression of BNIP3 between KTCN patients and the control group.⁴⁷ In our study, apart from great diagnostic efficacy in the training dataset, BINP3 also showed the best diagnostic effect in the validation dataset. What’s more, peripheral blood samples from KTCN patients and healthy controls provided further supporting evidence. Aberrant BINP3 expression might suggest the presence of mitochondrial dysregulation in KTCN patients and is expected to be a valid biomarker for KTCN. Taken together, these genes might be important markers in KTCN, and further work is needed to elucidate their potential molecular mechanisms.

To provide potential mechanistic studies for future interventions, mRNA-miRNA and mRNA-TF interaction networks have been constructed. Among these genes, UBC exhibited the highest number of interactions with TFs, and we hypothesize that it is a major regulatory target in the transcriptional regulation of autophagy-related genes, which is crucial for the development of KTCN. More research is required to explore the RNA targets and particular mechanisms of miRNAs and TFs involved in regulating DEARGs in KTCN.

Based on the vital function of immunity in KTCN, Wang’s study compared immune infiltrating cells in KTCN with controls using ssGEEA and identified 18 significantly different immune cells.⁵² However, there have been no systematic studies elucidating the interaction between autophagy and immune cells in KTCN. Therefore, we first performed a comprehensive assessment of ARDEGs and immune infiltrating cells in KTCN using CIBERSORTx and ssGSEA. Our findings revealed that the proportion of gamma delta T cells and M2 Macrophages was significantly higher in KTCN patients than in controls. However, compared to the controls, the proportion of M0 Macrophages, Tregs, pDC, NK cells, central memory CD8 T cells, and effector memory CD8 T cells was much lower. Among these cells, gamma delta T cells, Tregs, NK cells, and effector memory CD8 T cells expression was consistent with the results of Wang’s study.⁵² The difference in some results might be attributed to the different datasets we selected for analysis. Gamma delta T cells are central to local immunity in the conjunctiva,⁷¹ which highly expresses IL-17A. A previous study reported ocular surface gamma delta T cells proportion in patients with KTCN is much higher, consistent with our results.⁵⁰ Conversely, another study reported the proportion of NK cells in the ocular surface of KTCN is much higher, which requires further validation.⁵⁰ Except for gamma delta T cells and NK cells, the remaining sic immune cells must be further investigated in vivo and in vitro investigations since they have not been previously reported in KTCN. Then the correlation between hub genes and the infiltration abundance of the above eight immune cells was analyzed. Our results revealed a correlation between them, and pDC had a good positive correlation with TGFB1 and UBC respectively (r > 0.7, P < 0.001). These findings indicate that the interaction between autophagy and immune response might be vital in the pathogenesis of KTCN, and further researches need to elucidate the mechanism of their interaction.

Unavoidably, this study has some limitations. First, our study is a secondary mining and analysis of a previously published dataset, and prospective data are necessary to validate its clinical value. Second, the nomogram we created did not incorporate clinical parameters, and its predictive capacity may improve if clinical factors were added. Third, it is unclear whether the diagnostic model based on ARDEGs is specific to KTCN only. For example, it is not clear whether the model could differentiate between patients with KTCN and corneal dilatation. Fourth, the limited number of clinical samples we used and the fact that these ARDEGs’ expression was validated using patients’ blood, while corneal tissue derived from our bioinformatics analysis, might have confounded the factors in our analysis and need to be confirmed further in multicenter, large-scale clinical trials and prospective studies for further validation.

Conclusion

In conclusion, we identified four ARGs (CFTR, PLIN2, DDIT3, and BAG3) as potential characteristic biomarkers of KTCN. In addition, the KTCN diagnostic model constructed using the LASSO risk score might provide a reliable predictive method for the diagnosis of KTCN patients. Finally, immune cell infiltration and correlation analysis of KTCN suggested that the interaction between the immune response and autophagy might be involved in the pathogenesis of KTCN. These findings not only provide potential diagnostic biomarkers for KTCN patients but also enhance our understanding of the regulatory mechanisms of ARDEGs and immune cells in KTCN, which helps to develop more effective therapeutic strategies.