A Keratin 12 Expression-Based Analysis of Stem-Precursor Cells and Differentiation in The Limbal-Corneal Epithelium Using Single-Cell RNA-Seq Data

1. Introduction

The corneal epithelium (CE) constitutes the first line of defense of the vision system against pathogen invasion and it is also a critical component of optical visual focusing. The mature CE consists of 5-6 layers of stratified non-keratinized epithelium. Its differentiation plan and homeostasis resemble that of the other outer ectoderm-derived stratified epithelium such as the epidermis and various mucosae. Briefly, the surface-exposed cells at the upper layer, continuously exfoliate and this cell loss, in turn, is balanced by the proliferation of cells within the basal layer of the strata. Studies over the last 4 decades have shown that the short-term needs for cell replacement, whether homeostatic or acute, are fulfilled by a rapidly proliferating and partially differentiated cell type. Because the capacity for cycles of replication of these cells is limited, they have been classically referred to as transient amplifying (TA) cells. The resupply of this short-term response cell population depends on a population of stem/precursor cells which are normally in a semi-quiescent or slow cycling state but can increase their proliferation as needed to replenish a diminished TA cell population. These cells tend to reside within specialized niches of the underlying stroma and to maintain an extended replicative capacity [1,2,3].

Another common feature of most ectoderm-derived stratified epithelia is the generalized expression of the tonofilament-forming acidic-basic cytokeratin pair Krt14 and Krt5. However, the epithelia are distinguished from each other by the expression of distinct complements of acidic and basic cytokeratins, likely needed to generate tonofilaments that best suit the specific functional needs of each lining at each specific stage of differentiation [4]. The CE is characterized by the expression of the Krt12, not identified in any other tissue and Krt3, which is also expressed in other stratified epithelia [5]

In addition to its cytokeratin specificity. the CE presents a unique feature that distinguishes it from all other epithelia. Owing to the avascular nature of the CE, the precursor stem cells are segregated to the vascularized outer annulus that separates the cornea from conjunctival lineages. All basal limbal cells are negative for the tissue specific cytokeratins till the border with the transparent corneal zone and become strongly positive, in an abrupt manner, as the cell migrate into the avascular cornea [6]. From a health perspective, this dual-domain arrangement represents a biological weakness. Damage to the thin limbal rim interrupts the supply of cells for the maintenance of TAs at the adjacent peripheral corneal zone and opens the door for the centripetal invasion of the pro-vascular conjunctival epithelium [7]. This colonization leads to neovascularization, and, because the CNJE performs poorly as a fluid permeability barrier and pathogen refractor, to stromal scarring and pathogen penetration. The resulting blinding state is referred to as limbal stem cell deficiency (LSCD). The correction of an LSCD state requires the transfer of a suitable population of cells through a transplant operation. The success of these transplants depends on the regenerative capacity of the transferred cells. Hence, an in-depth characterization of the precursor compartment and the stage of early differentiation are valuable tools for the design of best clinical reparative strategies.

Microarrays and real-time polymerase chain reaction have been used to identify molecular and functional features that distinguish the conjunctiva, limbal and corneal gene expression features underpinning each phenotype [8,9,10]. More recently, the advent of the single-cell RNA seq (scRNA) methodology has led to the separation of the CE into clusters of cells exhibiting similar GE profiles, and through this clustering to the assignment of genes as primarily belonging to specific limbal, corneal, and conjunctival, basal or suprabasal domains [11,12,13,14]. Bearing in mind the intimate relationship between the expression of Krt12 and the differentiation process we have now used scRNA-Seq data to examine the correlation between the levels of Krt12 expression and the expression of other genes along the limbus to central cornea differentiation and migration axis.

2. Methods

2.1. Tissue Procurement

Unidentifiable cadaver donor corneas unsuitable for human transplantation were received in Optisol from the NDRI within 48 h of donor death and immediately processed. The Institutional IRB has determined that unidentifiable donated cadaver tissue is not subjected to human subject regulations. The donor were Caucasians males of ages, 68, 71 and 65.

2.2. Cell Preparation

We implemented a cell adhesion protocol designed to generate a cell population that harvested most basal epithelial cells, while excluding or drastically minimizing the recovery of supra-basal epithelial cells. Two strips of (approx. 0.5 mm width, one consisting of the whole limbal width and similar section of the corneal periphery (LiPe sample) and the other a strip of the adjacent peripheral cornea (Pe sample) as depicted in Figure 1A, were dissected under stereoscopic visualization. Maximal care was talking to ensure that the LiPe segment excluded conjunctival epithelium. Central cornea (Co) and conjunctival (CNJ) samples were also collected. The expression of Krt12 and the putative presence of non-epithelial cells within the epithelium is schematically described in Figure 1B.

The biopsies were incubated with 2 mg/ml Dispase II (Fisher, Waltham, MA) made in a bicarbonate-free 1:1 mix of Dulbecco’s modified minimum essential medium and Ham F-12 (DMEM-F12; Fisher) for 18 h, at 4° C under a 60 tilting/min motion. At the end of this treatment, sheets of epithelial cells were either free-floating in the Dispase solution or were lightly attached to the stroma, from which they were released by gently prodding with the tip of a jeweler’s forceps. The sheets were incubated in a 0.25 % Trypsin (Fisher) solution at 37° C for 5 min (2 ml/sample), the Trypsin medium was admixed with 4 volumes of DMEM-F12 – 10 % fetal bovine serum(Atlanta Biologics, Flowery Branch, GA) and single-cell suspensions were generated by 40 passes through a 5 ml pipette. Cells were spun down in a clinical centrifuge and resuspended, in 2 ml of DMEM/F12, triturated again, filtered through a 70 µ filter, and cultured for 3 ½ h in a 25 ml culture flask. After visually confirming the presence of single attached cells displaying side blebs indicating incipient spreading, the flask was set horizontally for 3-4 min to allow full draining of unattached cells. The accumulated medium was fully aspirated and the lightly attached cells were released by gentle streaming from a 1 ml pipette tip. This protocol recovery about 1/3 rd of the cells in the suspension. Eighty percent of the cells we used for scRNA-Seq analysis. One full set of CNJ, LiPe, Pe, and Co samples (Exp. 1) were derived from a single donor and two additional LiPe samples (Exp. 2 and 3) from another two corneas. The remaining twenty percent of the samples were complemented with 1 µg/ml propidium iodide (PI) and the forward light scattering, a relative measure of overall cell size, of the PI-negative cells was used to determine the fraction of Li and Pe derived cells within the population (Figure 1 C).

2.3. Single Cell RNA-Seq Reading

Libraries were prepared using the 10x Genomics Chromium Controller (Pleasanton, CA) in conjunction with the single-cell 3′ v2 kit. Briefly, the cDNA synthesis, barcoding, and library preparation were carried out according to the manufacturer’s instructions. Libraries were sequenced on a Novaseq 6000 instrument. (Illumina, San Diego). Sample demultiplexing, barcode processing, and unique molecular identifiers (UMI) counting were performed by using the 10x Genomics pipeline Cell Ranger v.2.1.0 with default parameters. Specifically, for each library, raw reads were demultiplexed using the pipeline command ‘cell ranger mkfastq’ in conjunction with ‘bcl2fastq’ (v2.17.1.14, Illumina) to produce two fastq files: the read 1 file contains 26-bp reads, each consists of a cell barcode and a unique molecule identifier (UMI), and the read 2 file contains 96-bp reads including cDNA sequences. Reads then were aligned to the human reference genome (GRCh38), filtered, and counted using ‘cell ranger count’ to generate the gene-barcode matrix. CSV files were derived from the feature-barcode matrix data following the 10x Genomics instructions. The CSV files were opened in Microsoft Excel and processed as described below.

3. Results

3.1. Data Pre-Processing

Figure 1C shows the flow cytometry FSC distributions for all four cell preparations submitted for scRNA-Seq. Limbal basal cells are remarkably and uniformly smaller than the adjacent peripheral cell and that the size transition occurs sharply at the Li-Pe interphase [15]. Basal limbal cells accounted for 53 % of the cells in the sample. The percent for the other two LiPe samples were 62% and 58 %, respectively. PI permeable cells amounted to less than 2 % of in all cases.

As a first step, since we seek to analyze only the gene expression of the bona fide epithelial cells within the overall population, we excluded any column containing markers for nonepithelial cells. Cell expressing any melanocyte markers (DCT, Tyr, Tyrp1, and/or Mlan-A), amounted to 6.1 % of the cells in the LiPe samples (Table 1). About 20 % of the melanocyte containing markers occurred in columns concurrently expressing keratin epithelial markers (Krt5; Krt14 or Krt 12). Hence it is reasonable to assume that these columns consist of epithelial-melanocyte hetero-cellular duplets. Consistent with the expected localization of melanocytes in the limbus, there were few cells expressing melanocyte markers in the pure peripheral or central corneal samples (Table 1). Finally, since the results from the CNJE sample and a previous study [16] indicated that this AQP5 is abundantly expressed in the CNJE and not expressed in the cornea, a few cells expressing this gene were also excluded. After the gene ID-based exclusions, all remaining cells expressed at least one of the two cytokeratins, Krt5 and Krt14, which are expressed in all stratified epithelia. The starting cell number, the number and percent of purged cells, and the number of cells remaining after the cell exclusion for each sample are detailed in Table 1.

Next, the mean gene count (MGC) per cell was calculated and plotted in decreasing order (Figure 2A). The resulting curves easily revealed an inflection at around an MGC of 0.5 in all cases. This accelerated count decrease is suggestive of cells in a decaying state. Thus, the analysis was limited to cells with MGCs above this value. In addition, we excluded all cells with an MGC larger than 2.5, as those are likely to include a large percentage of cell duplets. Post-trimming images are shown in a comparative format in Figure 2B. The columns were then normalized to an MGC of 1. All the calculations were performed in a ThinkStation computer controlled by an i9 8 core processor.

3.2. Domain Identification

The expression normalized spreadsheets remaining after the exclusion of cells from non-CE lineage and outlier cells were first used to examine the distribution of Krt12. Figure 2 C depicts the range of Krt12 expression in the LiPe, Pe, and Co samples for Exp. 1. Figure 1D incorporates the ranges for three LiPe samples each derived from a different donor cornea and processed at different times. The high similarity of the Krt12 plots for all three LiPe samples demonstrates the reproducibility of the scRNA-Seq protocol. Plots of all three LiPe samples in the lowest 0.5 % of the full expression range (0 – 550 A.U.s) highlighted the presence of 5 distinct domains of Krt12 increase,

D0 through D4. Figure 2E shows the distribution for the larger LiPe sample, LiPe-1. In fact, the pattern of Krt12 increase and the expression ranges were remarkably similar in all three LiPe samples, allowing their amalgamation into a single file. As described below this amalgamation was essential to augment the number of differentially expressed genes (DEGs). The Krt12 changes for the amalgamated files are depicted in Figure 2F. The Krt12 range and the cell number included in each domain are displayed in the inserted table.

In contrast with the Krt12 distribution in the LiPe samples, where 32 % of the cells did not expressed Krt12, the Pe and Co samples contained no Krt12-negative cells, and their Krt12 plots displayed a rapid raise starting from the lowest values (Figure 1C). In the Pe sample, only 20 cells or about 1 % of the total population, expressed Krt12 matching the expression range of D0-D3, 0-5.32 arbitrary units (A.U.) (Figure 2E). Thus, it can be reasonably concluded that the overwhelming majority of cells within the D0-D3 range derive from the limbus. They accounted for 60 % of the cells in the amalgamated sample (4003 out of 6667), roughly consisting with the percent of Li cells observed in the flow cytometry plot (Figure 1).

3.3. Gene Expression Differential Analysis

Differentially expressed genes (DEGs) were defined as genes with Benjamini–Hoechsberg-adjusted p-values (BHp values) lower than 0.01 in an expression comparison. Due to the high frequency of cells showing nil expression for the low expression genes, the identification of bona fide statistically significant genes decreased with the decrease in gene expression level. This is exemplified in Figure 3 where we examined the relationship between gene expression and the frequency of BHp < 0.01 values in a D0 to D4 comparison. Likewise, the size of the populations being compared can be expected to have a strong influence in the calculation of p values. This expectation was confirmed by comparing the BHp yields when only half of the D0 and D4 domains (e.g., by excluding even columns) was compared with the yield when using the full population sets; the list of identified genes organized by decreasing D0/D4 expression ratio were nearly identical, but BHp values were lower by more than 2 orders of magnitude, setting a large number of genes with BHp values outside our BHp > 0.01 DEG limit. Hence, comparisons were made using the amalgamated file and were limited to the highest 6647 expressing genes out of 21300 genes with a finite positive expression. The D4 domain we split into a D4i domain including the first 2100 cells and a D4ex domain comprising the last 564 cells. In the domain comparisons, it is likely that at the domains interfaces cells from both domains will intermingle. Hence, to sharpen the differential analysis we excluded 10 columns from each side of the D0-D1, D1-D23, and D0-D123 comparisons. Additionally, for the D123-D4i.1 comparison we excluded the first 100 columns from D4i and defined the next 1000 cells as D4i.1. The remaining 1000 cells of D4i were defined as Dri.2. This domain trimming is indicated in Figure 2F by the gaps in the Krt12 line. The D4ex subdomain (Figure 2F, table inset) was not used in the analyses. The Figure 2F inset provides an account of the range of Krt12 expression in A.U.s, the start and end of each domain in Supplemental Table S1 and the number of cells present in each domain (#). Table 2A provides an account of the identified DEGs for various domain comparisons. The D0-D1 comparison is particularly intriguing because it represents the expression changes associated with the initial gene expression of Krt12 that occur while cells are located within the limbal domain; it identifies the genes associated with the cell cohort enriched in the lineage stem cells. There were nearly twice as many genes down regulated genes as upregulated ones. Table 3 lists the top down regulated (D0/D1; R >1) and up regulated (D1/D0; 1/R >1) DEGs, sorted according to decreasing R or 1/R, respectively. The use of the signal ratios seems more relevant than the use of the BHp values, as discussed above these values are heavily influenced by the level of gene expression.

The next sequential comparisons, D1-D2 and D2-D3 yielded no DEGs. The combination of D2 and D3, to form the D23 domain, where the maximal Krt12 is three time as large as the maximal Krt12 in D1, yielded only 8 upregulated DEGs exceeding the 1/R >1.25x threshold. One of these genes, CRTAC1 was present within the upregulated D0-D1 DEG list. But CRTAC1 and all other seven genes, namely, TGFBI, CPXM2, CLU, IGFBP7, ALDH1A1, IGFBP2 and FTH1 were within the top genes in the upregulated D123-Di4.1 (see below) list, suggesting that the change in Krt12 rate of increase that establishes the D2 and/or D3 domains represent the earliest cell changes associated with the transition towards the gene composition of the corneal periphery cell. Comparing D0 against the full set of limbal domain Krt12-positive cells (D123) increased the >1.25x DEG number by about 30-40 % from the number of DEGs yielded by the D0-D1 comparison (Table 2). However, the average D1/D23 ratio for both down and up regulated DEGs yielded statistically identical ratios (Table 2C). Furthermore, 81 % of BHp values were lower in the D0-D123 set than in the D0-D1 set. Thus, notwithstanding the eight genes mentioned above, most of the increase in the DEG list when using D123 instead of D1 in the comparison with D0 is likely due to the effect of the sample size described above in the + p values.

Rationally, the next sequential comparison should be between D3 and the first Di4 subdomain Di4.1. However, given the minimal differences between D1, D2, and D3 and considering the small size of D3 and the strong effect of the population size in the chances for identifying DEGs, comparing the whole D123 with Di4.1 instead seemed a more efficient way for DEG identification. Consistent with the concept of a sudden gene expression pattern change at the Li to Pe transition this comparison yielded very large down and up regulated DEGs (Table 2). More than 40 % of the total cell population displayed an BHp > 0.01. Intriguingly, why differentiation is intuitively associated with new gene expressions, the same 1.5 excess of down to up regulated genes seen in the D0 to D1 transition was also a feature of the D123-D4i.1 DEG set. The top up and down DEGs, sorted according to descending R or 1/R values, respectively are presented in Table 4.

The comparison between D4i.1 and D4i.2 identified genes that undergo statistically significant change along with the increase in Krt12 GE. Table 5 lists the top genes of each category. The average expression values for Krt12 for the D4i.1 and D4i.2 were 37.1 and 127.3 A.U., respectively, equivalent to a 1/R of 3.43. Using this value as a basis for comparison only 3 genes displayed an upregulation of larger magnitude than Krt12. In contrast, 25 genes down regulated at a faster speed than the rate of Krt12 increase. Thus, down regulation is likely to be as critical for the drastic phenotype change through all the basal cell differentiation process. Having identified the DEGs associated with each transition, it was possible to identify DEGs that either increase or decrease in expression at each transition and those that display a differential expression on y at one of the

Krt12-defined domains. There were 151 genes that exhibit down regulation in each of the Krt12 transitions and 72 genes that exhibit continuous up re regulation. The top gene of these two categories are listed in the upper panels of Table 6. The ratio between the two extremes of our comparative analyses was used to sort genes in decreasing ratio. Most of the continuously up regulated genes are well known genes contributing to the corneal phenotype.

In addition, there were 179 genes that undergo down regulation only at the start of Krt12 expression in D1. These genes are strong candidates as genes that contribute to the stem/precursor cell phenotype and, thus, deserve further attention. Finally, there were 53 genes showing upregulation only in the Do-D1, but the increases were tenuous, only one gene showed an increase in excess of 1.25. The top genes of these two gene expression patterns are listed in the lower panels of Figure 6.

Two well establish notions of CE biology are that the basal limbal cells have low proliferation rates under steady state conditions and that the TA cells concentrate in the limbus-adjacent periphery. If so, one will expect that genes that undergo rapid increases of expression during the cell cycle. The expression pattern of one such gene, BARD1, conformed with this expectation. The BHp values for the D0-D1 and D0-D123 comparisons were 0.60 and 0.49., respectively, but the D123-D4i.1 comparison yielded an inverse ratio (1/R = 1.26) with a BHp value of 6.0 10^-9. Unfortunately, several other genes exhibiting cell cycle expression dependency were not present within the restricted gene list.

Finally, we used the single large 4700-cell central cornea sample to identify genes associated with the expression of high levels of Krt12 in the basal corneal cells. These cells are the most likely to undergo stratification. DEGs were evinced from a comparison of the lowest 1000 Krt12 (Krt12l^o) expression level cells vs the highest

1000 (Krt12^hi). This comparison identified 1563 downregulated and 3391 upregulated DEGs, i.e., 86 % of the genes in the sample. Focusing on the ratios for down to up regulated genes (Table 2) a dynamic pattern is noticeable.

The changes within the limbus and in the transition from limbus to Pe are dominated by gene down regulation. At the early stages of Pe Krt12 change within the Li- adjacent Pe down and up events are even. But at the center of the cornea, when considering the most extreme difference gene down regulation becomes the minor event. The results presented below in Figure 4 and Figure 5, suggests that the reduced down regulation reflects a completion of the down regulation for a large number of genes by the time cells reach the central epithelial zone. Table 7 lists the top down and up regulated genes. The Krt12^hi to Krt12^lo ratio for Krt12 was4.22. Thus all 160 genes listed in the Table 7, undergo fold-changes that exceed the change in Krt12.

3.4. Gene Ontology Analysis

To explore the functional and phenotypic source of the differential gene expression between the different we resorted to a differential gene ontology analysis. The down and up regulated whole genes lists for the Do-D1 and the D123-D4i1 sets were submitted to for Panther software classification at https://www.geneontology.org/. The resulting statistically significant domains overrepresented GO terms unique to either each of the domains were identified. The Table 8 upper panels list the top unique GO terms (UGOTs) of D0 in relation to those terms present in D1, and the respective D1 UGOTs, limited to those with a false discovery rate < 0.01 and ranked according to the extent of gene overrepresentation (FE). The Table 8 bottom panels list top unique D123 and D4i.1, in relation to each other. Due to the hierarchical structure of the gene ontology classification, lists contain multiple terms related to a single metabolic or bioenergetic activity. Thus, to allow inclusion of terms associated with different functions within the list length limitation, we have selectively removed the general category terms.

Four intriguing terms within D123, the domain representing all Krt12 positive cells within the limbus are positive regulation of protein localization to Cajal body, maturation of LSU-rRNA, negative regulation of stem cell differentiation, and regulation of stem cell population maintenance. For the Q4i.1 set it is possible to note lumen acidification of several organelles, desmosome organization, cellular response to arsenious substances and glycolytic processes.

3.5. Gene Expression Correlation Analysis

An alternative inquiry on gene expression within the cornea can be based on the degree of correlation between the change of any gene in comparison with the change in Krt12. Due to the very shallow Krt12 rate of change this approach is not effective for the small Krt12 increase within the limbal D0-D3 domains. The larger cell list for the central cornea, in contrast, provided a basis to obtain robust data. The first 4,600 cells of the Co sample were organized in ascending Krt12 expression levels and split into 20 230- cell quantiles (Q1-Q20) and the correlation of each gene expression change with the changes in Krt12 were calculated using Excel’s Correl function. The analysis evinced 507 genes with a correlation coefficient (CC) > 0.95 and 1268 with a CC > 0.9. Negative correlating genes amounted to 199 and 493 at the < -0.9 and <-0.9 CC levels, respectively. The top genes positively and negatively correlating with the Krt12 expression changes, respectively, are listed in Table 9.

Cytokeratin expression. Each stratified epithelium expresses a unique set of tonofilament-forming cytokeratin pair. This pattern is likely to represent an evolutionary adaptation to optimize each lineage to function in its own specific environment. In the CE, the tissue specific cytokeratins are Krt12 and Krt3. However, the CE expresses a large complement of other cytokeratins, in particular the universal stratified epithelia Krt5 and Krt14. If cytokeratin expression is related to optimization of function and the expression of Krt12 changes with the degree of differentiation, the question arises as to whether the expression of these other tonofilament-forming proteins change and if so, how the changes relate to the change in Krt12. The gathered data presents a unique opportunity to examine this issue. The 6647-gene LiPe list contains 15 cytokeratins within the D4 domain. The 2000-cell D4i subdomain was divided into 20 100- cell quantiles and the average keratin expressions were calculated and plotted from Q1 to Q20. As depicted in Figure 4. four of the keratins displayed meaningful changes with respect to the changes in Krt12 expression. Krt3 expression changes, as expected, tightly correlated with Krt12. The universal cytokeratins Krt14 and Krt5 displayed opposite changes, Krt 5 changes closely followed the Krt12 changes; Krt14 showed a gradual continuous decrease. Finally, Krt75, a low expression cytokeratin, showed a particularly intriguing distribution. Within the D0-D3 zone peaked at D3 region associated with the limbal zone (Figure 4 inset) but decrease sharply in D4i.

4. Discussion

Single cell RNA measurements are usually processed by clustering algorithms that calculate the similarity/ dissimilarity of data points. One of the main uses of this methodology is the identification and characterization of heretofore undetected phenotypically distinct minute cell populations residing within an organ or tissue, e.g., immunosurveillance or incipient transformed cells. In the ocular surface, several studies have employed this methodology to assign cells as primary belonging to one of various distinguishable domains of conjunctival and corneal lineages. Following this clustering, using differential analysis it becomes possible to associate each expressed gene with a specific cluster and calculate a statistical probability of the assignment. Given its significance in the management of limbal stem cell deficiency, the cluster of basal limbal cells, the site of the lineage stem/precursor cells is of particular interest.

The present study uses an alternative approach for the identification of corneal epithelial domains. It is based on the fact that the CE cells undergo a single linear differentiation path, coupled to the hypothesis that the degree of corneal epithelial cell differentiation within the basal cell compartment is reflected in the level of Krt12 expression. The results of the differential expression and correlation analyses performed borne out the validity of the approach. A graphic analysis of the rate of change of Krt12 in the LIPe sample led to the identification of four distinct Krt12 expression domains within the limbal zone. Approximately half of the cells were Krt12-negative and the rest showed Krt12 increasing in three discreet steps. The highest Krt12 level in the limbal domains, though, amount to well less than 1 % of the maximal expression level in the central corneal epithelium. Hence, this very low expression is unlikely to translate into detectable protein levels. Nevertheless, comparative analysis of gene expression in the D0 vs. D1 domains revealed that at the start of intra-limbal Krt12 expression very large number of genes undergo down or up regulation resulting for the first time in the segregation of sub-domains of differentiation within the basal limbus. The bona fide stem/precursor cells can be reasonably expected to reside within D0. This proposition seems supported by the D0 UGOTs (Table 8). Firstly, regulation of hemopoiesis and, embryonic organ development are direct indicators of a relationship to stemness. Secondly, the highest overrepresented term, positive regulation of core promoter binding combined with an overrepresentation of genes involved in the negative regulation of DNA-templated transcription and the negative regulation of RNA biosynthetic process and transcription by RNA polymerase II, yields a picture of cells with high potential for wide gene expression capacity which is prevent from strong expression by mechanisms aimed at slowing RNA and protein synthesis, these are features expected from the quiescent stem/precursor cell. Thirdly, an increased ability to deal with misfolded protein is an expected critical capacity of a cells surviving in the tissue for an extended period. The NADH to ubiquinone aerobic electron transport chain UGOT may reflect a higher preference for dependence of oxygen of cells closely associated with the blood circulations than in the D1 cells, initiating the differentiation towards the anaerobic-preferring feature of the avascular central cornea. The main D1 UGOTs, regulation of plasma membrane organization, intermediate filament cytoskeleton organization, cell-cell junction organization, keratinocyte differentiation, are indicative of the structural changes the differentiated phenotype, including a substantial increase in cell size.

The smaller D2 and D3 domains showed very little difference with the D1 domain, they seem to belong to cells that are mostly unchanged in gene structure from D1 except for a subtle upsurge in genes that undergo much stronger increases as cell transition from the limbal to the peripheral D4 domain. This feature allowed to define the effect of the cell migration from the vascular limbus to the avascular cornea using the whole D123 domains. Consistent with the multiple visible sharp phenotypic changes in the CE cell at the Li-Pe interface, the D123-D4i.1 comparison yielded a very large list of DEGs including about 40 % of the probed genes (Table 2). Probably more significant, the ratio of down to up regulated DEGs underwent a dramatic reversion. Whereas as in the D0-D1 or D0-D123 comparisons downregulated genes were twice as many as upregulated ones, with the transition to the corneal Pe domain, the dominant change was gene up regulation by a 2 to1 margin, i.e., a fundamental change not only in the rate of phenotypic change towards full maturation but in the nature of the cellular maturation processes.

In regard to the D123-D4i.1 comparison itself, a full interpretation of the UGOTs listed in Table 8B is not easily accomplished. From the top UGOTs for the combine D123 domains, it is clear, though, that abundant RNA processing in the nuclear Cajal body that occurs in the Li, ceases or drastically decreases once cells migrate to the Pe. The Cajal body function may be associated the UGOT, including regulation of telomere maintenance via telomerases, since the enzyme mRNA has been found to associate with the Cajal body and its telomere length regulation [19,20]. Other intriguing UGOTs, potentially reflecting the much less differentiated state of the D123 domain in comparison to D4i.1 are, regulation of hematopoietic progenitor cell differentiation, the regulation of stem cell differentiation and regulation of stem cell population maintenance. All these terms are present both the D0, D1 thus, do not show up in either the D0 or D1 Table 8A UGOTs. The D4i.1 UGOT list includes processes that increase organelle acidification and cellular response to arsenic-containing substances. Arsenic inhibits various mitochondrial enzymes leading to the uncoupling of oxidative phosphorylation. Thus, both the acidification and the arsenic response terms may simply reflect the genetic changes associated with the transition of the cornea cell from an aerobic to an anaerobic-able gene configuration.

After completing the determination of DEGs associated with each transition it was possible to categorize genes as continuously changing along the differentiation path or identify genes with selective differential expression at certain stages. This examination identified as many genes undergoing down regulation as up regulation. The latter cohort contains very well recognized genes associated with the corneal phenotype, such as NQO1 and the aldehyde dehydrogenases, both important detoxification genes. How the strong downregulation of many genes contributes to the limbal-corneal phenotype remains to be examined.

Probably the most intriguing cohort, though, is that of the genes that only undergo downregulation with the transition from Do-D1; they are likely to be critical genes for the stem/precursor cell or its survival. The top two genes in the list are NR2F2 and ID3. The first is a retinoid-responsive nuclear factor with a wide gene promotion pattern. At least for a specific transformed cell NR2F2 had been shown to act as a promoter of stemness [17,18]. ID3 is a repressor of basic helix-loop-helix transcription factors and had been shown to support human embryonic stem cell maintenance [20]. Both are particularly interesting subjects for further research.

The additional two comparisons, identify either the genes associated with CE cell maturation in relatively early stages in the peripheral zone next to the limbus (D4i-1-D4i.2 comparison; Table 6), or between the lower and higher extremes of Krt12 expression in the central cornea (Q1-Q20 comparison; Table 7), respectively. The D4i-1-D4i.2 yielded very similar number of down and up regulated genes. A clustering based scRNA-Seq study has led to the suggestion that the top down regulated gene in the D4i.1 to D4i.2 comparison, GPHA2, is a marker for limbal stem cells [14]. In this study, though, GPHA2 was clearly substantially expressed in the lower subdomain of the Pe. Furthermore, a D0-D4 plot for showed a GE distribution similar to that shown in Figure 4 for Krt7. Since Krt75 belongs in the same Table 8 list as GPHA2, it was intriguing to examine the distribution of the latter in the 20 quantiles of D4 and the relationship to other down regulating genes of the D4i.1-D4i.2 comparison. Figure 5 describes the changes in GPHA2 and for the 4 genes with the highest correlation coefficient ( > 0.95) to its distribution. The inset shows the combination of the previously defined D0-D3 domains with the first 5 quantiles of D4. It is clear that GPHA2 achieves it maximum expression at the very start of the cell transition to the Pe domain (D4Q1). The other genes in Figure 5 display a very similar pattern (not shown). The role of this expression pattern in the CE differentiation remains to be discovered.

In summary the present study demonstrates that the application of Krt12 expression as a gauge of the extent of differentiation is an efficient approach for the identification of the dynamics of gene expression changes underpinning the stem/precursor cell phenotype and the progress of CE differentiation. The Discussion provides a very limited example of the analytical possibilities afforded by the results. Future focus on individual genes may help establish a full representation of the coordination of growth and differentiation in the limbal-corneal epithelial lineage.

5. Conclusions

The degree of expression of Krt12, the corneal specific cytokeratin, in corneal basal cells subjected to single cell RNA sequence measurement identify 5 different domains characterized by the rate of Krt12 expression increase. Differential gene analysis between these domains demonstrate that they represent defined stages od differentiation. Four of these stages occur within the limbal zone, the seat of the limbal-corneal stem/precursor cells. These results combined with a study of the Krt12-gene correlation generates a whole picture of gene dynamics during all stages of differentiation.