Exploring Corneal Neovascularization: An Integrated Approach Using Transcriptomics and Proteomics in an Alkali Burn Mouse Model

Abstract

Purpose: Corneal neovascularization (CNV) impairs corneal transparency and visual acuity. The study aims to deepen our understanding of the molecules involved in CNV induced by alkali burns, facilitate a better grasp of CNV mechanisms, and uncover potential therapeutic targets.

Methods: Eighty-four mice were selected for establishing CNV models via alkali burns. On days 3, 7, and 14 after the burns, corneal observations and histological investigations were conducted. An integrated analysis of RNA sequencing (RNA-seq)-based transcriptomics and label-free quantitative proteomics was performed in both normal and burned corneas. Bioinformatics approaches, encompassing Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, were applied to discern differentially expressed genes (DEGs) and crucial signaling pathways. Four potentially CNV-related genes were validated using quantitative real-time PCR (qRT-PCR) and Western blot.

Results: Significant CNV was observed on the seventh day. Forty-one genes were differentially expressed in neovascularized corneas, with 15 upregulated and 26 downregulated at both mRNA and protein levels. Bioinformatics analysis revealed that these DEGs participated in diverse biological processes, encompassing retinol and retinoic acid metabolism, neutrophil chemotaxis, and actin filament assembly, along with significant enrichment pathways like cytochrome P450, tyrosine, and phenylalanine metabolism. The upregulation of lymphocyte cytosolic protein 1 (LCP1) and cysteine and glycine-rich protein 2 (CSRP2) genes and the downregulation of transglutaminase 2 (TGM2) and transforming growth factor-beta-induced (TGFBI) genes were confirmed.

Conclusions: We analyzed gene expression differences in mouse corneas 7 days after alkali burns, finding 41 genes with altered expression. The exact role of these genes in CNV is not fully understood, but exploring angiogenesis-related molecules offers potential for CNV treatment or prevention.

Corneal neovascularization (CNV) is a pathological condition caused by various factors, such as chemical injury, oxidative stress, infections, immune disorders, and other corneal injuries.¹ These factors trigger an inflammatory response in the cornea, leading to the formation of new blood vessels that can impair corneal transparency and visual acuity. Chemical injury is one of the most common causes of CNV, particularly in cases of corneal alkali burns.² The infiltration of neutrophils, depletion of epithelial cells, and stromal lysis are the primary mechanisms of CNV in chemical injuries.³ Additionally, other conditions, such as dry eye disease, infectious keratitis, graft rejection, and ocular tumors, also involve inflammation as a key factor in the development of CNV.⁴ 

Currently, the available treatments for CNV include topical corticosteroids, antivascular endothelial growth factor (anti-VEGF) agents, photodynamic therapy, laser photocoagulation, keratoplasty, and amniotic membrane transplantation.⁵ However, these treatments have limited clinical efficacy and may cause adverse effects. Prolonged use of corticosteroids can lead to sterile ulcers and corneal perforation, whereas anti-VEGF drugs may cause perforation of the ischemic area in severe burns.^5,6 Therefore, the existing therapeutic strategies for CNV do not adequately improve ocular function or patient outcomes. 

The murine corneal alkali burn model is a well-established method used in studying CNV.⁷ Most studies examining alkali-burned corneas have focused on RNA or protein levels. In previous related studies, various techniques, such as circRNA microarray analysis,⁸ transcriptome sequencing,⁹ and quantitative proteomics,¹⁰ were used to identify dysregulated molecules in corneal tissues after alkali exposure. However, studies on the integration and analysis of the transcriptomics and proteomics data of mouse alkali-burned corneas are limited. This has led to a growing trend in the use of integrative approaches aimed at better understanding the transcription-to-translation cascade in the context of CNV. Therefore, the interaction of the transcriptome and proteome in CNV pathogenesis needs to be further explored. 

In this study, we performed RNA sequencing (RNA-seq)-based transcriptomics and label-free liquid-chromatography tandem mass spectrometry (LC‒MS/MS) quantitative proteomics analyses on mouse corneal tissues on day 7 post-alkali burn, comparing these with unburned control tissues. The expression profiles of four genes linked to CNV were chosen and subsequently validated using quantitative real-time PCR (qRT‒PCR) and Western blot analysis. Our results may contribute novel perspectives towards comprehending CNV pathogenesis and identifying emerging therapeutic targets. 

Methods

Mice

Eighty-four female BALB/c mice that were 8 weeks old were housed in Tongji University's laboratory animal facility. They had access to a 12-hour light/dark cycle, typical mouse food, and enough water. The Association for Research in Vision and Ophthalmology's (ARVO) guidelines for the Use of Animals in Ophthalmology and Vision Research were followed in this study. The Tenth People's Hospital in Shanghai's Animal Ethics Committee gave its approval to all animal studies (ID: SHDSYY-2023-3503). In accordance with ethical guidelines, mice were euthanized by cervical dislocation under 5% isoflurane inhalation anesthesia to mitigate pain and distress. 

Alkali-Induced Corneal Injury Model

Except for the 12 mice used for monitoring and photographic analysis during alkali burns, the remaining mice were randomly assigned to 2 groups: a control group without any intervention and an NaOH group subjected to alkali burn treatment. Each experiment involved 18 mice, with 9 mice per group. On the seventh day after alkali burns, corneal tissues from the right eyes of three individuals were pooled to create one sample, and each group comprised three biological replicates for subsequent experiments. Mice were anesthetized with 5% isoflurane and maintained under anesthesia with 2% isoflurane. Before causing corneal alkali damage, a topical anesthetic was applied using 0.5% proparacaine hydrochloride (Alcon, Fort Worth, TX, USA). A circular filter paper (2.0 mm × 2.0 mm) soaked with NaOH (1 mol/L) was attached to the central cornea of the right eye for 40 seconds to induce an alkali burn. Afterward, the paper was quickly removed, and the conjunctival sac was washed entirely with a 0.9% sterile saline solution for 1 minute. After alkali burns, the mice in the experimental and control groups were both treated with ofloxacin eye drops twice daily for 3 days to prevent infection. 

In Vivo Slit Lamp Evaluation of Mice Corneas

Corneal observations, including frontal and lateral perspectives along with fluorescein staining, were conducted using a slit lamp biomicroscope at three specific time points (day 3, day 7, and day 14). Pictures were taken with a digital camera mounted on the microscope in dark conditions for ImageJ analysis. During fluorescein staining, 5 µL of sodium fluorescein was instilled into the conjunctival sac of the mice, followed by manually inducing three blinks of the eyes. Images were then captured using cobalt blue light in a darkroom. The formation of neovascularization was assessed by quantifying the pixel count within vessel-occupied regions using ImageJ software. Furthermore, the epithelial defect area and epithelial thickness were manually measured using ImageJ. 

Hematoxylin and Eosin Staining

The eye globes were harvested and fixed overnight in 10% neutral-buffered formalin to preserve the tissue. After that, slices of corneal tissue with a thickness of 5 µm were prepared using the fixed tissue that had been immersed in paraffin. A microscope (Nikon Eclipse Ci-L microscope) was used to view the stained sections after they had been subjected to the normal techniques for hematoxylin and eosin (H&E) and periodic acid-Schiff staining. 

RNA Extraction, cDNA Library Construction, and RNA-Seq

The RNA extraction procedure from selected tissues was conducted using the Total RNA Kit I (Omega, China). Using a Nanodrop (NanoDrop Technologies, Wilmington, DE, USA), the amount and quality of RNA were assessed. The NEBNext UltraTM RNA Library Prep Kit for Illumina (New England Biolabs (NEB), Ipswich, MA, USA) was used to create cDNA libraries utilizing a total of 1 µg of RNA with an optical density (OD) value between 1.8 and 2.0 per sample. Oligo(dT) magnetic beads were used to enrich polyA-tail mRNA, which was subsequently randomly fragmented with the help of divalent cations in NEB Fragmentation Buffer. 

The M-MuLV reverse transcriptase system was utilized to generate the first-strand cDNA using the fragmented mRNA as the template, and random hexamers were used as primers. RNaseH was used to breakdown RNA chains, and the DNA polymerase I system was used to create second-strand cDNA utilizing dNTPs as substrates. The end-repaired, A-tailed, and ligated double-stranded cDNA was purified before being joined with sequencing adapters. The cDNA fragments of approximately 250 to 300 bp were chosen using AMPure XP beads (Beckman Coulter, Brea, CA, USA), and the final library was created by PCR amplification, PCR product purification, and AMPure XP bead selection. The Qubit 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) was used to undertake the library's initial quantification once it had been built. An Agilent 2100 bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) was used to measure the insert size of the library after it had been diluted to 1.5 ng/µL. After the insert size was as anticipated, the quality of the library was confirmed by using qRT‒PCR to precisely measure the effective concentration of the library (which was greater than 2 nM). 

The Illumina NovaSeq platform (Illumina, San Diego, CA, USA) was used for sequencing, producing 150 bp paired-end reads. The sequencing method originated in “Sequencing by Synthesis.” The sequencing flow cells were supplemented with four fluorescently labeled dNTPs, DNA polymerase, and adaptor primers for amplification. The insertion of a fluorescently tagged dNTP causes a corresponding fluorescence to be released because each complementary strand in the sequencing cluster is expanded. The fluorescence signal was collected by the sequencer, and computer software transformed the light signal into a sequencing peak to retrieve the fragment sequence information. Gene differential expression analysis was performed using DESeq2 software, and the significance threshold for detecting differences in gene expression was set at fold change (FC) > 2 with a P value < 0.05. The raw RNA-seq datasets are available at National Center for Biotechnology Information (NCBI; PRJNA1046668). 

Protein Extraction and Label-Free Mass Spectrometry Analysis

This study was carried out in accordance with the Minimum Information about a Proteomics Experiment (MIAPE) guidelines of the Human Proteome Organization's Proteomics Standards Initiative (HUPO-PSI) for reporting proteomics investigations.¹¹ Protein extraction from tissues was carried out using a cold radioimmunoprecipitation assay (RIPA) lysis buffer (Beyotime, China) supplemented with phenylmethane sulfonyl fluoride (Beyotime, China) on ice for 15 minutes. Subsequently, the lysis mixtures were centrifuged at 14,000 × g for 10 minutes at 4°C. The resulting supernatants were then collected to determine the protein concentration using the bicinchoninic acid (BCA) protein assay kit (Beyotime, China). The proteins were stored at −80°C until preparation for analysis. For assessing protein quality, 20 µg of protein from each sample was separated on SDS‒PAGE gels (at a continuous voltage of 200 V for 40 minutes) and then stained with Coomassie blue. Each sample, comprising 100 µg of proteins, was first reduced with 100 mM DTT for 5 minutes at 100°C and then alkylated with 100 mM iodoacetamide for 30 minutes in the dark. Subsequently, the protein suspensions underwent digestion with 4 µg of trypsin (Promega) in 40 µL of 50 mM NH₄HCO₃ buffer overnight at 37°C, and the resulting peptides were collected as a filtrate. Finally, the peptide segment was desalted using a C18 column. The samples underwent chromatography using the Easy nLC system (Thermo Fisher Scientific, USA), followed by analysis with an Orbitrap Exploris 480 mass spectrometer (Thermo Fisher Scientific, USA). The peptides (2 µg) were loaded onto the C18 reverse-phase analytical column (Thermo Fisher Scientific, USA) in buffer A (0.1% formic acid), followed by separation with a linear gradient of buffer B (80% acetonitrile and 0.1% formic acid) at a flow rate of 300 nL/min. Mass spectrometry (MS) was acquired using a data-dependent top 10 method that dynamically chose the most abundant precursor ions from the survey scan (350–1200 m/z) for higher-energy collisional dissociation (HCD) fragmentation. MS1 scans were obtained with an AGC target of 300% and a maximum IT of 50 ms, acquired at a resolution of 120,000 at m/z 200. Operating in a data-dependent mode with a cycle time of 1.5 seconds, MS2 scans were obtained at a resolution of 15,000 at m/z 200, with an AGC target of 75%, a maximum IT of 35 ms, and an isolation width of 1.6 m/z. The normalized collision energy was set to 33%, and the dynamic exclusion for selected ions was set to 30 seconds. The MS raw files were analyzed using MaxQuant software (version 1.6.17.0) and searched against the database Uniprot_MusMusculus_17107_20220308_swissprot. Carba-midomethylation of cysteines was set as a fixed modification, whereas protein N-terminal acetylation and methionine oxidation were considered variable modifications for database searching. The global false discovery rate (FDR) cutoff for peptide and protein identification was set to 0.01. Protein abundance was calculated based on normalized spectral protein intensity (LFQ intensity). Proteins with FC > 2 and P value < 0.05 were considered differentially expressed proteins (DEPs). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE¹² partner repository with the data set identifier PXD047367. 

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Analyses

Based on the Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.kegg.jp/) databases, GO and KEGG analyses were conducted for differentially expressed genes (DEGs), DEPs, and the combined dataset of DEGs and DEPs using the clusterProfiler R package (version 3.4.4) (https://bioconductor.org/packages/clusterProfiler/).¹³ FDR-corrected P values below 0.05 (FDR < 0.05) were used to distinguish significantly enriched terms. The results of the enrichment analysis are presented using bar graphs. 

Protein Subcellular Localization and Domain Analysis

WoLF PSORT (https://wolfpsort.hgc.jp/) was used to predict the subcellular localization of DEPs. Additionally, the Interpro database (http://www.ebi.ac.uk/interpro/) was used to analyze the structural domain annotation of DEPs. The significance level of enrichment for functional structural domains was evaluated using Fisher's exact test. 

Association Analysis Between Transcriptomics and Proteomics

Gene expression is a complex process regulated at multiple levels. Most research indicates that there is not much expression consistency between mRNAs and the related proteins.¹⁴ Therefore, a comprehensive investigation integrating both the transcriptome and proteome can enhance our understanding of gene expression control. Similar to the approach described above, expression correlation between mRNA and proteins from these two distinct “omics” was quantified using R language, followed by enrichment analysis for GO and KEGG. We selected mRNAs with FC > 2 and proteins with FC > 1.2 (P < 0.05) to avoid having too few DEGs and DEPs, which could be deleterious for further functional analysis. 

Quantitative Real-Time PCR Analysis

Total RNA was extracted from animal samples, as described above, and then cDNA templates were synthesized using HiScript II reverse transcriptase (Vazyme, China). Subsequently, qRT‒PCR was performed using 2 × AceQ qPCR SYBR Green Master Mix (Vazyme, China) on ABI QuantStudio Dx Real-Time PCR instruments (Thermo Fisher, USA) in accordance with established protocols. The internal control gene for the primer sequences given in Supplementary Table S1 was GAPDH. The 2−ΔΔCt technique was used to calculate the relative expression of the target genes. 

Western Blot Analysis

Animal tissues were lysed, and proteins were extracted as described above. After SDS–PAGE, proteins were transferred to nitrocellulose membranes (Millipore, Burlington, MA, USA), which were subsequently blocked with a 5% BSA solution. The membranes were then incubated overnight at 4°C with primary antibodies, including anti-LCP1 (ab236280, 1:1000; Abcam, Cambridge, UK), anti-TGM2 (ab109121, 1:1000; Abcam, UK), anti-TGFBI (ab169771, 1:1000; Abcam, UK), and anti-CSRP2 (ab178695, 1:1000; Abcam, UK). Anti-β-actin (ab179467, 1:5000; Abcam, UK) and anti-GAPDH (ab181603, 1:10000; Abcam, UK) were used as protein loading controls. The specific bands were detected using the LI-COR Odyssey CLx scanner (LI-COR, USA) and further analyzed with ImageJ software. 

Statistical Analysis

Three independent replicates were used in each experiment. GraphPad Prism software (version 9.5.1, USA) was used to conduct the statistical analyses. Data were presented as the mean ± standard deviation (SD) and were compared using Student's t-test between two groups or 1-way analysis of variance among several groups. Statistics were deemed significant at P < 0.05. 

Results

Longitudinal Monitoring Reveals Elevated CNV and Epithelial Damage Following Alkali Burn

Subsequent to the establishment of the corneal alkali burn model, we performed longitudinal monitoring of the CNV over a course of 14 days, with 3 time points captured on days 3, 7, and 14 (Fig. 1A). The analysis of several parameters was performed using ImageJ software, including the neovascularization area, corneal epithelial defect area, and corneal epithelial thickness (Figs. 1B, 1C, and 1D). We detected a noteworthy elevation in the CNV area across all time points compared with normal corneas, with the peak occurring on day 7. Fluorescein staining and H&E staining revealed a central epithelial defect on day 3. A central scar then appeared on day 7, with residual defects in the periphery. By day 14, complete corneal epithelialization had been observed. H&E staining further showed histological changes, including thinning of the corneal epithelium, stromal edema, and infiltration of inflammatory cells, with day 7 exhibiting the most pronounced changes. Based on these findings, we used RNA-seq and label-free quantitative proteomics analysis to compare the transcriptomics and proteomics data of normal and neovascularized corneal tissues in mice, specifically focusing on the 7-day post-alkali burn period for further analysis. 

Figure 1.

 

Longitudinal monitoring of CNV on 3, 7, and 14 days after corneal alkali burns. (A) Representative slit lamp, fluorescein sodium staining, and H&E staining images of control and alkali-burned corneas. (B-D) Quantitative analysis of CNV area, corneal epithelial defect area, and corneal epithelial thickness. The measurements from the control were set to 1 (mean ± SD). *P < 0.05, **P < 0.01, ***P < 0.001, ****P <0.0001; ns, no significant difference.

Transcriptomic Analysis Reveals Dysregulated Gene Expression and Enriched Pathways in Alkali-Burned Corneas

Supplementary Table S2 presents a summary of the quality control data for the transcriptome in the normal and alkali-burned corneal libraries, which produced a total of 128,623,864 and 133,359,910 raw reads, respectively. After filtering the raw reads and assessing the error rate and GC content distribution, 124,505,990 and 129,864,990 clean reads, respectively, were available for subsequent analysis. To assess the reliability, the fragments per transcript kilobase per million fragments mapped (FPKM) approach was used to calculate gene expression for each sample, which was corrected for sequencing depth and gene length (Fig. 2A). Additionally, significant between-group differences and reliable within-group duplications were detected by principal component analysis (PCA; Fig. 2B). These results indicated that data from the control and NaOH groups were reliable and suitable for further analyses. 

Figure 2.

 

Through transcriptomic analysis, the DEGs between groups were determined, and functional analysis in normal corneas (the control group) and alkali-burned corneas (the NaOH group) was performed. (A) A box plot based on FPKM values was employed to visualize the distribution of gene expression for different samples, with sample names on the abscissa and log2 (FPKM+1) on the ordinate. Each box in the plot denotes 5 key statistical measures, including the maximum value, Q3 quartile, median, Q1 quartile, and the minimum value. C1/2/3 = normal corneal samples and Na7/8/9 = alkali-burned corneal samples. (B) Principal component analysis of corneal samples in control versus NaOH groups. Each colored point represents a sample. (C) Volcano plot of DEGs. Upregulated genes were denoted with red dots, whereas downregulated genes were denoted with green dots. GO (D) and KEGG Pathway (E) Enrichment Analysis of DEGs. The abscissa refers to GO categories or KEGG pathway terms, and the ordinate refers to the enrichment significance. The numbers indicated on each column represent the gene count for the respective categories.

The results of the differential analysis unveiled significant DEGs between normal and alkali-burned corneas, identifying a total of 1248 DEGs, which comprised 794 upregulated genes and 454 downregulated genes. A volcano plot was subsequently generated to visualize these DEGs (Fig. 2C). Details on the DEGs are listed in Supplementary Table S3. Further analysis of the DEGs was performed using GO function and KEGG pathway enrichment analysis. The top 10 significantly enriched GO terms in each of the biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were selected and represented by a bar plot (Fig. 2D). The BP category showed that DEGs were mainly enriched in leukocyte migration, myeloid leukocyte migration, leukocyte chemotaxis, and neutrophil migration. In the CC category, DEGs were mainly enriched in the extracellular matrix, proteinaceous extracellular matrix, extracellular matrix component, and collagen trimer. The MF category showed enriched terms in receptor ligand activity, receptor regulator activity, cytokine activity, and cytokine receptor binding. KEGG pathway analysis revealed significant changes in certain pathways, such as cytokine–cytokine receptor interaction, Staphylococcus aureus infection, malaria, osteoclast differentiation, extracellular matrix (ECM)–receptor interaction, tuberculosis, rheumatoid arthritis, the NF-kappa B signaling pathway, and the IL-17 signaling pathway (Fig. 2E). 

Label-Free Proteomics Analysis Reveals Dysregulated Protein Expression and Enriched Pathways in Alkali-Burned Corneas

Label-free proteomics yielded 26,545 unique peptides and 4228 protein groups. The protein identification results are shown in Figure 3A. A total of 145 DEPs were identified, with 68 upregulated and 77 downregulated in alkali-burned samples compared with normal samples. For the visual representation of these DEPs, a volcano plot was generated (Fig. 3B). Details of the DEPs are listed in Supplementary Table S4. 

Figure 3.

 

Protein identification and functional analysis of DEPs analyzed by proteomics. (A) Protein identification of normal and alkali-burned corneal tissues. (B) Volcano plot of DEPs. Upregulated proteins were represented in blue, downregulated ones in cyan, and proteins without differential expression in gray. (C) Subcellular localization. (D) Structural domain analysis. The abscissa represents the rich factor, and the ordinate represents the structural domain name. The circle color indicates the significance level of enrichment for the corresponding structural domain category (-log10 transformed P value), with darker shades representing lower P values. The bubble size corresponds to the number of DEPs. GO (E) and KEGG pathway (F) enrichment analysis of DEPs. The abscissa represents the enrichment significance, and the ordinate represents GO terms or KEGG pathways. The numbers indicated on each column correspond to the rich factor.

We conducted subcellular localization prediction and domain analysis for all DEPs in the two groups, as shown in Figures 3C and 3D, respectively. Subcellular localization showed that the DEPs in neovascularized corneas were mainly found in the extracellular region, followed by the cytosol. Structural domain enrichment analysis revealed that the ferritin-like superfamily and calgranulins were among the most significantly enriched structural domains in the DEPs, with the serpin family being the most enriched structural domain. 

GO enrichment analysis of DEPs is presented in a histogram (Fig. 3E). In the BP classification, the DEPs were mainly enriched in negative regulation endopeptidase activity, fibrinolysis, and response to bacterium. In the CC analysis, the highest proportion of DEPs was located in the extracellular space, followed by a specific granule and chylomicron. For MF, the predominant DEPs were related to serine-type endopeptidase inhibitor activity, cysteine-type endopeptidase inhibitor activity, and structural constituent of the eye lens. Through KEGG enrichment analysis, we found that complement and coagulation cascades, drug metabolism–cytochrome P450, metabolism of xenobiotics by cytochrome P450, drug metabolism–other enzymes, and tyrosine metabolism were the main enriched pathways for DEPs (Fig. 3F). 

Transcriptomics and Proteomics Analyses Reveal Consistent Genes and Enriched Pathways in Alkali-Burned Corneas

We used a Venn diagram and four-quadrant plots to investigate the correlation between transcriptomics and proteomics data of normal and alkali-burned corneas (Figs. 4A, 4B). Our investigation revealed that the expression of mRNA and proteins is not linearly correlated. Only 41 genes showed a consistent pattern of mRNA and protein expression, either being upregulated or downregulated simultaneously. Among these genes, 15 were upregulated, whereas 26 were downregulated in the CNV model groups (see the Table). 

Figure 4.

 

Conjoint analysis of transcriptome and proteome differences between normal and alkali-burned corneal tissues. (A) Venn diagram of the overlapping molecules between the proteome and transcriptome datasets from NaOH and control group. (B) Four-quadrant plot analysis. Green dots = significant expression changes at the protein level only. Purple dots = significant expression changes at the mRNA level only. Red dots = the opposite trends at the protein and mRNA levels. Blue dots = the same trends at the protein and mRNA levels. GO (C) and KEGG pathway (D) enrichment analysis of 41 DEGs. The abscissa represents the enrichment significance, and the ordinate represents GO terms or KEGG pathways.

Figure 5.

 

(A–D) Western blot analysis of LCP1, CSRP2, TGM2, and TGFBI, with the relative protein expression level of the control set to 1 (mean ± SD). (E–H) The qRT‒PCR analysis of LCP1, CSRP2, TGM2, and TGFBI, with the relative mRNA expression level of the control set to 1 (mean ± SD). *P < 0.05, **P < 0.01, and ***P < 0.001, ****P <0.0001, ns: no significant difference.

Table.

 

List of the 41 Significantly Changed Genes With 15 Upregulated and 26 Downregulated

 

Based on these 41 DEGs that exhibited a consistent trend, we carried out GO and KEGG analyses. Significantly, the most prominently enriched BP terms encompassed retinol and retinoic acid metabolic processes, along with leukocyte migration in inflammation. Subsequent classification by CC revealed that a majority of the DEGs were extracellular, including collagen-containing extracellular matrix, intermediate filament, and lipid droplet. Furthermore, examination of MF enrichment demonstrated that these genes were predominantly involved in actin filament binding, retinol dehydrogenase activity, and carboxylic ester hydrolase activity (Fig. 4C). The KEGG pathway analysis identified the top 10 enriched pathways that were significantly associated with the DEGs, including drug metabolism–cytochrome P450, metabolism of xenobiotics by cytochrome P450, chemical carcinogenesis, tyrosine metabolism, glycolysis/gluconeogenesis, drug metabolism–other enzymes, phenylalanine metabolism, tight junction, and beta-alanine metabolism, as depicted in Figure 4D. 

Expression Analysis of Angiogenesis-Related Genes and Proteins Using qRT‒PCR and Western Blot

Based on the findings from bioinformatic analysis, we selected 4 candidate DEGs involved in CNV, including LCP1, CSRP2, TGM2, and TGFBI, for further validation using qRT‒PCR and Western blotting. The qRT‒PCR analysis revealed a significant upregulation of LCP1 and CSRP2, whereas TGM2 and TGFBI exhibited significant downregulation in alkali-burned corneal tissues compared with the control group. Moreover, Western blot analysis demonstrated that LCP1 and CSRP2 expression levels were increased following alkali burns, whereas the expression levels of TGM2 and TGFBI were markedly reduced (Fig. 5). These findings confirm the accuracy and effectiveness of the transcriptomics and proteomics analyses conducted in the present study, as the expression trends of candidate genes remained consistent at both the transcript and protein levels. 

Discussion

Multi-omics analysis presents promising avenues for comprehending the intricate information flow underlying diseases, in contrast to a unidimensional focus on single omics studies. The integration of varied omics data is frequently used to clarify potential pathological alterations contributing to diseases or to identify therapeutic targets, which can subsequently be tested in molecular investigations.¹⁵ CNV is one of the most vision-threatening complications in patients with corneal alkali burns. Currently, effective treatments for CNV are lacking, and in-depth investigation is necessary to discern novel mechanisms and therapeutic targets. Recent advancements have indicated that integrated bioinformatics analysis can expedite and enhance this process. Our integrated transcriptomic and proteomic analyses offer a comprehensive view of the complex pathophysiology of corneas in response to alkali burns, surpassing the insights provided by individual proteomic or transcriptomic approaches. By combining data from various molecular levels, including mRNA and protein expression, bioinformatics analysis can reveal critical pathways and networks that may otherwise remain elusive to traditional methods. 

Our findings imply that CNV is a time-dependent phenomenon in a mouse model of alkali-induced CNV and that early therapeutic interventions may prevent its progression. For further investigation, we conducted a comprehensive analysis of the transcriptomics and proteomics data obtained from mouse corneas collected on day 7 after alkali burn injury. Consistent with previous reports,¹⁶ our data have demonstrated a nonlinear rather than a simple one-to-one correspondence between mRNA levels and protein expression. Remarkably, we identified 41 DEGs with FC >2 at the mRNA level and > 1.2 at the protein level, which are potentially involved in CNV progress. 

Multiomics analyses have shown that inflammatory cell migration and ECM remodeling are sequential pivotal events during the development and progression of CNV. Moreover, metabolic processes and cellular signaling undergo complex changes in response to CNV. The significant enrichment of retinol and retinoic acid metabolism suggests their essential roles in corneal repair and the inhibition of CNV following alkali injury. A study by Comptour et al.¹⁷ demonstrated the specific interaction between all-trans retinoic acid (the active form of vitamin A) and LOXL4, a gene involved in ECM dynamics, thereby contributing to epithelial repair. Furthermore, the enrichment of molecular functions associated with retinol dehydrogenase suggests that corneal injury and CNV might influence the visual cycle. The cytochrome P450 pathway, which is responsible for the detoxification and elimination of various compounds, was identified as an important pathway, indicating significant metabolic changes in the cornea during CNV occurrence that could affect the pharmacokinetics and toxicity of topically administered drugs.¹⁸ Tight junctions that regulate the transport of ions, water, and solutes across the corneal epithelium while maintaining its barrier function were also found to be enriched in neovascularized corneas.¹⁹ Additionally, pathways involved in amino acid metabolism, such as tyrosine and phenylalanine metabolism, exhibited enrichment, suggesting that alterations in protein synthesis, degradation, and turnover may play critical roles in CNV following alkali burns. 

LCP1 is an actin-binding protein that is specifically expressed in leukocytes.²⁰ LCP1 participates in integrin activation and facilitates leukocyte adhesion²¹ while promoting the assembly of NLRP3 inflammasomes.²² Recent investigations have also revealed the potential of LCP1 as a therapeutic target for cerebral ischemia–reperfusion injury.²³ In our study, LCP1 exhibited low expression in normal corneas but showed significant upregulation in CNV corneas, possibly because of the low density of immune cells at immune-privileged sites.²⁴ However, following corneal alkali burn, the influx of immune cells, particularly neutrophils, migrated to the injured area, potentially influencing LCP1 expression.²⁵ 

CSRP2 is known for its role in angiogenesis and cytoskeleton formation.²⁶ Studies indicate its association with collagen secretion in pulmonary artery smooth muscle cells and its upregulation in hypoxia-induced models of pulmonary hypertension.²⁷ Our results indicated that alkali exposure led to increased CSRP2 expression, highlighting its potential contribution to CNV. 

TGM2, a multifunctional enzyme, was downregulated in alkali-injured corneas. Although the understanding of TGM2 in the cornea is limited, it has been implicated in various diseases. In gastric cancer vascular endothelial cells, TGM2 acts as the candidate receptor for the antiangiogenic GX1 peptide.²⁸ Conversely, in colorectal cancer, TGM2 promotes tumor progression, and inhibiting it attenuates cell viability and angiogenesis while inducing cell apoptosis.²⁹ In glaucoma, TGM2 contributes to increased outflow resistance in the trabecular meshwork through ECM protein cross-linking.³⁰ Moreover, TGM2 enhances the adhesion, spreading, and migration in corneal epithelial cells.³¹ Our study suggests that the downregulation of TGM2 in alkali-burned corneas might serve as a protective mechanism against excessive inflammation and scarring, contrasting its roles in promoting angiogenesis in certain cancers and glaucoma. 

TGFBI emerged as a regulator of angiogenesis in different tissue contexts. Studies have revealed its role in modulating angiogenesis during the critical window of vascular development in early alveolarization, where it promotes pulmonary endothelial angiogenesis and coordinates alveolar growth via NF-κB signaling.³² However, emerging evidence suggests that TGFBI may also have a distinct role in regulating angiogenesis in adipose tissue. Specifically, TGFBI deletion can promote adipose angiogenesis in mice, leading to an increase in vessel density and a CD31-positive region in adipose tissue.³³ These findings suggest a dual role for TGFBI in modulating angiogenesis in different tissue contexts. Our study found that TGFBI was downregulated in the alkali-burned corneas of mice, indicating the potential role for TGFBI in regulating CNV. 

In summary, the innovative integration of RNA-seq and label-free quantitative proteomics analysis allowed for a comprehensive investigation of normal and neovascularized corneal tissues in mice, resulting in the identification of several DEGs that may be closely associated with CNV. The expression patterns of 4 key genes—LCP1, CSRP2, TGM2, and TGFBI—were successfully validated using qRT‒PCR and Western blot techniques. Furthermore, through bioinformatics analysis, we elucidated the biological functions attributed to these DEGs, aiding in a deeper understanding of the pathophysiology of CNV. It is imperative to acknowledge that this study primarily characterizes the transcriptomic and proteomic landscape of alkali-burned corneas at the peak of neovascularization (7 days). However, the mechanisms underlying the differential expression of various molecules during the dynamic development of CNV remain to be explored. Additionally, the inconsistent patterns between mRNA and protein levels may result from factors like translation efficiency, protein stability, posttranscriptional modifications, or post-translational modifications affecting protein expression.³⁴ In the current study, we investigated genes that showed concordant expression patterns between the transcriptomics and proteomics data. Furthermore, genes displaying inconsistent trends at both levels are also of significant research interest. Further research is necessary to fully comprehend the intricate spatial and temporal dynamics of these molecular mechanisms, which will ultimately pave the way for exploring their therapeutic potential in clinical settings. 

Acknowledgments

The authors gratefully thank Dexi Bi for his experimental guidance and Yefei Shi and Chang Liu for their help in revising the article. 

Supported by Project of Shanghai Science and Technology (20ZR1443600). 

Disclosure: W. Wang, None; M. Deng, None; M. Li, None; L. Liu, None; J. Zou, None; Y. Qian, None