Phenotypic and Transcriptomics Analyses Reveal Underlying Mechanisms in a Mouse Model of Corneal Bee Sting

Corneal bee sting (CBS) is one of the most common ocular traumas and can lead to blindness. The ophthalmic manifestations are caused by direct mechanical effects of bee stings, toxic effects, and host immune responses to bee venom (BV); however, the underlying pathogenesis remains unclear. Clinically, topical steroids and antibiotics are routinely used to treat CBS patients but the specific drug targets are unknown; therefore, it is imperative to study the pathological characteristics, injury mechanisms, and therapeutic targets involved in CBS. In the present study, a CBS injury model was successfully established by injecting BV into the corneal stroma of healthy C57BL/6 mice. F-actin staining revealed corneal endothelial cell damage, decreased density, skeletal disorder, and thickened corneal stromal. The terminal-deoxynucleotidyl transferase mediated nick end labeling (TUNEL) assay showed apoptosis of both epithelial and endothelial cells. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that cytokine–cytokine interactions were the most relevant pathway for pathogenesis. Protein–protein interaction (PPI) network analysis showed that IL-1, TNF, and IL-6 were the most relevant nodes. RNA-seq after the application of Tobradex® (0.3% tobramycin and 0.1% dexamethasone) eye ointment showed that Tobradex® not only downregulated relevant inflammatory factors but also reduced corneal pain as well as promoted nerve regeneration by repairing axons. Here, a stable and reliable model of CBS injury was successfully established for the first time, and the pathogenesis of CBS and the therapeutic targets of Tobradex® are discussed. These hub genes are expected to be biomarkers and therapeutic targets for the diagnosis and treatment of CBS.


Introduction
According to multiple surveys, hymenopterous insects, including ants, bees, and wasps, cause over 100 million injuries annually worldwide, among which bee stings animal model, which provides new therapeutic strategies for the clinical treatment of CBS and even skin bee stings.

Results
To verify the validity and stability of the corneal bee sting (CBS) model, we made a comparison with clinical cases. As shown in Figures 1 and 2, significant ocular changes occurred in the BV group as compared with the Norm and VE groups. The ocular changes included corneal edema, pupil abnormalities, corneal epithelial defects, anterior chamber exudates, and cataract formation. Symptoms were persistent and resulted in corneal blindness and severe endophthalmitis, which are consistent with clinical cases and indicate successful construction of a stable and reliable CBS model. surface irritation experiments. As the primary refractive medium of the eye, damage caused by blindness is self-evident. Consequently, CBS deserves attention in the form of positively exploring scientific and effective therapy methods. In the present study, we investigate the pathological characteristics and pathogenesis of CBS by establishing a stable and reliable animal model, which provides new therapeutic strategies for the clinical treatment of CBS and even skin bee stings.

Morphological Changes in the Cornea
The general morphological changes in whole cornea induced by CBS were evaluated. Corneal edema was severe in the BV group and mild in the VE, VT, and BV groups; however, the difference was not statistically significant ( Figures 1A,C and 2A,B). Corneal thickness increased significantly in the BV and BT groups; however, the corneas in the BT  Moreover, corneal edema and epithelial defects in the VT and BT groups were significantly improved as compared with the VE and BV groups, indicating successful establishment of a model of Tobradex ® treatment after a corneal bee sting.

Morphological Changes in the Cornea
The general morphological changes in whole cornea induced by CBS were evaluated. Corneal edema was severe in the BV group and mild in the VE, VT, and BV groups; however, the difference was not statistically significant ( Figures 1A,C and 2A,B). Corneal thickness increased significantly in the BV and BT groups; however, the corneas in the BT group were thicker than those in the BV group, which is an interesting phenomenon (Figures 2A,C and 3C-Side view). Both corneal edema and thickness tended to increase and then decrease in the BV group, but only decreased over time in the BT group. These symptoms were not fully recovered at 72 h in the BV group but were absent at 72 h in the BT group. Moreover, observation of epithelial cell nuclei and cytoskeletal fluorescence staining demonstrated that corneal structure was significantly changed at 24 h after CBS. Furthermore, the typical corneal structure was no longer visible at 72 h after CBS in the BV group, with perfuse disordered nonspecific staining, while the corneal structure in the BT group was significantly changed but remained visible ( Figure 3C,D view).
Next, changes in each layer of the cornea were evaluated. Corneal epithelium, the outermost layer of the cornea, can be directly observed and the changes quantitated by sodium fluorescein staining using slit-lamp microscopy. Epithelial defects only occurred in the BV and BT groups, and the area tended to decrease over time. These defects were not fully recovered at 72 h in the BV group but were absent at 72 h in the BT group ( Figure 1B,D). Moreover, epithelial structure was significantly changed at 24 h after CBS in the BV group, and the typical epithelial structure was no longer visible at 72 h, with diffuse nonspecific staining. In the BT group, corneal epithelial cell structure was still visible in the areas with no defects ( Figure 3A). Comparison with the VE group suggests that the corneal injury caused by CBS was due to both mechanical and toxic damage, which has been widely discussed in previous studies.
Corneal stroma, which occupies most of the thickness of the cornea, plays a major role in the thickening of corneal edema. Corneal speckle distribution and corneal crosssectional area are often used to access the level and thickness of corneal edema. Corneal stroma had marked edema and was significantly thickened in the BV group, tending to increase and then decrease over time; however, these symptoms were not fully recovered at 72 h. Unexpected results were seen in the BT group, which are discussed later in combination with the whole-mount observations (Figure 2A-C). Observation of stromal nuclei and cytoskeletal fluorescence staining shows that typical stromal structure was no longer visible at 24 h in the BV group, and became more severe at 72 h, with perfuse disordered nonspecific staining. In the BT group, typical corneal stromal structure was still visible despite the altered stroma ( Figure 3A), suggesting that corneal thickening and edema had different mechanisms in the BV and BT groups.
Corneal endothelium, the innermost layer of the cornea, has been widely studied because it plays a key role in corneal functional homeostasis. Corneal endothelium had disappeared completely at 24 h in the BV group, and had not recovered at 72 h. In the BT group, nuclear pyknosis was present at 24 h, and the cytoskeleton had disappeared; however, typical endothelial structure was observed at 72 h ( Figure 3B). Moreover, the corneal endothelium is separated from the corneal stroma by a dense and impenetrable Descemet's membrane, which is composed of collagen elastic fibers, indicating that BV can easily penetrate Descemet's membrane and damage the corneal endothelium.
then decrease in the BV group, but only decreased over time in the BT group. These symptoms were not fully recovered at 72 h in the BV group but were absent at 72 h in the BT group. Moreover, observation of epithelial cell nuclei and cytoskeletal fluorescence staining demonstrated that corneal structure was significantly changed at 24 h after CBS. Furthermore, the typical corneal structure was no longer visible at 72 h after CBS in the BV group, with perfuse disordered nonspecific staining, while the corneal structure in the BT group was significantly changed but remained visible ( Figure 3C,D view). Next, changes in each layer of the cornea were evaluated. Corneal epithelium, the outermost layer of the cornea, can be directly observed and the changes quantitated by

Morphological Changes in the Anterior Chamber
Aqueous humor in the anterior chamber plays a key role in the structural and functional stability of the anterior segment and can also reflect whether the function of the anterior segment is in the normal state. Anterior chamber exudate had accumulated by 24 h in the BV group, tending to increase and then decrease, and disappearing by 72 h after CBS, which indicates that CBS may induce acute iris-ciliary body dysfunction via the strong penetration ability of BV components. Although the accumulation of anterior chamber exudate also tended to increase and then decrease in the BT group, the degree of exudation was slightly weaker than that seen in the BV group at 6 h and 12 h; however, after 12 h, the degree of exudation was slightly stronger than that of the BV group and still . This suggests that the application of Tobradex ® slowed down the process of BV damage.

Apoptosis in the Cornea
In the BV group, signs of apoptosis were seen at 24 h. The structure of the three tissues had been destroyed by 72 h, indicating that the apoptotic process was complete. In the BT group, apoptosis in the epithelium and endothelium was seen at 24 h, suggesting that Tobradex ® slowed down the progression of damage elicited by BV ( Figure 3A,B). Apoptosis in the VE and VT groups at 24 h may have been caused by mechanical damage, further indicating the double-injury factor of CBS. As a result of the data obtained from the BV group, statistical analysis of the number of apoptotic cells was not possible.

Transcriptomics Analysis in the BV Group
To ensure stability of the bee sting injury model, we performed correlation analysis between the VE and BV groups. The highest correlation coefficient was 0.990 within the sample groups and 0.723 between the groups, and the results ( Figure 4A) indicate that our model was in a relatively stable state. To understand the intrinsic pathogenic mechanism of a bee sting, we next performed RNA-seq and analyzed the differentially expressed genes (DEGs) between the VE and BV groups. Principal component analysis (PCA) was performed to detect the degree of aggregation in the corneal samples. As shown in Figure 4B, 89.2% of the differences between the six corneal samples were from PC1; therefore, both the VE and BV samples were separated by PC1. According to the selection criteria of DEGs, the selected threshold was |log2(FC)| > 1 and p.adj < 0.05, and the number of individuals meeting this threshold was 3594, of which 1886 had upregulated expression (log2(FC) > 1) and 1708 had downregulated expression (log2(FC) < -1). DEGs are displayed as a volcano plot ( Figure 4C). Subsequently, the hierarchical cluster analysis (HCA) of significantly up and downregulated genes was performed. A heatmap was plotted ( Figure 4D) by https://www.bioinformatics.com.cn (accessed on 7 May 2022), a free online platform for data analysis and visualization showing the top 20 genes for each significant difference between the two groups.
To gain further insight into the differentially expressed genes, all 1886 upregulated DEGs were classified into GO functions: cellular components (CC) and molecular functions (MF). As shown in Figure 4E, BV injection triggered a series of biological processes such as leukocyte migration, T cell activation, ribosome biogenesis, and ribosomal RNA metabolic processes. The analysis of cellular components showed that the metabolic reactions occurred mainly in the membrane microdomain. Upregulated genes related to molecular function were mainly involved in integrin and DNA catalytic activity.

KEGG Pathway Enrichment Analysis in the BV Group
The KEGG metabolic pathways associated with environmental information processing and metabolism were significantly altered after CBS ( Figure 5A). Further KEGG enrichment plots showed that the top 5 enriched pathways related to upregulated DEGs were cytokine-cytokine receptor interaction, viral protein-cytokine receptor interaction, osteoclast differentiation, small cell lung cancer, and DNA replication ( Figure 5C). Interestingly, chemokine-related (CXCL3, CXCL5, CXCL11) and interleukin-related (IL-24, IL-1b, IL-12a) genes were consistent with the GO enrichment stress and chemokine responses ( Figure 5B). Moreover, KEGG PATHVIEW showed that these chemokines were upregulated in the cytokine-cytokine receptor interaction pathway ( Figure S1). viduals meeting this threshold was 3594, of which 1886 had upregulated expression (log2(FC) > 1) and 1708 had downregulated expression (log2(FC) < -1). DEGs are displayed as a volcano plot ( Figure 4C). Subsequently, the hierarchical cluster analysis (HCA) of significantly up and downregulated genes was performed. A heatmap was plotted (Figure 4D) by https://www.bioinformatics.com.cn (accessed on 7 May 2022), a free online platform for data analysis and visualization showing the top 20 genes for each significant difference between the two groups.

RNA-Seq Analysis following Tobradex ® Treatment of CBS
In clinical practice, the preferred treatment of bee stings is topical Tobradex ® application. To further study the underlying mechanism of action of Tobradex ® , RNA-seq analysis of the BV and BT groups was performed. As shown in Figure 6A, 73.8% of the differences between the six corneal samples were from PC1; therefore, both the BV and BT samples were separated by PC1. A total of 2117 DEGs were identified by differential gene expression analysis, including 1639 upregulated genes and 478 downregulated genes ( Figure 6A,B). Subsequently, the HCA of significantly up-and downregulated genes was performed. A heatmap was plotted ( Figure 6C) by https://www.bioinformatics.com.cn (accessed on 7 May 2022), a free online platform for data analysis and visualization showing significant differences in gene expression patterns between the two groups. enrichment plots showed that the top 5 enriched pathways related to upregulated DEGs were cytokine-cytokine receptor interaction, viral protein-cytokine receptor interaction, osteoclast differentiation, small cell lung cancer, and DNA replication ( Figure 5C). Interestingly, chemokine-related (CXCL3, CXCL5, CXCL11) and interleukin-related (IL-24, IL-1b, IL-12a) genes were consistent with the GO enrichment stress and chemokine responses ( Figure 5B). Moreover, KEGG PATHVIEW showed that these chemokines were upregulated in the cytokine-cytokine receptor interaction pathway ( Figure S1).

RNA-Seq Analysis following Tobradex ® Treatment of CBS
In clinical practice, the preferred treatment of bee stings is topical Tobradex ® application. To further study the underlying mechanism of action of Tobradex ® , RNA-seq analysis of the BV and BT groups was performed. As shown in Figure 6A, 73.8% of the differences between the six corneal samples were from PC1; therefore, both the BV and BT samples were separated by PC1. A total of 2117 DEGs were identified by differential gene GO enrichment analysis shows that upregulated DEGs were highly concentrated in visual perception, perception of light stimulation, regulation of transmembrane ion transport, organization of extracellular matrix, and extracellular structure. In terms of cell components, the extracellular matrix and collagen were enriched. Molecular functions mainly involved the activities of gated ion channels ( Figure 6D,E). After Tobradex ® treatment of CBS, the main signaling pathways changed significantly, for example GABAergic synapses, morphine addiction, light signal transduction, nicotine addiction, and dopaminergic synapses ( Figure 6F,G).

Genome-Wide Gene Set Enrichment Analysis (GSEA) and PPI Network
GO and KEGG analysis mainly investigated which pathways were enriched for DEGs and tend to miss some genes that are not significantly differentially expressed but are biologically important. Therefore, genome-wide GSEA enrichment analysis was performed to further investigate the functional alterations associated with BV. Figure 7A-C shows that the genome was most significantly enriched for reactome signaling by interleukins, reactome neutrophil degranulation, and reactome GPCR-ligand binding.
Hub genes are those that play a crucial role in biological processes and are influenced by the regulation of other genes in related pathways. These genes are often important targets in disease pathogenesis and treatment; therefore, the String database and Cystoscope plug-in Cyto-Hubba were used to construct protein-protein interaction network diagrams (PPI) composed of DEGs between the two groups. The degree of protein-protein interaction determines the size and color depth of each node. The top 30 and top 10 hub genes were screened with degree as 300 standards. As shown in Figure 7D-I, TNF, IL-6, and IL-1b were hub genes in the development of corneal injury caused by bee stings, and GFAP, MAG, and MBP played key roles in the mechanism of action of Tobradex ® . expression analysis, including 1639 upregulated genes and 478 downregulated genes (Figure 6A,B). Subsequently, the HCA of significantly up-and downregulated genes was performed. A heatmap was plotted ( Figure 6C) by https://www.bioinformatics.com.cn (accessed on 7 May 2022), a free online platform for data analysis and visualization showing significant differences in gene expression patterns between the two groups.  targets in disease pathogenesis and treatment; therefore, the String database and Cystoscope plug-in Cyto-Hubba were used to construct protein-protein interaction network diagrams (PPI) composed of DEGs between the two groups. The degree of protein-protein interaction determines the size and color depth of each node. The top 30 and top 10 hub genes were screened with degree as 300 standards. As shown in Figure 7D-I, TNF, IL-6, and IL-1b were hub genes in the development of corneal injury caused by bee stings, and GFAP, MAG, and MBP played key roles in the mechanism of action of Tobradex ® .

Discussion
Although bee stings are a common trauma seen in the emergency department, CBS is rare; thus, there are no in-depth studies or standardized therapy strategies for CBS to date. CBS leads to various ocular complications, such as corneal opacity, edema, corneal hyperalgesia, iris atrophy, secondary glaucoma, cataracts, endophthalmitis, and optic neuropathy, which can critically alter vision. Accordingly, early diagnosis and treatment of these complications are essential to preventing permanent blindness [37].
In the case of CBS, most patients present with corneal edema, corneal epithelial defects, photophobia, pain, and other inflammatory responses during the acute phase [38]. In the present study, 12 h after stromal injection of PBS, corneal edema subsided and the epithelial defects recovered rapidly; however, mice injected with BV still suffered to some extent from corneal epithelial defects and turbidity after 72 h. Moreover, even following treatment with Tobradex ® , the cornea was thicker than that seen in the VE group, indicating incomplete healing [29]. Unexpectedly, residual corneal scarring and clouding of the anterior lens capsule is found in clinical cases several years after conventional treatment [39], and the density of endothelium in the affected eye is significantly reduced [29]. Therefore, it is imperative to investigate the pathogenesis of CBS with a view to exploring novel drug targets and biomarkers.
Corneal epithelium is the initial protective barrier of the eye against physical and chemical injury and pathogenic infection, thus maintaining transparency of the cornea and integrity of the visual system. Therefore, we hypothesized that even if the bee venom was extremely toxic, it could cause serious sustained damage to the cornea only after breaking through the epithelium and Bowman's membrane. Then, the bee venom in the corneal stroma breaks through the Bowman's membrane and enters the eye, causing damage to the whole anterior segment and even the whole eye. Therefore, we hypothesized that the penetration of bee venom through the Bowman's membrane plays a key role in the process of CBS. Therefore, we used a needle to simulate the penetration of bee venom through the Bowman's membrane and inject the bee venom into the corneal stroma to simulate CBS more realistically.
Corneal epithelial injury triggers cytokine-mediated interactions among epithelial, stromal, and immune cells. Damage to corneal epithelial cells leads to the release of a large number of cytokines (such as IL-1, IL-6, and TNF), which attract immune cells to migrate from limbal blood vessels to the corneal stroma. Destruction of Bowman's membrane accelerates the diffusion of cytokines into the stroma and activates target cells. Mechanical damage to intraocular tissues caused by bee stings provokes a non-inflammatory cascade reaction, while melittin, PLA2, hyaluronidase, and other active substances in BV alter the permeability of the plasma membrane and accelerate the diffusion of BV, and histamine and MCD peptides trigger allergic reactions. The immune response in the cornea increases the release of chemokines, leading to the accumulation of inflammatory cells and ultimately resulting in corneal cell death.
CBS injuries are commonly associated with pain. The cornea is the most densely innervated tissue in the body and is supplied by the ciliary nerve from the ophthalmic branch of the trigeminal ganglion. Different types of sensory neurons are present in the cornea, which are functionally divided into those expressing multimodal nociceptive receptors, cold and heat receptors, and selective mechanonociceptive receptors [40,41]. Most of the sensory nerve fibers that innervate the cornea are of the multimodal type. Therefore, we hypothesized that the corneal pain caused by CBS might be attributed to the stimulation and destruction of corneal nerves by the bee venom.
Melittin activates transient receptor potential (TRP) channels through the PLA2 cascade, sensitizing primary pain receptors and causing pain [16]. According to GO and KEGG analysis of the differentially expressed and hub genes, we suggest that melittin and PLA2 work together on PLA2R1, thereby upregulating the expression of PLA2G4C, NOS, and ALOX5AP to generate metabolites such as NO and arachidonic acid. These metabolites, in turn, selectively activate the family of TRP channels (TRPV2, TRPV4, TRPM2, and TRPM6) through the phospholipase A2-lipoxygenase (PLA2-LOX) pathway, leading to Ca 2+ influx and membrane depolarization, sensitization of primary nociceptive neurons, and ultimately pain [42,43].
Simultaneously, melittin causes tissue damage and activates pain receptors. Together with MCD, melittin induces mast cells to release mediators, such as histamine, bradykinin, and ATP, which activate G protein-coupled receptors (GPCRs) located on the terminal cell membrane of nociceptive neurons, mediating ligands of one or more signaling pathways (e.g., protein kinase A (PKA), protein kinase C (PKC), phospholipase C (PLC), or extracellular signal-regulated kinase (ERK) [44,45]). Other chemical mediators, including growth factors acting on tyrosine kinase receptors, can activate phospholipase C (PLC), phosphoinositide 3-kinase (PI3K), and mitogen-activated protein kinase (MAPK). Kinases activated by these signaling pathways, in turn, phosphorylate existing protein targets expressed in nociceptive neurons.
Apamin, another critical peptide component of BV, selectively inhibits Ca 2+ -dependent K + channels in the central nervous system. Hence, we speculate that apamin selectively in-hibits SK channels (KCNN4, KCNK1, and KCNQ5) on corneal neurons, decreasing delayed hyperpolarization of cells and enhancing sustained neuronal firing, thereby increasing cellular sensitivity to pain. Likewise, it has been reported that MCD blocks Ca 2+ -activated K + channels, thereby increasing neural excitability.
Chen et al. [46] reported that an imbalance between excitatory amino acids (EAAs) and inhibitory amino acids (IAAs) in the spinal cord is associated with the maintenance of persistent pain-related behaviors in inflammatory pain states. PLA2 activates and enhances glutamatergic excitatory synaptic transmission in glial neurons. Spinal substantia gelatinosa neurons not only accept glutamatergic excitatory transmission but also GABAergic and glycinergic inhibitory transmission.
It has been reported [47] that corneal injury can cause corneal cells and locally infiltrated immune cells to release a large number of pro-inflammatory cytokines (such as TNF, IL-6, and IL-1β), resulting in ocular inflammation [48]. In addition, peritoneal injection of BV has been demonstrated to induce neuronal cell death in mice. Furthermore, PLA2 not only directly induces neuronal death in vitro but also indirectly by activating cytokines, such as TNF and IL-1β [49]. Therefore, we speculate that in the CBS model, the injured cornea releases a large number of pro-inflammatory factors, especially TNF, under the synergistic action of melittin and PLA2, causing corneal nerve injury and activating primary sensory neurons, which leads to localized demyelination and axonal degeneration. In the BV group, the expression of injury-related genes FOS, SOCS3, and ATF3 was upregulated, indicating that BV injection not only causes a series of immunoinflammatory reactions in the cornea but also injury of primary neurons under the synergism of cytokines.
In the clinic, conventional treatment is generally preferred for CBS injuries, with immediate removal of the venomous stinger and topical steroid medication to control inflammation, with a view to protecting vision; however, some patients experience serious complications and require complex surgery. Here, RNA-seq transcriptomics analysis was performed for the first time on the cornea following the bee sting and Tobradex ® application. The results show that the expression levels of pro-inflammatory factors, such as IL-1, IL-6, IL-24, CXCL2, and CXCL3, were significantly downregulated after treatment, indicating that Tobradex ® can alleviate a series of immunoinflammatory responses caused by corneal trauma via the cytokine-cytokine receptor interaction by downregulating these pro-inflammatory factors [50].
A PPI network was constructed to identify the hub genes in the corneal healing process following a bee sting, which revealed that GFAP, MAG, PLP1, and MBP became hub genes after Tobradex ® therapy. GFAP is associated with astrocyte differentiation, MAG and MBD are involved in the formation of myelinated Schwann cells, and PLP1 is a major component of myelin sheaths. It was discovered that the astrocyte-specific marker Aldh1l1 and the unmyelinated Schwann cell markers SCM7A, MATN4, NCAM1, and GFRA3 were considerably upregulated in the BT group. KEGG enrichment analysis of the hub genes also indicates that Tobradex ® can protect corneal nerves by repairing axons. It has been reported that treatment of sciatic nerve injury with dexamethasone results in the thickening of nerve fibers and myelin sheathing, thereby accelerating axonal regeneration [51]. Accordingly, we hypothesize that Tobradex ® may play a role in corneal nerve repair by "activating" Schwann cells to promote the formation of new myelin sheathing.
In addition, this study had some limitations that could not be ignored. CBS is one of the most common ocular traumas and may lead to blindness because of the complications arising due to the toxic effects of bee venom. However, we investigated only the changes in the anterior segment, especially in the cornea, due to our experimental conditions. Additionally, we did not further verify and explore the results of bioinformatics of highthroughput sequencing.

Conclusions
In the present study, a stable CBS animal model was successfully established for the first time. Based on this model with or without Tobradex ® therapy, we demonstrated that the phenotypic changes and the active components of BV sensitized the nociceptive receptors at the corneal nerve endings and damaged the corneal nerves to some extent. The treatment with Tobradex ® exerted a protective effect on the corneal nerves by restoring damaged axons through the upregulated expression of GFAP, MAG, PLP1, and MBP, besides demonstrating a potent anti-inflammatory action.

Preparation of the BV Solution
Western honeybee (Apis mellifera)'s BV was purchased as powder from Yimin apiary (China). The powder was dissolved in 1× PBS to prepare a 0.5 mg/mL BV solution, which is the minimum concentration that can cause a severe corneal bee sting.

Establishment of the Corneal Bee Sting Model
A total of 45 mice were randomly assigned to 5 groups with 9 in each: the normal group (Norm), vehicle group (VE), vehicle with therapy group (VT), BV injection group (BV), and BV injection with therapy group (BT). Mice were anesthetized by intraperitoneal injection of 1% pentobarbital sodium solution (5 µL/g body weight), following which 0.5% proparacaine hydrochloride eye drops were applied to numb the cornea, and an eye speculum was used to expose the ocular surface. As shown in Figure 8, under the stereomicroscope, a 32G needle was used to make a tunnel through the central cornea of the right eye into the corneal stroma. In the BV group, 2 µL 0.5 mg/mL BV solution was injected into the corneal stroma with a 34G needle, and the conjunctival sac was washed three times with 1× PBS immediately afterwards. In the BT group, Tobramycin Dexamethasone eye ointment (Tobradex ® , Alcon Laboratories Inc., Fort Worth, TX, USA) was applied based on the same process as the BV group. The VE group was injected with 2 µL 1× PBS solution, and dexamethasone eye ointment was applied to the VT group based on the same process as the VE group. The eyelids were closed to protect the ocular surface. Finally, the anesthetized mice were placed on a thermostatic electric blanket to recover from the anesthesia. ulum was used to expose the ocular surface. As shown in Figure 8, under the stereomicroscope, a 32G needle was used to make a tunnel through the central cornea of the right eye into the corneal stroma. In the BV group, 2 μL 0.5 mg/mL BV solution was injected into the corneal stroma with a 34G needle, and the conjunctival sac was washed three times with 1× PBS immediately afterwards. In the BT group, Tobramycin Dexamethasone eye ointment (Tobradex ® , Alcon Laboratories Inc., Fort Worth, TX, USA) was applied based on the same process as the BV group. The VE group was injected with 2 μL 1× PBS solution, and dexamethasone eye ointment was applied to the VT group based on the same process as the VE group. The eyelids were closed to protect the ocular surface. Finally, the anesthetized mice were placed on a thermostatic electric blanket to recover from the anesthesia.

Slit-Lamp Microscopic Examination
At 6 h, 12 h, 24 h, and 72 h after the injection, the ocular changes (including corneal opacity, iris pigment changes, and pupil size) and sodium fluorescein staining (corneal epithelial defects) were observed using a slit-lamp microscope (BQ-900; Haag-Streit AG, Koeniz, Switzerland). The fluorescent area in the image was analyzed using Fiji [52] to assess epithelial injury.

Optical Coherence Tomography
Optical coherence tomography (OCT) was applied to the biological measurement of the cornea, anterior chamber angle, lens, and other anterior segment structures (RT-100;

Slit-Lamp Microscopic Examination
At 6 h, 12 h, 24 h, and 72 h after the injection, the ocular changes (including corneal opacity, iris pigment changes, and pupil size) and sodium fluorescein staining (corneal epithelial defects) were observed using a slit-lamp microscope (BQ-900; Haag-Streit AG, Koeniz, Switzerland). The fluorescent area in the image was analyzed using Fiji [52] to assess epithelial injury.

Optical Coherence Tomography
Optical coherence tomography (OCT) was applied to the biological measurement of the cornea, anterior chamber angle, lens, and other anterior segment structures (RT-100; Optovue, Fremont, CA, USA) at 6 h, 12 h, 24 h, and 72 h after injection. OCT images of the cornea and anterior chamber were analyzed using Fiji to assess the changes.

Preparation of Corneal Whole-Mounts
Mice were sacrificed by cervical dislocation at 24 h and 72 h after the establishment of the model. The eyeballs were quickly fixed in 4% paraformaldehyde solution at 4 • C for 1 h. Subsequently, the whole cornea was cut into four petals with a scalpel under a surgical microscope. The cornea was placed in an EP tube with 1× PBS solution and rinsed three times in a shaker (at slow speed) at room temperature for 10 min each. The cornea was fixed in cold acetone for 2 min, rinsed three times with 1% TD-Buffer for 10 min each, and blocked in 2% BSA solution for 1 h.

F-Actin Staining of Corneal Whole-Mounts
Alexa Fluor™ 555 Phalloidin (A34055, ThermoFisher Scientific, Waltham, MA, USA) was used to observe the overall shape and structure of the cells in the corneal whole-mount according to the method provided by ThermoFisher Scientific. Nuclei were labeled with Hoechst33342 (62249, ThermoFisher Scientific, USA). Whole-mounts were observed using a ZEISS LSM 880 confocal microscope. The images were used to create a 3D reconstruction of the cornea with IMARIS 9.2.0.

TUNEL Assay of Corneal Whole-Mounts
A colorimetric TUNEL apoptosis assay kit (C1088, Beyotime, Shanghai, China) was used to observe the level of apoptosis in corneal whole-mounts according to the method provided by Beyotime. Nuclei were labeled with Hoechst33342. Whole-mounts were observed using a ZEISS LSM 880 confocal microscope.

Whole-Genome mRNA Sequencing and Bioinformatics Analysis
Total RNA was isolated as described previously [53]. RNA was segmented by interrupting buffer, the random N6 primers were reverse-transcribed, and cDNA was synthesized. The synthetic double-stranded DNA was flattened and phosphorylated at the 5 end to form a sticky end protruding an "A" at the 3 end, and then connected to a bulbous connector with a protruding "T" at the 3 end. The ligand products were amplified by PCR using specific primers. PCR products were thermally denatured into single strands, following which the single-stranded DNA was cycled with a bridge primer to obtain a single-stranded circular DNA library. The cDNA library of all specimens in this experiment was collected and sequenced with an Illumina Hi Seq X-Ten (Huada Gene, Shenzhen, China). Hub genes were screened out comprehensively by the Cyto-scape software 3.7.2 in combination with the cyto-Hubba plug-in through 12 topological analysis methods.

Statistical Analysis
All the experiments were conducted using the randomized double-blind control method, and each experiment was repeated at least three times. Experimental data are represented as the mean ± standard deviation. GraphPad Prism 9 (San Diego, CA, USA) was used for statistical analysis and data expression. PCA was performed to detect the degree of aggregation in the samples. HCA was performed to cluster genes with similar expression patterns. Statistical analysis and visualization were performed in R version 3.6.3 and ggploT2 packages (for visualization).

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/toxins14070468/s1, Figure S1: Differences in the most significant signaling pathways. Gene expression in cytokine-cytokine receptor interaction pathways was mapped using PATHVIEW. Red: upregulated genes; green: downregulated genes.