Identification and Verification of Hub Biomarkers and Immune-Related Pathways Participating in the Trabecular Meshwork after using Corticosteroid

Abstract

Background: The corticosteroids is associated with increased Intraocular Pressure (IOP), especially in the way of topical application. However, the cause and potential molecular events are not clearly explained. It was reported that immune cells may impact on matrix metalloproteinase pathway and IOP. The purpose of this study was to identify the key biomarkers and immune-related pathways involved in the Trabecular Meshwork (TM) after corticosteroids use.

Methods: The gene expression omnibus database was used to retrieve the expression profile for GSE124114 and GSE37474. Based on differential expression analysis (DEGs), hub markers for the possible molecular pathways in the TM following the use of corticosteroids were mined. The hub gene modules linked to higher IOP were found using Weighted Gene Co-Expression Network Analysis (WGCNA), and the immune cells presence of the TM was assessed using CIBERSORT.R (version 3.6.1) was used to carry out enrichment analysis on DEGs. The STRING database created the Protein-Protein Interaction (PPI) network of DEGs. The combined GSE6298 and GSE65240 datasets was used to confirm the expression of hub genes, and receiver operating characteristic curve analysis was also carried out.

Results:A total of 30 DEGs were recognized. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses demonstrated that these DEGs focused primarily on positive regulation of cytokine production and phenylalanine metabolism relevant signaling pathways. Two hub modules were enriched on rheumatoid arthritis pathway and the AGE-RAGE signaling pathway in diabetic complications. The 2 most firmly related hub genes (TSC22D3 and FKBP5) among 24 overlapped hub genes were identified from the PPI network. The immune infiltration results revealed the most significant relationship was macrophages M0. TSC22D3 was strongly related with macrophages M0 (R=0.75, p=0.018). ROC curve analysis demonstrated FKBP5 gene was important in TM treated with steroid hormone. FKBP5 gene was verified through the consolidated GSE6298, GSE65240 database.

Conclusion: Two key genes (TSC22D3, FKBP5) are involved in the comprehension of the molecular mechanisms in corticosteroids-induce ocular hypertension. TSC22D3 was strongly related to macrophages, which was associated with the pathogenesis of TM. FKBP5 could be a potential novel diagnostic marker in blood samples from patients with higher IOP.

Keywords

Glaucoma; Trabecular meshwork; Ocular hypertension; Biomarker; Immune cell infiltration; Weighted gene co-expression network analysis

Introduction

Glaucoma is the leading cause of global irreversible blindness, affecting approximately 76 million people in the year 2020 and more than 111 million worldwide in 2040 [1]. Ocular hypertension is an independent risk factor in glaucoma development and progression. It is classified into two group’s primary and secondary ocular hypertension. Primary ocular hypertension usually happens due to the special anatomy, while there are many reasons for secondary ocular hypertension, including ocular inflammation and trauma, pigment dispersion and exfoliation, neovascularization, dense cataract formation, corneal pathologies and the use of corticosteroids [2]. Ocular inflammation is the most important cause of secondary elevated Intra Ocular Pressure (IOP). Topical glucocorticoid is one of the most commonly immunosuppressive or anti-inflammatory treatments. However, it was reported that one third of the general population showed an increased IOP after using topical corticosteroids for two weeks or more [3], which was explained that corticosteroids may cause a decrease regulated outflow since a corticosteroid receptor has been detected in the trabecular meshwork [4,5]. It was widely accepted that corticosteroids would cause elevated IOP, but the explicit mechanism has not been clarified.

Corticosteroids is commonly used in clinical ophthalmology. For example, uveitis, as an immune-related disease, require long-term corticosteroid treatment. However, prolonged steroid application usually causes secondary glaucoma at the expense of elevated IOP [6]. Be that as it may, there was a limitation of monitoring mode during long-term use of corticosteroids-related eye drops until IOP spike happened. To decrease IOP, stopping using hormones seems helpful, but the disease may reoccur. So, it is important to balance corticosteroid use with the prevention of complications. Therefore, it is critical to study the potential molecular mechanisms of corticosteroids-induced ocular hypertension and consequently identify more valid diagnostic techniques and more reliable molecular markers for early prevention of increasing IOP during hormone therapy time. Since glucocorticoids mainly play a role through glucocorticoid receptors in regulating the expression of target genes, it seems important to explore the interaction with glucocorticoid application, changes in gene expression and secondary ocular hypertension. Gene expression microarrays are commonly used to research gene expression profiles, which is a relatively novel approach to examining genes and has a wide range of applications for drug-based molecular targeting. Gene Expression Omnibus (GEO) [7] has just published a considerable quantity of data. Integrating these datasets can also help researchers learn more about molecular mechanisms.

Numerous novel computational strategies have been created to aid in the identification of disease-specific biomarkers as next-generation sequencing technology has advanced. Weighted Gene Co-Expression Network Analysis (WGCNA), a novel systems biology method, facilitates the creation of free-scale gene co-expression networks and the identification of gene modules and hub genes [8]. We can determine which modules are connected to disease phenotypes by analyzing the relationship between modules and clinical factors [9].

In this study, we downloaded four original microarray datasets (including GSE124114, GSE37474, GSE6298, and GSE65240) from the GEO database which involved a total of 41 samples, with 16 healthy controls and 25 corticosteroids co-incubated trabecular meshwork samples. GSE124114 and GSE37474 were consolidated to obtain Differentially Expressed Genes (DEGs), according to log2|FC| size as the most significant DEGs, analysis utilizing packages in R (version 3.6.1) software. In addition, an analysis of the pathways that DEGs were enriched for was carried out using the Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Weighted Gene Co-Expression Network Analysis (WGCNA) was used to identify hub gene modules associated with increased IOP. The key module with the highest level of significant correlation with corticosteroids-induced was identified, and the genes with gene significance >0.2 and module membership>0.8 were selected as hub genes in the key module. Then, DEGs and hub gene modules were cross-linked to further screen the hub gene. After that, the Protein-Protein Interaction (PPI) network was used to analyze the anticipated associations for a particular group of proteins through the STRING online database. Because GSE37474 was the unique database for the angle of the tissue (trabecular meshwork and corneoscleral tissue), after data batch correction and standardization, we chose GSE37474 to do immune infiltration by performing the CIBERSORT algorithm. GSE6298 and GSE65240 were consolidated and did batch correction for verifying key genes. In addition, the accuracy of the hub markers in monitoring corticosteroids-induced was subsequently evaluated using the Receiver Operating Characteristic curve (ROC). In this way, we primarily explored the mechanism of corticosteroid-induced ocular hypertension, and provide clues to pathological mechanisms, diagnostic and therapeutic strategies for glucocorticoid induced glaucoma.

Materials and Methods

Search strategy

We used the keyword “corticosteroids” and “glaucoma” to search the GEO data-base (https://www.ncbi.nlm.nih.gov/geo/), and there was a sum of 4 results in the GEO database. By restricting the entry type (series), study type (expression profiling by array) and tissue sources (Homo sapiens), 4 pieces of items that were not related to the purpose of this study were excluded. After that, 4 series from 4 platforms of GPL6244, GPL17077, GPL570, GPL4564 were included, and gene expression profiles of GSE37474, GSE6298, GSE65240 and GSE124114 were downloaded.

Microarray data information

The platform for GSE37474 was GPL570   [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, which included 5 tissues (trabecular meshwork and corneoscleral) from corticosteroids-induced sample and other 5 from health samples. The platform for GSE124114 was GPL6244 [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version], which included 9 trabecular meshwork cell samples induced by corticosteroids and other 9 from the control sample. The platform for GSE65240 was GPL17077 Agilent-039494 Surprint G3 Human GE v2 8x60K Microarray 039381 (probe name version), which included 3 trabecular meshwork cell samples induced by corticosteroids and other 3 from the control sample. The platform for GSE6298 was GPL4564 stanford human cDNA SHEW, which included 9 trabecular meshwork cell samples. Platform and series matrix file(s) were downloaded from the GEO and saved as TXT files. R software (version 3.6.1) was used to process the downloaded files. Figure 1 showed the details of the selection process.

Figure 1: Serials selection process.

Integration of microarray data and DEGs

GSE124114 and GSE37474 raw datasets were incorporated for the analysis. The coordinated microarray datasets were batch normalized applying R software (https://bioconductor.org/biocLite.R) and then saved as a TXT file using limma packet (biocLite "limma") analysis. The DEGs in control and experimental samples were examined using the limma programme, which processed the operating instruction codes by using R software. The instruction code is run using R programme. Genes that were upregulated or downregulated were collected separately and used for further investigation. The R programme was used to transform the downloaded data (which included the platform and series of the matrix).

Data processing and identification of DEGs

Using the same limma programme, a volcano plot was constructed to depict all the up-regulated and down-regulated differentially expressed genes. The hierarchical clustering was then conducted using the heatmap programme and limiting the elements to 30 differentially expressed genes. A heatmap was used to visualize the findings. The ID associated with the probe name was translated to a gene symbol using the Perl programming language (version 5.30.0) and a genome-wide annotation library (http://www.bioconductor.org/packages/release/ data/annotation/html/org.Hs.eg.db.html) before being stored in a TXT file. DEGs were defined as P-values adjusted to 0.05 and log fold change (logFC)>1.5.

GO and KEGG pathway enrichment analyses of DEGs

The differentially expressed genes were used to perform GO and KEGG pathway enrichment analysis using the clusterProfiler package [10], with Q (Bonferroni-corrected p-value) <0.05 set as the statistically significant threshold to evaluate the biological functions and pathways affected by differential gene sets affecting corticosteroids-induced sample. GO keywords with p<0.05 were considered statistically significant, including Biological Processes (BPs), Cellular Components (CCs), Molecular Functions (MFs), and KEGG pathways. The "ggplot2" R tool was used to depict the GO and KEGG analysis findings.

Construction of co-expression modules by WGCNA of datasets

The gene co-expression networks were built using the GSE124114 and GSE37474 aggregate expression profile datasets [8]. We used GO and KEGG pathway enrichment studies of the target modules to learn more about the probable biological activities and pathways in the corticosteroids-related module. Meanwhile, hub gene sets with gene importance >0.2 and module membership >0.8 were examined. The hub gene sets were then intersected with DEGs for additional investigation.

GO enrichment and KEGG pathway analysis of genes in each module

The possible biological processes, cellular components, and molecular functions of genes in each module were determined using GO enrichment and KEGG pathway analysis. The hypergeometric test was also used to perform significant KEGG pathway and GO enrichment analyses, with an adjusted p-value (q value) <0.05 deemed significant. The bar graphs were created using the top 5 GO enrichment findings, while the GO and KEGG pathway analyses were shown using the "ggplot2" R programme.

PPI network and the identification of intersection genes between hub gene and DEGs

The search tool for the retrieval of interacting genes (STRING, https://string-db.org/cgi/input.pl) is an online database resource that includes over 24.6 million proteins and over 3.1 billion interactions from 5090 organisms [11]. We created PPI networks of inter-section genes between hub genes and DEGs using STRING (version 11.0). The minimum needed interaction score was set to the middle confidence value (0.700).

Expression verify and receiver operating characteristic curves of the key genes

According to consolidate GSE124114 and GSE37474, Receiver Operating Characteristic curves (ROCs) were plotted using the R to assess the levels of key genes distinguishing between corticosteroids-induced trabecular meshwork and control sample. Further-more, the two key genes expression levels and diagnostic value of the key genes were validated with a separate external dataset (consolidated GSE6298 and GSE65240).

Immune infiltration by CIBERSORT analysis

After data batch correction and standardization, we chose GSE37474 to perform immune infiltration using the CIBERSORT method since it was the sole database for trabecular meshwork and corneoscleral tissue. The immune cells infiltration matrix was obtained using the CIBERSORT software [12]. Twenty-two immune cells were studied in this study: Macrophages M2, plasma cells, neutrophils, mast cells activated, T cells CD8, macrophages M1, T cells gamma delta, B cells memory, monocytes, B cells naive, T cells follicular helper, NK cells activated, dendritic cells resting, T cells CD4 memory activated, T cells CD4 naive, NK cells resting, T cells regulatory (Tregs), dendritic cells activated, eosinophils, macrophages M0, T cells CD4 memory resting and mast cells resting.

Analysis of transcription factors, and microRNA interaction

To analyze target gene transcription factors, we built a network between transcription factors and genes using the miRTarBase [13]. ChEA (https://dev.networkanalyst.ca/NetworkAnalyst/Secure/AnalysisOverview.xhtml) was also used to examine the targets of empirically proven microRNA (miRNA)-mRNA interactions.

Results

Compared with the control, an aggregate of 30 DEGs was acquired, of which 20 were up-regulated and 10 were down-regulated respectively (Figure 2). Up and down differentially expressed genes were shown in Table 1 separately.

Table 1. Up and down-regulated differentially expressed genes.

The purple spot represents the down regulated genes, while the red ones mean up regulated genes. R-heatmap software was utilized to draw a heatmap of the 20 up and 10 down-regulated DEGs, as shown in Figure 3.

Figure 2: Volcano plot of the differentially expressed genes between corticosteroids-induced and normal sample. Black points represent the adjusted P-value>0.05. Blue points represent adjusted P-value<0.05 and down-regulated genes. Red points represent adjusted P-value <0.05 and the up-regulated genes.

Figure 3: Heatmap of the top 30 DEGs according to the adjusted P-value and logFC. Red indicates higher gene expression and green indicates lower gene expression

GO term and KEGG pathway enrichment analysis of DEGs

The GO term and KEGG pathway enrichment analysis of up and down-regulated genes with an adjusted P value of <0.05 were obtained respectively. The results of the GO term in the corticosteroids group were shown in Table 2 and Figure 4. The visual analysis results of the KEGG enrichment of DEGs in the corticosteroids group were shown in Table 2 and Figures 4A and 4B. The DEGs genes were mainly enriched in positive regulation of cytokine production, regulation of extracellular matrix organization, actin filament bundle assembly, actin filament bundle organization, T-helper 2 cell cytokine production, phenylalanine metabolism. In the KEGG analysis, the DEGs were only enriched in the chemokine signaling pathway.

Table 2. TOP 5 GO and KEGG differentially expressed genes.

Figure 4A: GO enrichment analysis of differentially expressed genes.

Figure 4B: KEGG enrichment analysis of differentially expressed genes.

WGCNA for Identification of Hub Genes and PPI for key genes

The sample clustering dendrogram of consolidated GSE124114 and GSE37474 (Figure 5A) showed no obvious discrepancy between the samples incorporated into the WGCNA. We selected 8 as the optimal soft threshold power based on the scale-free topology model and the mean connectivity (Figure 5B). The gene cluster dendrogram is shown (Figure 5C), where each leaf and branch on the tree represents a gene and co-expression module, respectively. The heat map (Figure 5D) illustrates the correlation between different modules and corticosteroids-induced traits. We obtained 2 modules, in which the yellow and blue consensus module was the most relevant module with corticosteroids (correlation value=0.65, -0.74; significance level p<0.05). 86 hub genes were selected in yellow co-expression modules with gene significance >0.2 and module membership >0.8. And there were 3 genes overlapped with DEGs. 130 hub genes were selected in blue co-expression modules with gene significance >0.2 and module membership >0.8. And there were 21 genes overlapped with DEGs (Figure 6). The intersection gene between the hub module and DEGs were INHBA, RGCC and CHI3L1 in the yellow module, and ABCA6, ANGPTL7, APOD, FKBP5, GALNT15, HSD11B1, ITGBL1, LIF, LUM, MAOA, MRO, MYOC, NEDD9, OLAH, PDLIM1, PMCH, PRUNE2, SERPINA3, SLPI, TSC22D3 and ZBTB16. The minimum needed interaction score was set to the middle confidence value (0.700). After we created PPI networks of intersection genes, TSC22D3, FKBP5 were defined as the most important key gene (Figure 5E).

Figure 5A: Sample clustering dendrogram of GSE124114 and GSE37474 to detect outliers.

Figure 5B: Analysis of the scale-free fit index (left) and the mean connectivity (right) for selecting various soft thresholding powers.

Figure 5C: Clustering dendrogram for genes in corticosteroids-induced traits; each color below represents one co-expression gene module.

Figure 5D: Heatmap depicting correlations between the module and corticosteroid-induced traits.

Figure 5E: Intersection genes n revealed by PPI.

Figure 6: Venngene between the yellow module and DEGs and venngene between the blue module and DEGs.

Top 5 GO term and KEGG pathway enrichment analysis of yellow and blue module

The results of the GO term in the yellow and blue modules were shown in Figures 7A-7D. The visual analysis results of the KEGG enrichment of the yellow and blue modules were shown in Figures 7A-7D. We used the statistical method of hyper-geometric test and selected p-values <0.05 with corresponding Q-values <0.05 as the significantly enriched pathways. The yellow module genes were mainly enriched in response to hypoxia, response to decreased oxygen levels, the vascular process in the circulatory system, response to oxygen levels, and regulation of tube diameter. In the KEGG analysis, genes of the yellow module were enriched in African trypanosomiasis, AGE-RAGE signaling pathway in diabetic complications, complement and coagulation cascades. The blue module genes were mainly enriched in defense response to the virus, defense response to the symbiont, response to the virus, cell-substrate adhesion, and extracellular matrix organization. In the KEGG analysis, genes of the blue module were enriched in Rheumatoid arthritis, AGE-RAGE signaling pathway in diabetic complications, influenza A, measles, and TNF signaling pathway.

Figure 7A: GO enrichment analysis of blue module.

Figure 7B: KEGG enrichment analysis of blue module.

Figure 7C: GO enrichment analysis of yellow module.

Figure 7D: KEGG enrichment analysis of yellow module.

Validation of expression levels of hub genes in consolidated GSE6298 and GSE65240

FKBP5, and TSC22D3 were identified as the key hub genes for corticosteroid-induced status and were selected for subsequent analysis. We further visualized the expression of FKBP5 genes in the consolidated GSE6298, GSE65240 dataset and found that the expression levels of these genes were significantly higher in the corticosteroids-induced group than in the normal group (p <0.05, Figure 8). But TSC22D3 was not expressed in those datasets. Among 24 intersection genes, ANGPTL7, APOD, MAOA, and SLPI were significantly higher in the corticosteroid-induced group than in the normal group (p <0.05, Figure 8). INHBA was significantly lower in the corticosteroids-induced group than in the normal group (p <0.05, Figure 8).

Figure 8: Validation of expression levels of two hub genes in related GEO datasets.

ROC analysis of the selected two hub genes

Next, we selected FKBP5 to do Receiver Operating Characteristic (ROC) analysis. Based on the expression data of consolidated GSE124114 and GSE37474, the ROC curve was plotted to assess the monitor value of the genes for corticosteroids (Figure 9). Figure 9 depicts the results of ROC validation analysis of the gene FKBP5 (AUC=91.30%).

Figure 9: The ROC curve of five hub genes in consolidated GSE124114 and GSE37474

Immune infiltration analyses and key gene colleration with immune cell

The Heatmap (Figure 10A) result showed that macrophages M0 and macrophages M1 had a negative correlation (value=-0.71), macrophages M0 and macrophages M2 had a negative correlation (value=-0.46), macrophages M0 and T cells CD4 naive had a negative correlation (value=0.89). After GSE37474 through data batch correction and standardization, the correlation Heatmap (Figure 10B) summarized the results obtained from 10 gene expression matrix and the relative percent of the 22 immune cells was shown in Figure 10C. Compared with normal tissue, the violin plot of the immune cell showed that macrophages M0 infiltrated statistically less in corticosteroids-induced sample (P=0.032) (Figure 10D). TSC22D3 was positive related to macrophages M2 (R=0.75, p=0.018) (Figure 11).

Figure 10A: Results of CIBERSORT analysis of gene expression omnibus database. Correlation matrix of infiltration degree of immune cells in corticosteroids samples. Red indicates trends consistent with the positive correlation and blue indicates trends consistent with the negative correlation between two immune cells. The bigger size of the number of statistics data represents the more positive or negative correlation.

Figure 10B: Results of CIBERSORT analysis of gene expression omnibus database. The landscape of immune cell infiltration.

Figure 10C: Results of CIBERSORT analysis of gene expression omnibus database. The distribution of 22 immune cells in 10 gene matrix. Red indicates higher immune infiltration expression and green indicates lower expression.

Figure 10D: Results of CIBERSORT analysis of gene expression omnibus database. Violin diagram of immune cell proportions in two groups. The blue fusiform fractions on the left represent the normal group and the red fusiform fractions on the right represent the corticosteroids samples. .

Figure 11: The relationship between TSC22D3 and macrophages M2.

Construction of transcription factors, and miRNA regulatory network of the two hub genes

The network between transcription factors and genes using the miRTarBase13 showed in Figure 12A. The network between microrna (miRNA)-gene interactions derived from ChEA showed in Figure 12B.

Figure 12A: Construction of interaction network maps with transcription factors miRNAs for two hub genes TSC22D3, and FKBP5. Interaction network diagram of the two hub genes and transcription factors, where red nodes represent key genes and blue nodes represent the transcription factors corresponding to the two hub genes.

Figure 12B: Interaction network diagram of key genes and miRNAs, where red nodes represent key genes and yellow nodes represent miRNAs corresponding to the two hub genes.

Discussion

IOP is the independent risk factor in the development and progression of glaucoma [14], which is a leading cause of irreversible blindness globally. The increased IOP is associated with the application of corticosteroids, in the way of topical or system application [15]. However, when treating eye illness, corticosteroids are typically and regularly used, leading to the inevitable ocular hypertension. To control corticosteroid-related ocular hypertension, it’s an emergency to figure out the mechanism of corticosteroid-induced ocular hypertension, which has not been clearly explained. Some corticosteroid receptors have been found in the trabecular meshwork, indicating that corticosteroids may lead to decreased regular aqueous humor outflow [4]. However, there is yet no precise chemical mechanism to demonstrate how corticosteroids alter aqueous fluid outflow causing increased IOP. Some changes in gene expression and structure were previously reported after treatment of corticosteroids [4]. Furthermore, it has been demonstrated that Selective Laser Trabeculoplasty (SLT) works well to reduce cortico-steroid-induced IOP spike [5]. A biologic theory suggests that SLT treatment induces trabecular and endothelial cells to release chemotactic and vasoactive agents. These cytokines stimulate the accumulation of macrophages. Macrophages may impact on matrix metalloproteinase pathway, resulting in phagocytosis and extracellular matrix remodeling, which results in lower IOP [6]. So, it was reasonable that the immune system may influence aqueous flow. Therefore, we aimed to explore the potential pathogenesis of corticosteroid-induced molecular changes in the immune system of the trabecular meshwork, and screening hub gene for early warning of IOP increasing or exploring a potential way to lower the increased IOP.

In this work, we identified 30 DEGs that might be involved in IOP increasing after using corticosteroids. The enriched GO biological process terms assigned to these genes were mainly enriched in positive regulation of cytokine production [16,17], regulation of extracellular matrix organization [18], actin filament bundle assembly, actin filament bundle organization [19], and T-helper 2 cell cytokine production [18]. KEGG pathway analysis showed significant enrichment of phenylalanine metabolism [20]. To identify the modules that are most strongly associated with genetic changes in the corticosteroid-induced change of trabecular meshwork thereby causing the decrease of aqueous humor outflow and finally leading to ocular hypertension, all above genes were divided into two separate modules (the yellow one and the blue one) according to WGCNA algorithm analysis. The yellow module gene was mainly enriched in a series of hypoxia responses. It was reported that hypoxia stimulated a robust elevation in vasorin (a kind of glycoprotein) expression, and vasorin suppressed TNF-α-induced cell death in trabecular meshwork cells [21]. It reveals that corticosteroids may trigger a response to hypoxia and have an influence on the function of trabecular meshwork cells, which can keep the flow of aqueous humor [5]. In the KEGG analysis, these genes were mainly enriched AGE-RAGE signaling pathway in diabetic complications. One study suggested that AGE-RAGE signaling pathway in diabetic complications may be the major signaling pathway under conditions of reactive oxygen species-induced damage in TM cells [22]. According to the blue module gene, these GO seemed strongly related to the immune system, such as defense response to the virus, defense response to the symbiont, response to the virus, cell-substrate adhesion [23], and inflammatory cytokines were discovered to be increased in glaucoma patients' aqueous humor by several studies [24-27]. Furthermore, the blue module gene is enriched in an extracellular matrix organization. According to Zhou’s study, glucocorticoid modulated extracellular matrices in the trabecular meshwork and reported extracellular matrix played an important role in the development of glaucoma [28]. In the KEGG analysis, the gene of blue module were also enriched in AGE-RAGE signaling pathway in diabetic complications [22].

Further gene analysis is based on the following principles: (A) Belongs to differential expression genes. (B) Gene significance >0.2 and module membership >0.8 in WGCNA. (C) Co-occurrence in the PPI network (middle confidence value=0.700). Two hub genes were identified as the key genes for the corticosteroid sample, namely, TSC22D3 and FKBP5. According to gene expression in consolidated GSE6298 and GSE65240, FKBP5 was also significantly up-regulated. Then, we performed ROC curve analysis to evaluate the monitoring value of the FKBP5. The results demonstrated that these genes might serve as monitoring markers for the corticosteroids sample, because the AUC of these five genes was >0.9 in the consolidated GSE124114 and GSE37474.

We had also predicted transcription factors, and regulatory microRNAs associated with these two burn diagnostic markers, and these predicted molecules provide new targets for the corticosteroid-induced treatment of increasing IOP.

FKBP5 is a 51-kDa immunophilin that belongs to the family of FK506-binding proteins, originally named after their ability to bind the immunosuppressant FK506 [29]. Zhang et al. reported that FKBP5 was involved in nuclear transport of the dominant-negative receptor for glucocorticoid (GRb), suggesting protective effects of FKBP5 against glucocorticoid-induced glaucoma [30]. The magnitudes of DNA demethylation in the FKBP5 gene upon corticosteroid stimulation were variable which might affect the ability of negative feedback against glucocorticoid-induced signals [31]. Ewald et al. found a correlation of the FKBP5 DNA methylation status in response to glucocorticoid exposure between peripheral blood and the brain, the FKBP5 genes DNA methylation status in peripheral blood could be a useful biomarker for glucocorticoid-induced glaucoma [32]. Our immune infiltration analysis showed macrophage M0 was significantly lower in the corticosteroid-induced sample, meanwhile, macrophage M1 and M2 increased. Biologic effects may be more important and include some immediate responses involving the release of chemotactic and vasoactive agents, such as the cytokines interleukin-1a, interleukin-1b, and tumor necrosis factor-a. These factors are involved in the release of gelatinases, macrophage recruitment, and other activities affecting aqueous outflow directly or indirectly [33,34]. In our PPI network result, TSC22D3 was associated with FKBP5, and TSC22D3 was strongly related to macrophage M2. Exogenous glucocorticoids can up-regulate the expression of TSC22D3, a transcriptional suppressive molecule with strong immunosuppressive effects [35]. Macrophages present as different subtypes (M1 macrophages and M2 macrophages) in different tissues according to changes in their environment. M1 macrophages are capable of pro-inflammatory responses and produce pro-inflammatory cytokines (including interleukin 1 (IL-1), interleukin 12 (IL-12), tumor necrosis factor-α (TNF-α), chemokines (CXCL9-11, CCL2-5), Reactive Oxygen Species (ROS) and Nitric Oxide (NO) leading to the elimination of the pathogen. In contrast, M2 macrophages are important during anti-inflammatory responses and tissue regeneration. More specifically, M2 macrophages can polarize toward several subpopulations (M2a, b, and c), among these subpopulations, M2c macrophages polarize upon treatment with IL-10 or with stress-induced glucocorticoid. M2c macrophages produce interleukin 10 (IL-10) and transforming growth factor-β (TGF-β), which are associated with immune suppression and tissue remodeling [36,37]. Glucocorticoids may cause changes in TSC22D3 gene expression, which affect macrophage function, leading to microstructure-related trabecular meshwork and causing ocular hypertension. There was no study report the relationship of TSC22D3 in glucocorticoids-induced higher IOP, it may be a therapeutic target to suppress its high expression in glucocorticoids application in ocular disease, thus preventing elevated IOP in future.

Conclusion

In conclusion, Using multiple data set analyses, we identified genetic changes in trabecular meshwork tissue after corticosteroids-use and the biological functions of this process, and screened the module most relevant to trabecular meshwork change based on WGCNA analysis in the corticosteroids-induce sample dataset. Then, we further screened two key genes in the module within that module and then performed expression and ROC analysis to validate their monitoring efficacy. Then we found that TSC22D3 was highly associated with macrophage M2, according to results of immune infiltration, macrophages played an important role in the outflow of aqueous humor outflow. We proposed a new hypothetical pathogenetic mechanism for corticosteroid-induced ocular hypertension, shedding important new light on the question of the prevention of glucocorticoid-induced glaucoma.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article.

Author Contributions

Liwen Wang developed the research question and revised the manuscript before submission. Di Song wrote the first draft of the manuscript. All authors contributed to the development of the review protocol, data analysis, and refining of the manuscript, and approved the final manuscript.

Funding

Huzhou science technology bureau (2021GYB53).

Conflict of Interest

The authors have no conflicts of interest regarding the paper.

Acknowledgments

This work benefited from previous cohort studies. The author sincerely thanks all relevant researchers for the data shared and published.

References

1. Tham YC, et al. Global prevalence of glaucoma and projections of glaucoma burden through 2040: A systematic review and meta-analysis. Ophthalmology. 2014;121:2081-2090.

   [Crossref][Google Scholar][Pubmed]

2. Agarwal HC, et al. Secondary glaucoma. Indian J Ophthalmol. 1982;30:121-124.

   [Google Scholar]

3. Moss EB, et al. A randomized controlled trial to determine the effect of inhaled corticosteroid on intraocular pressure in open-angle glaucoma and ocular hypertension: The ICOUGH Study. J Glaucoma. 2017;26:182-186.

   [Crossref][Google Scholar][Pubmed]

4. Southren AL, et al. 5 beta-dihydrocortisol: Possible mediator of the ocular hypertension in glaucoma. Invest Ophthalmol Vis Sci. 1985;26:393-395.

   [Google Scholar][Pubmed]

5. Pizzimenti JJ, et al. Selective laser trabeculoplasty for intraocular pressure elevation after intravitreal triamcinolone acetonide injection. Optom Vis Sci. 2006;7:421-425.

   [Crossref][Google Scholar][Pubmed]

6. Thorne JE, et al. Risk of cataract development among children with juvenile idiopathic arthritis-related uveitis treated with topical corticosteroids. Ophthalmology. 2020;127:S21-S26.

   [Crossref][Google Scholar][Pubmed]

7. Barrett T, et al. NCBI GEO: Archive for functional genomics data sets-update. Nucleic Acids Res. 2013;41:D991-D995.

   [Crossref][Google Scholar][Pubmed]

8. Langfelder P, et al. WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

   [Crossref][Google Scholar][Pubmed]

9. Zhu T, et al. SDF4 is a prognostic factor for 28-days mortality in patients with sepsis via negatively regulating ER stress. Front Immunol. 2021;12:659193.

   [Crossref][Google Scholar][Pubmed]

10. Yu G, et al. ClusterProfiler: A R package for comparing biological themes among gene clusters. OMICS. 2012;16:284-287.

    [Crossref][Google Scholar][Pubmed]

11. Mering CV, et al. STRING: A database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31:258-261.

    [Crossref][Google Scholar][Pubmed]

12. Newman AM, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37:773-782.

    [Crossref][Google Scholar][Pubmed]

13. Huang HY, et al. miRTarBase 2020: Updates to the experimentally validated microRNA-target interaction database. Nucleic Acids Res. 2020;48:D148-D154.

    [Crossref][Google Scholar][Pubmed]

14. Weinreb RN, et al. The pathophysiology and treatment of glaucoma: A review. JAMA. 2014;311:1901-1911.

    [Crossref][Google Scholar][Pubmed]

15. Clark AF, et al. Glucocorticoid-induced formation of cross-linked actin networks in cultured human trabecular meshwork cells. Investig Ophthalmol Vis Sci. 1994;35:281-294.

    [Google Scholar][Pubmed]

16. Ying Y, et al. Activation of ATF4 triggers trabecular meshwork cell dysfunction and apoptosis in POAG. Aging. 2021;13:8628-8642.

    [Crossref][Google Scholar][Pubmed]

17. Takase H, et al. The presence of macrophage migration inhibitory factor in human trabecular meshwork and its upregulatory effects on the T helper 1 cytokine. 2002;43:2691-2696.

    [Google Scholar][Pubmed]

18. Li H, et al. Identification of the potential biological target molecules related to primary open-angle glaucoma. BMC Ophthalmol. 2022;22:188.

    [Crossref][Google Scholar][Pubmed]

19. Bachman W, et al. Glucocorticoids preferentially influence expression of nucleoskeletal actin network and cell adhesive proteins in human trabecular meshwork cells. Front Cell Dev Biol. 2022;10:886754.

    [Crossref][Google Scholar][Pubmed]

20. Tanwar M, et al. Mutation spectrum of CYP1B1 in north indian congenital glaucoma patients. Mol Vis. 2009;15:1200-1209.

    [Google Scholar][Pubmed]

21. Choi JA, et al. Role of vasorin, an anti‐apoptotic, anti‐TGF‐β and hypoxia‐induced glycoprotein in the trabecular meshwork cells and glaucoma. J Cell Mol Med. 2022;26:2063-2075.

    [Crossref][Google Scholar][Pubmed]

22. Li Q, et al. Label-free LC-MS/MS proteomics analyses reveal proteomic changes in oxidative stress and the SOD antioxidant strategy in TM cells. Clin Proteomics. 2022;19:12.

    [Crossref][Google Scholar][Pubmed]

23. Stepp MA, et al. Immune responses to injury and their links to eye disease. Transl Res. 2021;236:52-71.

    [Crossref][Google Scholar][Pubmed]

24. Du S, et al. Multiplex cytokine levels of aqueous humor in acute primary angle-closure patients: Fellow eye comparison. BMC Ophthalmol. 2016;16:6.

    [Crossref][Google Scholar][Pubmed]

25. Huang W, et al. Inflammation-related cytokines of aqueous humor in acute primary angle-closure eyes. Invest Ophthalmol Vis Sci. 2014;55:1088-1094.

    [Crossref][Google Scholar][Pubmed]

26. Wang Y, et al. Inflammatory cytokine profiles in eyes with primary angle-closure glaucoma. Biosci Rep. 2018;38:BSR20181411.

    [Crossref][Google Scholar][Pubmed]

27. Kuchtey J, et al. Multiplex cytokine analysis reveals elevated concentration of interleukin-8 in glaucomatous aqueous humor. Invest Ophthalmol Vis Sci. 2010;51:6441-6447.

    [Crossref][Google Scholar][Pubmed]

28. Zhou L, et al. Glucocorticoid effects on extracellular matrix proteins and integrins in bovine trabecular meshwork cells in relation to glaucoma. Int J Mol Med. 1998;1:339-385.

    [Google Scholar][Pubmed]

29. Wiederrecht G, et al. Characterization of high molecular weight FK-506 binding activities reveals a novel FK-506-binding protein as well as a protein complex. J Biol Chem. 1992;267:21753-21760.

    [Crossref][Google Scholar][Pubmed]

30. Zhang X, et al. FK506-binding protein 51 regulates nuclear transport of the glucocorticoid receptor beta and glucocorticoid responsiveness. Invest Ophthalmol Vis Sci. 2008;49:1037-1047.

    [Crossref][Google Scholar][Pubmed]

31. Matsuda A, et al. DNA methylation analysis of human trabecular meshwork cells during dexamethasone stimulation. Invest Ophthalmol Vis Sci. 2015;56:3801-3809.

    [Crossref][Google Scholar][Pubmed]

32. Ewald ER, et al. Alterations in DNA methylation of FKBP5 as a determinant of blood-brain correlation of glucocorticoid exposure. Psychoneuroendocrinology. 2014;44:112-122.

    [Crossref][Google Scholar][Pubmed]

33. Latina MA, et al. Selective laser trabeculoplasty: A new treatment option for open angle glaucoma. Curr Opin Ophthalmol. 2002;13:94-96.

    [Crossref][Google Scholar][Pubmed]

34. Pizzimenti JJ, et al. Selective laser trabeculoplasty for intraocular pressure elevation after intravitreal triamcinolone acetonide injection. Optom Vis Sci. 2006;7:421-425.

    [Crossref][Google Scholar][Pubmed]

35. Yang H, et al. Stress-glucocorticoid-TSC22D3 axis compromises therapy-induced antitumor immunity. Nat Med. 2019;25:1428-1441.

    [Crossref][Google Scholar][Pubmed]

36. Yunna C, et al. Macrophage M1/M2 polarization. Eur J Pharmacol. 2020;877:173090.

    [Crossref][Google Scholar][Pubmed]

37. Maciuszek M, et al. Cortisol metabolism in carp macrophages: A role for macrophage-derived cortisol in M1/M2 polarization. Int J Mol Sci. 2020;21:8954.

    [Crossref][Google Scholar][Pubmed]