Transcriptomic Profiling of Human Limbus-Derived Stromal/Mesenchymal Stem Cells—Novel Mechanistic Insights into the Pathways Involved in Corneal Wound Healing

Limbus-derived stromal/mesenchymal stem cells (LMSCs) are vital for corneal homeostasis and wound healing. However, despite multiple pre-clinical and clinical studies reporting the potency of LMSCs in avoiding inflammation and scarring during corneal wound healing, the molecular basis for the ability of LMSCs remains unknown. This study aimed to uncover the factors and pathways involved in LMSC-mediated corneal wound healing by employing RNA-Sequencing (RNA-Seq) in human LMSCs for the first time. We characterized the cultured LMSCs at the stages of initiation (LMSC−P0) and pure population (LMSC−P3) and subjected them to RNA-Seq to identify the differentially expressed genes (DEGs) in comparison to native limbus and cornea, and scleral tissues. Of the 28,000 genes detected, 7800 DEGs were subjected to pathway-specific enrichment Gene Ontology (GO) analysis. These DEGs were involved in Wnt, TGF-β signaling pathways, and 16 other biological processes, including apoptosis, cell motility, tissue remodeling, and stem cell maintenance, etc. Two hundred fifty-four genes were related to wound healing pathways. COL5A1 (11.81 ± 0.48) and TIMP1 (20.44 ± 0.94) genes were exclusively up-regulated in LMSC−P3. Our findings provide new insights involved in LMSC-mediated corneal wound healing.


Introduction
The cornea is the transparent and highly specialized tissue located in the anterior portion of the eye. In addition to its function as a protective barrier, the cornea is largely responsible for the transmission of light onto the retina, accounting for two-thirds of the eye's refractive power [1][2][3]. Anatomically, the cornea is made of three major layers: epithelium, stroma, and endothelium. The epithelium is a 4-6 layered outermost structure made of non-keratinized stratified squamous cells. Stroma is the middle layer comprising 90% of the corneal thickness and contributing to most of the structural framework. It is made of an extensive network of collagen fibrils with interstitially embedded cells called keratocytes, and proteoglycans such as lumican, keratocan, and decorin. The stroma is followed by endothelium, the innermost layer. Endothelium is majorly responsible for the ( Figure 1A) and late stages ( Figure 1B,C) of the primary culture (>80% confluence, days [8][9][10]. Further subcultures were observed to show the elevated number of stromal cells and the gradual decrease in the epithelial cells from passages P1 ( Figure 1D) to passages P2 ( Figure 1E). A pure population of the limbal stromal cells without any presence of the epithelial cells was derived in passage P3. This population of stromal cells also featured the presence of few myofibroblasts in their undifferentiated state, ( Figure 1F; orange arrow), dendritic cells ( Figure 1F; black arrow), and few quiescent fibroblastic cells ( Figure 1F; white arrow), similar to the native corneal stroma.

Stem Cell and Ocular Biomarkers
The immunostaining analysis has revealed a similar pattern of the expression (positive) of the stem cell (ABCG2, p63-α) and ocular biomarkers (PAX6) at both LMSC−P0 and LMSC−P3 stages of the culture. However, ABCB5 was found to be expressed in a low number of cells in P0 relative to P3. Additionally, the number of cells positive for p63-α were high in P0 relative to P3 (Figure 2). ABCB5 plays a vital role in the differentiation of limbal stem cells and is essential for corneal repair [34].
However, the RT-PCR data have revealed significant up-regulation of ABCG2 in LMSC−P3 with respect to that of LMSC−P0 and native limbus ( Figure 2). ABCB5 was found to be significantly down-regulated in both LMSC−P0 and LMSC−P3 and PAX6 was found to be significantly up-regulated relative to limbal tissue. P63α was found to be down-regulated by 3-fold in LMSC−P3 relative to native limbal tissue. On the contrary, the level of P63α was found to be up-regulated 3-fold in LMSC−P0 relative to the control. , and PAX6 genes quantified using qRT-PCR in limbal epithelial (LMSC−P0) and stromal (LMSC−P3) cells, relative to native limbal tissue (n = 5). P63α was found to be down-regulated in LMSC−P3 where all the other stem cell genes ABCG2, ABCB5 and PAX6 were found to follow same pattern of expression in both early and late passages of the culture. The results were plotted as mean log 2-fold change ± SD. The statistical analysis was performed using Kruskal-Wallis one-way ANOVA test. # p > 0.05, ** p < 0.01, *** p < 0.001.

Mesenchymal Stem Cell Markers
The limbal stem cells at LMSC−P0 were found to be positive in relatively low numbers for the mesenchymal stem cell (MSC) biomarkers CD90, CD105, and VIM (Vimentin). However, most of the cells at LMSC−P3 were found to be positive for the above markers ( Figure 3). The qRT-PCR analysis revealed a down-regulated expression of markers CD90, VIM (1-fold), and CD105 (6-fold) in LMSC−P0 relative to the control, limbal tissue. On the contrary, the levels were found to be up-regulated in LMSC−P3 by 2-fold of CD105, 3-fold of VIM, and~20-fold of CD90 relative to the control ( Figure 3).
The transmembrane proteins NCAD (N-cadherin) and ECAD (E-cadherin) were observed to express positively in LMSC−P0 through the immunostaining analysis. However, Ecad was found to be negative in LMSC−P3, and Ncad showed positive expression. The qRT-PCR analysis showed that NCAD was found to be up-regulated in both LMSC−P0 (8-fold) and LMSC−P3 (26-fold) compared to the native limbus. ECAD was found to be down-regulated in LMSC−P3 (3-fold), while it was found to have an increased expression in LMSC−P0 (3-fold) relative to the control ( Figure 3). Ncad (N-cadherin) were positive (red) in LMSC−P3 cells and Ecad (E-cadherin) did not show any expression in LMSC−P3. Level of expression of VIM, CD90, CD105, NCAD, ECAD genes quantified using qRT-PCR in LMSC−P0 and LMSC−P3 relative to native limbal tissue (n = 5). Except ECAD remaining genes were found to be up-regulated in LMSC−P3 with fold-change ranging between 2 to 20, which were down-regulated in LMSC−P0. Scale: 50 µm. The results were plotted as mean log 2-fold change ± SD. The statistical analysis was performed using Kruskal-Wallis one-way ANOVA test. * p < 0.05, *** p < 0.001.

Transcriptome Overview Using Principal Components Analysis Plot
To visualize the overall similarities or differences between gene expression patterns in different cell types, the counts were analyzed through Principal Component analysis (PCA). The counts data was subjected to Box-Cox transformation to stabilize the skewness in the data before PCA analysis. This analysis has showed the overall differences in the expression patterns of the samples in terms of the distances between them, which indicates the similarity between their expression profiles. It was found that sclera and cornea clustered together and are quite distant from the other samples. This indicates that the differences in their gene expression is not as heterogeneous ( Figure 4A) as compared to the rest of the analytes (Limbus, LMSC−P0, LMSC−P3 and ESC (embryonic stem cell)). The count data from all the samples were transformed using Box-Cox transform to compensate for skewness before PCA analysis. The closer proximity of the samples indicates similarity of the expression profiles of those samples. Sclera and cornea were found to show high similarity in whole transcriptome of expression and were clustered together. Similarly, LMSC−P0 and native limbal tissue were in close proximity, indicating similar transcriptomic signature. LMSC−P3 and ESC were found to be further away from one another, indicating an altered or different expression profile relative to the rest of the analytes. (B) The heat map representing the DEGs in all 5 samples relative to the control (scleral tissue). The rows indicate the genes and columns indicate the samples (cells or tissues). The color intensity represents the level of changes in expression. All significantly up-regulated genes are indicated in green and all significantly down-regulated genes are indicated in red. p < 0.05 was considered to be a statistically significant change in the gene expression. The limbus tissue and LMSC−P0 were observed to form an isolated cluster away from LMSC−P3 and ESC. The altered transcriptomic signature of LMSC−P3 may possibly be the result of repeated passaging and de-epithelialization. However, it was not very distinct from ESC, indicating possible shared/similar gene expression patterns such as pluripotent nature and dedifferentiation.

Visualizing the Asymmetry in Gene Expression of Various Tissues
Around 28,000 genes were detected via RNA-Seq in all the analytes. Among them, 7800 genes were differentially expressed (either up-regulated or down-regulated) against scleral tissue as a control ( Figure 4B). In limbal tissue, a total of 1036 genes were upregulated and 1093 genes down-regulated. LMSC−P0 had 1570 genes up-regulated and 1838 genes down-regulated, wherein LMSC−P3 774 and 1530 genes were up-regulated and down-regulated, respectively.
The asymmetry in gene expression by cornea, limbus, LMSC−P0, LMSC−P3, and ESC was visualized by plotting their transcriptome through volcanic plots. Volcanic plots provide a visual representation of the DEGs, showing their statistical significance (p values) versus the magnitude of change (fold-change). These scattered plots have shown that the transcriptome of corneal tissue had a major proportion of the down-regulated genes ( Figure 4C). On the other hand, limbal tissue had a distinct asymmetry, with a major proportion of the genes significantly up-regulated ( Figure 4D). LMSC−P0 showed a near symmetry in the plot ( Figure 4E), while LMSC−P3 has a smaller proportion of up-regulated genes ( Figure 4F). The ESC had shown a large number of down-regulated genes ( Figure 4G).

Tissue-Specific Differential Expression and Pathway Enrichment Analysis
The differential expression of the genes that were either specific to a particular type of cell or tissue was analyzed using 5-way Venn diagrams ( Figure 4H,I). The information on the number of genes commonly expressed in one or more cells/tissues was obtained. In addition, the number of genes that were either exclusively up-regulated or exclusively down-regulated in one particular type of cell/tissue was also obtained. The number of genes exclusively up-regulated in LMSC−P0 and LMSC−P3 was 459 and 223, respectively, while the exclusively down-regulated ones were 465 and 387, respectively ( Figure 4H,I).
Pathway-specific Gene Ontology (GO) enrichment analysis using the Enrichr tool provided insights into various cellular processes and pathways such as apoptosis, cell motility etc., where one or two particular cells/tissues were playing a major role. This was evident from the statistically significant gene expression with respect to the control (sclera). The relative median change indicated the up-regulation or down-regulation of such genes with respect to the basal expression levels of the whole transcriptome ( Figure 5A). Few of the prominent or majorly observed cellular processes were plotted against the relative expressions of DEGs specific to each of these processes. Figure 5. Interpretations from Gene Ontology enrichment analysis. (A) Gene ontology pathwayspecific analysis: Each row represents the genes belonging to a particular pathway/biological process in the GO database. The dot color/size represents the difference in median expression between genes of a particular pathway and rest of genes in whole transcriptome. A positive value indicates up-regulation (blue) and a negative value indicates down-regulation (red). The difference was tested using Mann-Whitney-Wilcoxon test and statistically insignificant ones were denoted with crosses. Inflammatory response is down-regulated more strongly in P3 versus P0, which may reflect why the use of P3 stage cells does not cause fibrosis in corneal stromal transplants. The stronger downregulation in ESC may reflect on the immune privilege of embryonic stem cells. The cell motility pathways are active in P3, cornea, and ESC. This is easily explained for the stromal stem cells in P3 and the ESC in terms of their proliferation and migration activity before differentiation, but for cornea may be representative of continuous cell migration required to replace lost corneal tissue. (B) The heatmap of 254 genes belonging to wound healing pathway. The rows indicate the genes and columns indicate the samples (cells or tissues). The color intensity represents the level of changes in expression. All significantly up-regulated genes are indicated in green and all significantly downregulated genes are indicated in red. p < 0.05 was considered to be statistically significant change in the gene expression. (C,D) Genes of wound healing pathway exclusively expressed in LMSC−P3: The levels of COL5A1 and TIMP1 in LMSC−P0 and LMSC−P3 assessed through RNA-Seq (C) and qRT-PCR (D). The LMSC−P3 has~10-fold (n = 3. p < 0.001) high expression of COL5A1, which was down-regulated in LMSC−P0, evident from both techniques. The levels of TIMP1 were~6-fold high (n = 3, p < 0.001) and 4-fold higher (n = 5, p < 0.01) in LMSC−P3, when assessed through RNA-Seq and qRT-PCR respectively. ** p < 0.01, *** p < 0.001. The results were plotted as mean log 2-fold change ± SD. The statistical analysis was performed using Kruskal-Wallis one-way ANOVA test for qRT-PCR and two-tailed T test for RNA-Seq analysis.

Interpretations from Gene Ontology Enrichment Analysis
A total of 6634 unique genes belonging to 16 relevant biological processes were found to be expressed by the corneal and limbal tissues and cells LMSC−P3, LMSC−P0 and ESC. The relative comparison of DEGs specific to various cell processes expressed by the cells of interest in this study-LMSC−P3 and LMSC−P0 cells-has provided interesting results.

Other Signaling Pathways
The AmiGO gene ontology analysis of the total DEGs has revealed that 211 genes playing a role in the Wnt signaling pathway were differentially expressed by cornea, limbus, ESC, LMSC−P3, and LMSC−P0. The relative expression levels of these genes from the RNA-Seq by each cell/tissue were plotted in Supplementary Figure S3B, tabulated in Supplementary Table S6. A total of 85genes belonging to the TGF-β signaling pathway were found to be differentially expressed by one or more cells/tissues. Among them, COL3A1 was found to be exclusively up-regulated in LMSC−P3 alone.
Of the 23 genes available in the GO database which belong to stem cell pathway, 13 DEGs (Supplementary Table S6) were found to be expressed by the cells or tissues analyzed in this study (GO consortium accession CL 00000034). The plot of their relative expression levels has shown (Supplementary Figure S3A) that the majority of genes were found to be significantly down-regulated in LMSC−P0.

Discussion
Various groups across the globe have attempted to understand the basic biology and the mechanisms involved in the healing or repair of the corneal wounds resulting from trauma [49][50][51] or the regeneration of the corneal epithelium and stroma lost due to the regular wear and tear [30,52]. The epithelial homeostasis is achieved primarily through the LESCs residing in the limbal crypts [53], in which the systematic synthesis and the degradation of the collagens in the stromal extracellular matrix (ECM) released by the native keratocytes helps in the maintenance of the corneal stromal integrity and homeostasis [54]. The interactions between the epithelial and stromal cells affect the repair of the cornea after an injury. Earlier studies have shown that the communication or the interaction between the LESCs and the stromal cells through their cytokines and other secretory molecules is essential for maintaining the corneal integrity [28,37,53] and thereby its transparency. IL-1 and its isoforms (IL-1α and IL-1β), produced by the epithelial cells during corneal injury, promote the production of TNF-α, KGF, and HGF [55,56]. Together with TNF-α, IL-1 also modulates the production of growth factors (PDGF and family) that modulate the chemotaxis and proliferation of corneal fibroblasts [57]. They also enhance the levels of cytokines such as G-CSF, neutrophil-activating peptide, IL-3 precursor, IL-4, IL-6, IL-7, IL-8, IL-9, and IL-17 [58]. These cytokines trigger the entry of inflammatory cells to the site of injury [59,60]. HGF and KGF released by the stromal fibroblasts, along with bFGF, IGF, and EGF, modulate the interactions between epithelial and stromal cells, regulating the migration and differentiation of damaged epithelial cells [61][62][63][64][65]. IL-6, a multifunctional cytokine, modulates the repair of the cornea in many ways. It enhances the epithelial wound closure, and low levels of IL-6 delays the healing [66][67][68]. Additionally, it reduces the levels of IL-1 and TNF-α, lowering inflammation [69]. A study by Samaeekia et al. [37] has shown that the exosomes isolated from the corneal and peripheral limbal MSCs enhance the migration and proliferation of corneal epithelial cells in vitro. The co-culture of corneal epithelial cells and corneal stromal cells has been shown to reduce the levels of pro-inflammatory cytokines and enhance the number of viable epithelial cells following an injury [70].
However, in the cases of corneal injuries (limited to the layers of epithelium and its surface, as well as the stroma), the healing process that follows involves one or more factors such as the native cells, growth factors, genes, cytokines, and antigen-presenting cells and even lipids [71][72][73][74][75]. The healing/repair process could involve just one of the above factors or a cascade of multiple events and reactions based on the site/location and the severity of the wound. Additionally, not all of them can be favorable towards the transparency of cornea. Mechanisms such as corneal fibrosis result in opaque/scarred cornea obscuring the visual pathway. The LMSCs were proven to be one of the promising intervention which could prevent and repair the corneal wound without needing a whole corneal replacement [22,38]. These cells are capable of differentiating into the native keratocyte phenotype [22,23]. Recent studies by Orozco Morales et al. [70], Hertsenberg et al. [76], Weng et al. [77] Chameettachal et al. [78] and Chandru et al. [79] have shown the potential of these cells in healing the cornea both in vitro and in vivo in animal models. However, the underlying mechanisms of how these cells achieve the scarless wound healing is not clearly studied. The current study aimed in uncovering the pathways and genes or other factors involved in the corneal wound repair by the LMSCs.
The LMSCs were isolated from cadaveric donor corneo-limbal rims and cultivated in a GMP-certified clean-room facility. Cells at the primary cultures (P0) where both mesenchymal/stromal stem cells of limbus and LESCs were obtained and cell population at the third passage where a pure population of the limbal mesenchymal/stromal cells were obtained (Figure 1), were subjected to RNA sequencing and immunostaining analysis. The outcomes of these two methods were further validated through the qRT-PCR. The mix population of the cells at primary culture were termed LMSC−P0 and the latter was termed LMSC−P3. The digestion of limbal tissue with collagenase alone and maintenance of low serum levels may possibly have led to the propagation of limbal mesenchymal/stromal cells only. The complete removal of serum may lead to the generation of fibroblastic cells with reduced keratocyte phenotype [80]. Conversely, a low quantity of serum (2%) [22] after digestion with collagenase alone [81] would allow stromal cells to proliferate with gradual loss of epithelial islands in the culture. Cells in both populations were found to express the stem cell ocular biomarkers positively ( Figure 2). However, the number of cells positive for collagens and mesenchymal biomarkers was more in LMSC−P3 (Figure 3). The collagens of corneal stromal ECM also followed a similar trend, with more expression in LMSC−P3 (Figure 4).
The principle component analysis plot revealed an altered transcriptomic signature of the LMSC−P3 from the rest of the clusters. Of the 28,000 genes detected, nearly 7800 were found to DEGs from all the samples, with LMSC−P0 having more number of DEGs. The asymmetry of the up-regulated or down-regulated genes visualized through volcanic plots revealed a near symmetry in LMSC−P0 ( Figure 5). The gene ontology enrichment revealed 6344 unique genes with functions in more than 16 biological processes ( Figure 6 and Supplementary Table S4). Genes belonging to signaling pathways such as Wnt (211 DEGs), TGF-β (85 DEGs), stem cell (23 DEGs), and wound healing pathways (254 DEGs) were also obtained (Supplementary Figure S1, and Supplementary Tables S5 and S6). Many studies have proven the anti-inflammatory and immunomodulatory properties of the LM-SCs [70,76,77,82]. The findings of the current study also support the anti-inflammatory nature of these cells. The overall genes of the inflammatory response (734) were downregulated in LMSC−P3 relative to LMSC−P0 ( Figure 5A). The pro-fibrotic gene IL-13 (Figure 7), and inflammatory genes C3 and CXCL8 which may lead to corneal neovascular-ization, etc., were down-regulated in LMSC−P3. Additionally, the anti-inflammatory gene 1L-10 ( Figure 7) was up-regulated in LMSC−P3 relative to LMSC−P0.
COL5A1 is a prominent and vital regulator of fibrillogenesis [83], the levels of which were reported to be high during the healing of scars [84,85]. During wound healing, the fibroblasts recruited to the site of injury produce collagens type I and V for extracellular matrix regeneration and restoration of the corneal thickness. In our study, the levels of COL5A1 were found to be higher in LMSC−P3 when compared to LMSC−P0 and native cornea. A similar finding was reported by Ruggiero et al. [86], who have shown that the amount of type V collagen produced by corneal fibroblasts in vitro is higher than that of the native cornea. Moreover, studies by McLaughlin et al. [87] and DeNigris et al. [88] have reported that the altered fibroblasts affect the level of collagen V in vitro. This also justifies and explains the levels of COL5A1 being proportionate to the number of fibroblasts in cells/tissues analyzed, i.e., LMSC−P3, followed by LMSC−P0, cornea, and limbus (Supplementary Table S5 and Figure S5C). The number of fibroblasts was also relatively high in LMSC−P3 compared to LMSC−P0 ( Figure 1A-C,F). These findings were similar to the study by Z.H. Guo et al. [89], who provided insights into the molecular mechanisms of differentiation and stemness maintenance by limbal stem cell niche in mice. The collagen genes of corneal stroma are responsible for collagen synthesis, which is predominantly regulated by COL5A1 [90]. The exclusive up-regulation of the COL5A1 by the LMSC−P3 cells evidently shows their ability and makes them an ideal source for repair and regeneration of corneal tissue through collagen fibrillogenesis ( Figure 5C,D).
The other gene exclusively up-regulated in LMSC−P3 was TIMP-1, an inhibitor of the matrix metalloproteinases (MMPs), the genes responsible for cleaving collagens. The tissue inhibitor of metalloproteinases (TIMPs), inhibit these MMP genes, highly regulating the corneal ECM. The binding of TIMPs to the MMPs prevents the degradation of the ECM. TIMP-1 inhibits all active MMPs, except membrane type matrixins (MT1-MMP), whereas TIMP-2 inhibits MMP-2, in particular [91,92]. These two groups of genes, i.e., the MMP family and the TIMP family, also a play vital role in the development of cornea [93]. In this study, we have found that TIMP-1, MMP-1, and MMP-9 were found to be up-regulated and that TIMP-2 and MMP-3 were down-regulated in LMSC−P3. A similar trend was observed in the LMSC−P0, except for the levels of MMP-3. Although MMP-9 in LMSC−P3 was up-regulated, the levels of TIMP-1 were much higher in terms of fold-change. Unlike the earlier studies [94,95], the positive correlation between the levels of TGF-β and TIMP1 was also not observed in our study ( Figures 5C,D and 7), which did not involve disease condition or the altering of their concentration in culture. Assessing all these genes in a disease condition may provide a better understanding of their respective roles in corneal regeneration.
The exclusively elevated genes on LMSC−P3 interact through various genes and biological processes. The network functional enrichment analysis performed to understand their interactions has revealed a set of interleukins, matrixins, chemokine receptors, and growth factors. Most of these were up-regulated in LMSC−P3 relative to LMSC−P0 (Figure 7). Corneal ECM genes such as Keratocan, Lumican, and SMA were expressed significantly higher in LMSC−P3 relative to LMSC−P0 and native limbus. Lumican and keratocan belong to the SLRP (small leucine-rich proteoglycan) family, which is critical for corneal clarity. They are responsible for the fibrillar organization of the collagens in the ECM of the corneal stroma [96,97]. Both these proteoglycans play a crucial role in corneal wound healing and regulate inflammation by localizing the macrophages to the site of injury and recruiting neutrophils [97]. The levels of lumican and keratocan were reported to decrease during the scarring of cornea [98]. Unlike the studies [99,100] that reported low expression of keratocan by keratocytes in vitro, we observed relatively high levels of keratocan in LMSC−P3. However, when compared to LMSC−P0, where there is no chance of differentiating the expression of keratocyte markers by a diverse set of cell populations and the relatively less number of stromal cells, the high number of stromal cells in LMSC−P3 could attribute to the high levels of keratocan and lumican. The down-regulation of TGF-β could also be attributed to the keratocan levels, as shown by Kawakita et al. [101], with decreasing levels of TGF-β maintaining the levels of keratocan. This indicates the strong keratocyte-like nature of the cells in LMSC−P3 with respect to LMSC−P0. The increased expression of SMA in LMSC−P3 relative to LMSC−P0 could be attributed to the relatively high number of myofibroblastic cells in LMSC−P3 than in LMSC−P0.
We have also found that the expression of VEGFA, a proangiogenic factor, was also significantly high in LMSC−P3. The continuous maintenance of corneal avascularity is important for optimal visual acuity. Angiogenesis is one of the many vital processes in wound healing for the successful repair of damaged tissue. The balance between the proangiogenic and anti-angiogenic factors is mandatory for maintaining corneal avascularity [102]. To assess the levels of angiogeneic factors that regulate the formation of vasculature on corneal surface, certain genes were quantified through qRT-PCR. VEGFA is one of the proangiogenic factors which, besides FGF-2 [102], plays a role in multiple processes such as immune-modulation, epithelialization, collagen deposition, and cell migration [103]. It decreases the duration of wound healing [104]. Although MSCs were reported to potentially lower angiogenesis [105], the surprisingly high expression of the VEGFA in LMSC−P3 (Figure 7) is questionable due to the fact that elevated vasculature over the surface of cornea can potentially affect the visual acuity [106,107]. However, the elevated levels of VEGFA (growth factor-induced or transfected or topically applied) in the wound bed were reported to enhance the wound repair of dermis/skin [108][109][110][111], but not many studies on corneal surface were reported. These elevated levels of VEGFA also contradict decreased expression of MMP9, the factors reported to feedback regulation mechanism [112]. However, other proangiogenic factors such as PDGF and its family (PDGFB, PDGFC, PDGFD, PDGFRA, and PDGFRB) are either unexpressed or down-regulated in LMSC−P3 (Supplementary Table S5). Although native corneal epithelial tissue is reported to have detectable levels of VEGFA and sflt-1 [113,114], not much information is available regarding the levels of VEGFA in native limbal tissue. However, the levels of VEGF expression occurs differently in different cells in vitro. The limbal epithelial cells were earlier reported [115] to be anti-angiogenic in nature and the limbal fibroblasts proangiogenic in nature. The corresponding high and low levels of the limbal fibroblasts in LMSC−P3 and LMSC−P0, respectively, could possibly explain the increased levels of VEGFA. However, a contradicting observation was reported much later in a study by Eslani et al. [116], who have shown that the LMSCs are anti-angiogenic. Low levels of VEGFA and high levels of the anti-angiogeneic factors SFLT-1 and PDGF were observed in the secretome of LMSCs. In the current study, determining the levels of SFLT-1, MMP-2, MMP-14, and CTGF genes in the cell populations/tissues tested could have provided a better answer to this conundrum. Further studies to explore/evaluate the levels of VEGF in a corneal wound model treated with LMSCs and monitoring of the progress of healing may be required.

Ethics Approval and Tissue Collection
Human donor corneas (donor age ranged between 18-60 years) were collected from the Ramayamma International Eye Bank (RIEB), LV Prasad Eye Institute, Hyderabad. Overall, 21 therapeutic-grade donor corneas, unutilized for surgical purposes, were used in this study (n = 21). The corneas were collected with informed consent and in compliance with the guidelines of the Declaration of Helsinki. Ethical approval was obtained from the Institutional Ethics Committee (Reference number LEC 05-18-081) and the Institutional Committee for Stem Cell Research, of LV Prasad Eye Institute, Hyderabad, India (IC-SCR-Ref No: 08-18-002).

Establishment of Limbal Stem Cell Culture
The tissue processing was done using a stereomicroscope (SZX10, Olympus, Japan) to set up the limbal stromal stem cell culture, as described previously [117], and for total RNA extraction. Briefly, cadaveric corneas were washed with 1X PBS (14190250, Thermo Fisher Scientific, Waltham, MA, USA) fortified with 2× antibiotics (15240062, Thermo Fisher Scientific, Waltham, MA, USA) and were stripped of endothelium and iris. Full thickness limbus was excised in 1× PBS and then fragmented to small pieces in plain DMEM/F12 media (BE04-687F/U1, Lonza, Basel, Switzerland). These tissue fragments were minced for 1-2 min. The dissected limbal tissue is then enzymatically digested by Collagenase type IV (17104019, Thermo Fisher Scientific, Waltham, MA, USA) enzyme (200 IU per one donor rim), added to 1 mL of plain DMEM/F12 media and then incubated for 16 h. The digested tissue is sedimented twice at 1000 rpm/3 min in PBS. The pellet is then suspended in complete media, i.e., DMEM/F12 fortified with 2% Fetal Bovine Serum (SH30084.03, Cytiva Life Sciences, Shrewsbury, MA, USA) and supplemented with human recombinant Epidermal Growth Factor (PHG0311L, Thermo Fisher Scientific, Waltham, MA, USA) and human recombinant Insulin (12585014, Thermo Fisher Scientific, Waltham, MA, USA). This suspension was plated and cultured for 3 generations. Cells upon confluence, at the stages of primary culture (LMSC−P0) and passage 1 (LMSC−P1) and passage 3 (LMSC−P3), were used for analysis.

Immunofluorescence Assay
Cells were grown on the surface of glass coverslips in complete media until confluence. The cells were then fixed with 4% Formaldehyde (30525-89-4-500G, Fisher Scientific, Bangalore, India) for 10 min and washed twice with 1× PBS before permeabilization with 0.  Table S1) diluted in 1% BSA. This was followed by a wash with PBS thrice for 10 min and incubation with secondary antibody (Supplementary Table S1) (diluted in 1% BSA) for 45 min, which was further washed thrice and mounted onto a glass slide using Fluoroshield Mounting Medium with DAPI (ab104139, Abcam, San Francisco, CA, USA) for nuclei counterstaining. Staining of negative controls was done by omitting the primary antibody. Images were documented using an inverted fluorescence microscope (Axio Scope A1, Carl Zeiss AG, Oberkochen, Germany). The biomarker panel of the MSC phenotype was chosen in accordance with the minimal criteria set for multipotent mesenchymal stem cells [18].

RNA Isolation
Total RNA was isolated from tissues (sclera, limbus, and cornea) and limbal stem cells (LMSC−P0 and LMSC−P3) and embryonic stem cell line (ESC) using TRIzol™ reagent (15596018, Thermo Fisher Scientific, Waltham, MA, USA). The spent medium was removed from the 80% confluent cell culture. Cells were then washed with 1× PBS (prepared with DEPC-treated distilled water for RNA isolation) (AM9920, Thermo Fisher Scientific, Waltham, MA, USA), and an appropriate volume of TRIzol™ reagent was added to the cells. The cell lysate was mixed several times through a pipette and transferred to a sterile 1.5 mL micro-centrifuge tube. To the lysate, 0.5 mL of Chloroform (96764, Sisco Research Laboratories, Mumbai, India) was added per every 1mL of TRIzol™ reagent and incubated at room temperature for 15 min. This was followed by centrifugation at 12,000 rpm for 15 min at 4 • C. The aqueous phase was collected in a fresh tube and 1 mL of Isopropanol (Q13825, Thermo Fisher Scientific, Waltham, MA, USA) was added (in equal volumes with TRIzol™ reagent) and incubated at room temperature for 3 min followed by centrifugation at 12,000 rpm for 3 min at 4 • C. RNA pellet was washed with 75% Ethanol (24102, Sigma-Aldrich, St. Louis, MO, USA), air dried, and dissolved in 25 µL of nucleasefree water (AM7020) (volume dependent on size of RNA pellet). RNA was quantified by measuring the absorbance using a spectrophotometer along with the purity evaluation by the ratio of A260/280 (NanoVue™ Plus, 28956058, GE Healthcare Bio-Sciences AB, Chicago, IL, USA). Further confirmation was done through gel electrophoresis, using 1% agarose gel (50004, Lonza, Basel, Switzerland) stained with Ethidium Bromide (93079, Sisco Research Laboratories Private Limited, Mumbai, India). The RNA was treated with DNase I (AM2222, Thermo Fisher Scientific, Waltham, MA, USA) according to manufacturers' protocol. Briefly, a 30 µL reaction volume containing 30 µg of total cellular RNA, 1× reaction buffer, 6U of DNase I (RNase free), and nuclease-free water. The reaction mix was incubated at 37 • C for 30 min. After incubation, 70 µL DEPC water was added to the reaction mix and the RNA was purified by adding 100 µL TRIzol™ reagent. The RNA was quantified by measuring the absorbance using a spectrophotometer, as previously described, and 1µg each of the RNA from the analytes was used for the RNA-Seq study.

Next Generation RNA Sequencing (RNA-Seq) and Library Preparation
One microgram each of the total RNA from limbus, cornea, sclera, LMSC−P0, LMSC−P3, and embryonic stem cells (ESC) were subjected to RNA sequencing via Illumina platform using the reagents provided in the Illumina® TruSeq® Stranded Total RNA Sample Preparation Ribo-Zero™ kit (RS-122-2201, Illumina, San Diego, CA, USA). The first step involves the removal of ribosomal RNA using Ribo-Zero™ rRNA removal beads provided in the kit. The Ribo-Zero™ rRNA reagent depletes samples of cytoplasmic ribosomal RNA. Following purification, the RNA was fragmented into small pieces by heat digestion using divalent cations (magnesium or zinc) under elevated temperature. The cleaved RNA fragments were copied into first strand cDNA using reverse transcriptase and random primers. This is followed by second strand cDNA synthesis using DNA polymerase I and RNase H. These cDNA fragments then have the addition of a single 'A' base and subsequent ligation of the adapter. The products have been purified and enriched with PCR to create the final cDNA library. This sample preparation protocol provides the advantages of (i) strand information on RNA transcript and (ii) library capture of both coding RNA and multiple forms of non-coding RNA. The processed cDNA library of all 6 samples was used for paired end sequencing run (50 × 2 cycles) on the Illumina HiSeq 2500 platform (SY-401-2501, Illumina, San Diego, CA, USA).

Pre-Processing of the RNA-Seq Data for Data Analysis
The Fastq file was obtained from sequencer after trimming the adapter sequences using bcl2fastq program. Fastq data was used for alignment with the hg19 version of the human genome using the TopHat program with options provided as transcript annotation file. The alignment data has been used for guided transcript assembly using the Cufflinks program. After that, we merged transcripts across samples using the Cuffmerge program to make a reference transcript assembly. This merged transcript assembly has been used as a reference to compare for differential gene expression between a pair of samples with the use of Cuffdiff program. The resultant Cuffdiff output file has provided the normalized expression of genes/transcript in the form of counts, and the fold differences converted into log2 values. The details of the reference links of all the software/programs/bioinformatics tools used in analysis of the RNA-Seq data were provided in the Supplementary Table S3.

Differential Expression Analysis
The counts obtained for each sample were analyzed by using the EBSeq tool (Supplementary Table S3) for differential expression by considering scleral tissue as the control. A list of DEGs was obtained for the tissues with pairwise comparison to sclera, and multiple testing corrections were applied at a False Discovery Rate (FDR) of 0.05 percent. The heatmaps were generated using R software.

Delineating Cell-Specific Gene Expression Patterns and Testing for Pathway Enrichment
To delineate the DEGs in different tissues and cells according to their cellular specific expression, 5-way Venn diagrams were used to find the genes which were exclusively upregulated and down-regulated. Two types of the gene expression patterns were analyzed.
The cell type-specific gene expression to find genes that were exclusively differentially regulated in only one type of the cell/tissue and which may therefore serve as transcriptomic markers to identify the unique cell type. Pairwise overlapping genes that are differentially regulated only in two cell types may indicate a shared functionality between the two cell types. To obtain pathway-level insights into the significance of the exclusively differentially regulated genes, we have conducted pathway enrichment analysis through the Gene Ontology studies using Enrichr tool. This analysis indicates statistically significant groups of genes that are belong to various biological/signaling pathways.

Gene Ontology Pathway-Specific Gene Expression Changes
Using the ontology keywords derived from the pathway enrichments obtained in the Enrichr analysis, the lists of genes specific to the pathway-keywords were obtained from the gene ontology database using the AmiGO tool. These gene lists were used to examine the differences in the frequency distribution of fold-change values of pathway-specific genes as compared to that of all the other genes in the transcriptome. These differences in the distributions were tested for statistical significance using the nonparametric Mann-Whitney-Wilcoxon U test. The median difference between the distributions was used to detect the direction of the shift in expression value. The values of median shift of the pathways across different samples were plotted against the crossed out the values which were not statistically significant (p ≤ 0.05).

Reverse Transcriptase PCR
One microgram each of the RNA was reverse transcribed to cDNA using the Super-script™ III First-Strand Synthesis System (18080051, Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions.

qRT-PCR
Quantitative PCR was performed using 200 ng of cDNA in a final volume of 25 µL reaction mix (K0221, Maxima SYBR Green/ROX qPCR Master Mix (2X), Thermo Fisher Scientific, Waltham, MA, USA) and a 0.2 µM primer concentration. The reaction was carried out using Step One (Applied Biosystems, Life technologies) hardware and software. The reactions were run in triplicates. The gene expression data were normalized to control the variability in expression levels to the geometric mean of the housekeeping gene. The expression level of target genes was represented as a relative expression by using 2 −∆∆Ct formula and the graphs were plotted using their Log2 fold-change values. The primer sequences are listed in the Supplementary Table S2.

Statistical Analysis
All the experiments were repeated at least thrice with biological triplicates. Statistical analysis was performed using the Graphpad Prism 6 software. The tests employed were Student's two-tailed t-tests, Kruskal-Wallis test, and a nonparametric one-way ANOVA test with p values ≤ 0.05 to assess the statistical significance. The results are presented as the mean ± standard deviation.

Conclusions
In the current study, we report the genes, biological processes, and pathways involved in the limbal stromal/mesenchymal stem cell-mediated corneal wound healing by employing RNA-Sequencing in human LMSCs (LMSC−Passage-0 and LMSC−Passage-3), for the first time. Differential expression of the genes (7800) belonging to the following pathways, namely, apoptosis, cell motility, dedifferentiation, inflammatory response, stem cell maintenance, tissue remodeling, and wound healing pathways, etc., were found. The interactions between the DEGs exclusively up-regulated by the LMSC−P3 in the wound healing pathway (COL5A1, COL1A1 and TIMP1) have revealed the processes involved in tissue remodeling and repair (collagen fibril reorganization, collagen biosynthesis, regula-tion of the metallopeptidase activity etc.) and the cytokines and other key genes regulating these processes. However, this study is limited by the small sample size, and further comprehensive studies needed to explore and understand all the DEGs and their biological relevance in corneal wound healing. On the whole, the findings of this study provide a brief glimpse into the molecular basis of tissue repair, and the remodeling of the cornea by human LMSCs and the therapeutic potential of this.