Spatial Analysis of Single Fiber Cells of the Developing Ocular Lens Reveals Regulated Heterogeneity of Gene Expression.

The developing eye lens presents an exceptional paradigm for spatial transcriptomics. It is composed of highly organized long, slender transparent fiber cells, which differentiate from the edges of the anterior epithelium of the lens (equator), attended by high expression of crystallins, which generates transparency. Every fiber cell, therefore, is an optical unit whose refractive properties derive from its gene activity. Here, we probe this tangible relationship between the gene activity and the phenotype by studying the expression of all known 17 crystallins and 77 other non-crystallin genes in single fiber cells isolated from three states/regions of differentiation, allowing us to follow molecular progression at the single-cell level. The data demonstrate highly variable gene activity in cortical fibers, interposed between the nascent and the terminally differentiated fiber cell transcription. These data suggest that the so-called stochastic, highly heterogeneous gene activity is a regulated intermediate in the realization of a functional phenotype.

Introduction

There is extraordinary heterogeneity in gene expression in single cells; whether it is noise or consequential activity remains unclear (Symmons and Raj, 2016). It is essential to understand the biological consequences of transcriptional heterogeneity at the single-cell level (Snijder and Pelkmans, 2011), but the very isolation of a cell without its neighbors (Klein et al., 2015, Macosko et al., 2015) brings into question the status of the cell as a relevant unit that contributes to the phenotype of a tissue. It is, however, important to recognize that single-cell transcriptomics has allowed the identification of previously unknown functional states and phenotypes (Glotzbach et al., 2011) (Klein et al., 2015, Macosko et al., 2015, Wucherpfennig and Cartwright, 2016, Yamanaka et al., 2013). At the current state of our knowledge the utility of this molecular cataloging of individual cells or cell types, nonetheless, remains unrealized because of the absence of a concurrent morphological, developmental, and/or spatial context.

A rational link between transcriptional heterogeneity and tissue phenotype therefore remains tentative (Lein et al., 2017). A serious shortcoming has been the general lack of correspondence between the single-cell phenotype and the functional organization of the tissue. Keeping in view all the technical limitations of preserving the status of a single cell as a component of a multicellular, functional entity (tissue), we present the developing eye lens as a paradigm for spatial transcriptomics.

Contrary to the perception of an inanimate bag filled with protein, the ocular lens is a highly organized cellular structure, which is remarkably simple in terms of its cellular composition; it contains only two cell types, the progenitor anterior epithelium and the fiber cells derived from it (Figure 2A) (Bassnett et al., 2011, Cvekl and Zhang, 2017). The epithelium differentiates at its edges (lens equator, Figure 2A), adding new fiber cells on an already existing fiber mass; its youngest fiber cells are thus at the periphery, and the oldest (terminally differentiated), at the center or in the nucleus, which makes the visual axis of the lens. In the adult lens, 95% of the volume of this tissue is occupied by fiber cells (Bassnett et al., 2011, Costello et al., 2016), which synthesize crystallins to achieve transparency and maintain a specific refractive index (RI) (Iribarren, 2015).

Figure 2

Unsupervised Clustering Reveals Pervasive Heterogeneity in the Ocular Lens Fiber Cells

(A) A cross section of the postnatal day two (PND02) mouse lens. We developed a procedure based on temporal release of fiber cells from the lens that allowed us to collect individual fibers (Figure 1) and assign them to one of the specific states/stages of the developing lens (indicated by colored bars; each bar, 200 μm; red, equatorial; green, cortical; blue, nuclear (see Transparent Methods in the Supplemental Information).

(B) Heatmap generated from the gene expression data (446 single fiber cells and 94 genes including 17 crystallins and 77 non-crystallins). Notice two clusters of high expression, open arrow #1 (13 crystallins) and open arrow #2 (58 fiber cells). The arrow #2 clusters fiber cells, which show very low expression of crystallins but high expression of non-crystallins; these fibers are predominantly cortical (green, follow sample groups on the x axis). Note that Hsf4, the major heat shock transcription factor of the postnatal lens is predominantly expressed in these fibers (indicated by a dotted arrow). Note also that in the left side of the figure (open arrow # 3, n = 93), the fibers show very poor expression of non-crystallins. Single fiber cells are annotated at the bottom of the heatmap in alphanumeric values (S1–S460); 388 single fiber cells (87%) express all 17 crystallin genes; ∼10 genes (including Mmp2, Rab5a, Lamp2, Atp1a3, Myl1, Atp1a2, Mcpt4, Fgf7, Slamf8, and Opn1sw) show significantly low expression in all fiber cells. Long broken arrows with gene names (Hsf4, Hspb1, Hspb2, Hspb6, Itgb1, Igf2, Id2, Bfsp2, Rab27b, and Cd84) indicate examples of heterogeneity. Predominant expression of a gene in a specific group of fiber cells gives those fiber cells their molecular identities. Other genes of interest are indicated in bigger fonts on the right. The sample group colors, red (equatorial), green (cortical), and blue (nuclear), provide spatial information. See also Figures S1 and S2.

A critical feature of this phenotype (transparency), which is set up very early in development (Peetermans et al., 1987) is that it is characterized by a gradient of RI, which follows the gradient of protein concentration, low in the periphery and highest in the center/visual axis (Campbell and Hughes, 1981, Chakraborty et al., 2014, Philipson, 1969, Pierscionek et al., 1987) (Pierscionek and Regini, 2012). A smoother gradient of RI assures pinpoint focus without spherical aberration. Importantly, therefore, every fiber cell becomes an optical unit that must reconcile to its spatial status and molecular uniqueness, attained through specific crystallin expression (Mahendiran et al., 2014) to refract light into the eye to make vision possible.

Given this developmental history, the lens lends itself to the isolation and study of single fiber cells, region by region in a spatially contextual fashion (Figures 1 and 2A). For example, the gene activities in single fiber cells can be followed spatially from the equator to the cortex to the nucleus (center) of the lens; the center of the lens has the highest protein concentration and the highest RI as its phenotype. At postnatal day 2 (PND02), differentiation from the overlaying epithelium has already generated an equatorial region (nascent fibers), a cortical region (differentiating fibers), and a nuclear region (terminally differentiated fibers), which makes the future visual axis. Importantly, the lens at this stage is pliable and therefore allows us to peel off fiber cells, one at a time, by gentle shaking of the fiber mass. The younger superficial fiber cells come off quicker, whereas the deeper fibers take longer. Based on the time it takes for fiber cells to come off, we have cataloged them as equatorial fibers, cortical fibers, and nuclear/central fibers (see Transparent Methods section).

Figure 1

Isolation of Single Fiber Cells from Two-Day Old Mouse Lens

(A) A picture of the whole postnatal day two (PND02) mouse lens.

(B and C) (B) Lens (fiber mass), without the epithelium and the capsule. (C) The schematic (blue) depicts three distinct regions (outer equatorial, middle cortex, and inner lens nucleus). The fiber mass was processed to collect individual fiber cells by incubating it with shaking (Transparent Methods). The outer equatorial fiber cells (length: 100–120 μm) were collected in first 10 min; the cortical fiber cells (length: 220–240 μm) were collected in next 10 min, and finally the nuclear fiber cells (length: 600–640 μm) were collected in the last 10 min. Scale bar, 500 μm. In this study, we did not see a direct correlation between fiber cell length and gene activity with respect to 94 genes (including crystallins).

In this investigation, we have studied single fiber cells of the PND02 mouse lens (Figure 1) using microfluidic quantitative qRT-PCR in a Biomark microfluidics system (Fluidigm, Inc.). The analytical capability and low technical noise of this system allows decisive assessment of transcript levels. We have focused on known crystallin genes (all 17 of them) and other seventy seven genes of known relevance to the lens phenotype (Table S1).

Results

There Is Unexpected Heterogeneity of Gene Expression in Individual Fiber Cells in the PND02 Lens

Single fiber cells (n = 446), isolated from PND02 mouse lenses were interrogated with probes for 94 genes. Unsupervised clustering of the data reveals remarkable heterogeneity (Figure 2B). As expected, the prevalent expression of crystallins is unmistakable (open arrow #1). A cluster of 13 crystallins is seen (Figure 2B, dotted rectangle); two other crystallins seem to associate with other genes, for example, Cryba2 with Uba52 and Crygn with Cd24a (Figure 2B). Similar associations involving same genes are also seen in region-specific clustering analyses (Figures S2A–S2C); they are not altered by subtle lens-to-lens (animal to animal) variations (Figure S1). Euclidean distances can change based on scaling and average expression levels, yet we see consistency in these gene associations, in clustering done with 446 fiber cells (Figure 2B) and in clustering done with about 150 each, region-specific fiber cells (Figures S2A–S2C). Another important observation about the data in Figure 2B is that the cluster of 13 crystallins is retained in the equatorial (Figure S2A) as well as in the cortical fiber cells (Figure S2B), but it is reduced to a cluster of nine in the nuclear fibers (Figure S2C; see also Figures 7A–7C).

Figure 7

Crystallin Expression in Single Fiber Cells Reflects the Phenotype of the Tissue

(A–C) Correlation matrix plots (Pearson R) for single fiber cells are presented region by region, starting from (A) equatorial to (B) cortical and finally to (C) nuclear, following the temporal and molecular progression of differentiation. The expression of nine crystallin genes (indicated on the top left in each plot) are positively correlated with each other and negatively correlated with the majority of non-crystallin gene activity in the equatorial and cortical regions (A, B, dark blue). Scale bars, +1 to −1 on the right. Similar data were obtained when all regions are analyzed together (Pearson correlations are not affected by scaling).

(D) Represents average expression (Log2) values in three regions as line plots. Nine crystallin genes, same genes as seen in the correlation plots (Crygf, Cryaa, Crygb, Cryga, Cryge, Cryba1, Crybb1, Crybb3, and Crygc) show increased expression in the nuclear/central fibers, which are part of the future visual axis, whereas other eight crystallins (Crygd, Cryab, Cryba4, Crybb2, Crygn, Cryba2, Crygs, and Crybg3) show a decrease.

(E) A schematic representation of gene activity in three spatial states of differentiation in the developing mouse lens. Note that the molecular progression starts in the equatorial fibers (red) and progresses through the cortical fiber cells (green) to finally express as terminal differentiation in the nuclear (blue) region. See also Figure S9.

Among the non-crystallins, remarkably, Uba52, a hybrid gene with an ubiquitin 5′ domain and a ribosomal Rpl40 domain on the C terminus (Kobayashi et al., 2016), is one of the highly expressed genes (Figure 2B); comparatively, ubiquitin B is very poorly expressed. Fgf2 is expressed in almost all fiber cells. This is significant because Fgf2 is known to be intimately associated with lens morphogenesis (McAvoy and Chamberlain, 1989). Fgf7 (Weng et al., 1997), on the other hand, is generally poorly expressed (here it is expressed in about 30% of the fiber cells). Interestingly, Rab27b (an important member of Rab GTPase family, involved in exosome secretion [Ostrowski et al., 2010]) and Cd84, a surface glycoprotein (a member of the signaling lymphocyte activation molecules [SLAMs] and used as a differentiation marker on hematopoietic progenitor cells (Zaiss et al., 2003), show elevated expression only in specific group of fiber cells (Figure 2B, dotted arrow). A number of other interesting non-crystallin gene activities are indicated in Figure 2 in a bigger font (on the right side of the figure).

What, however, commands attention is the absence of high crystallin expression within a group of fiber cells, which predominantly come from the cortical (green) region (Figure 2B, open arrow #2). These cells make 13% (n = 58) of all fiber cells analyzed. Interestingly, this group of cells expresses a large number of non-crystallin genes, at a much higher level, as if to compensate for the lack of robust crystallin expression (Figure 2B, open arrow #2). Regional clustering of the equatorial, cortical, and nuclear fiber cells further confirms that these fiber cells are predominantly from the cortical region (Figures S2A–S2C, black arrowheads with line brackets, at the bottom of the heatmaps). It is important to note that this group of fiber cells is made of cortical fibers from all five lenses used in this study (see Figure S1). There are many gene activities, which are present only in a few fiber cells, as indicated by dotted arrows with gene names on them (Figure 2B), suggesting the presence of many “cell types,” but lack of crystallin expression in lens fiber cells is unknown. Fiber cells in this group express some genes highly and others very poorly, generating fiber cells with individual identities. For example, we may simply assign a molecular tag to a fiber cell based on the expression of Hspb1, or that of Hspb2 or Hspb6 (Figure 2B). It is interesting to note that both calpain 3 (Capn3) (De Maria et al., 2009) and DNase 2b (Nishimoto et al., 2003), genes known to be expressed in deeper layers of fiber cells are highly expressed in these cells. Interestingly, these fiber cells do not express crystallins including Crygf and Cryge (De Maria and Bassnett, 2015).

Cortical Fiber Cells Show Highest Variability of Gene Expression

To ascertain the variability of gene expression in the three contiguous regions of differentiation (nascent fiber cells, cortical fiber cells, and terminally differentiated nuclear fiber cells), we calculated the variance of crystallin and non-crystallin genes (Figures 3A and 3B). The gene activity for both the crystallins and non-crystallins shows highest variability in fiber cells of the cortical region (Figures 3A and 3B, green bars). Obviously, the group of 58 fiber cells (referred to above) contributes to the variability substantially.

Figure 3

Gene Expression in the Developing Cortical Fibers Is Highly Variable

(A) Variance in crystallin gene expression.

(B) Variance in non-crystallin gene expression. Cortical fiber cells (green bars) show highest variance in >95% of genes analyzed compared with equatorial (red) and nuclear fiber cells (blue). The variance values are plotted based on highest to lowest in the nuclear fiber cells.

(C) Violin plots of gene expression in individual 446 fiber cells from equatorial (red), cortical (green), and nuclear (blue) regions. The distributions in equatorial and nuclear regions are similar and mostly unimodal (y axis, expression (log2); x axis, number of fiber cells). The cortical fiber cells (green) stand out (one of these columns has been indicated with * on top). These cells show bimodal/multimodal distributions indicating expression in more than one population of cells. For region-specific violin plots, see also Figures S3–S5.

Interestingly, this variability contributes to variability in the fiber cell populations as demonstrated by the violin plots (probabilistic population distributions of gene activity) shown in Figure 3C. The violin plots for all the 446 fiber cells reveal bimodal or multimodal distributions in cortical fibers (Figure 3C, green, one of the columns of these distributions is indicated with an asterisk on the top) suggesting high variability and differential expression of the same gene in different sub-populations (Shalek et al., 2013). There is lower heterogeneity in the equatorial (red) and nuclear (blue) fiber cells, which are mostly unimodal (for violin and boxplots of fiber cell gene activity in individual regions see Figures S3–S5).

Principal component analysis (PCA) reveals that the three-dimensional space occupied by the expression levels of 94 genes sifts the 446 fiber cells into roughly five clusters (Figure 4A). Cluster #1 contains 58 fibers, predominantly cortical; it contributes >80% variance. This cluster is characterized by very low crystallin expression (see Figure 4D heatmap and the corresponding partial violin plot in Figure 4E). Partial violin plots for clusters #1–3 are presented in Figures 4E, 4G, and 4I, respectively. In cluster #1, there is very little Cryaa (the most highly expressed crystallin in the ocular lens) and Cryba1 (Figure 4E, asterisks), yet we see Cryab, Cryba2, Cryba4, Crygd, Crygn, and Crygs, suggesting possibly the non-refractive activities of these crystallins (Bhat, 2004).

Figure 4

PCA Reveals Multiple Populations of Fiber Cells

(A) PCA compartmentalizes the expression of 94 genes in 446 fiber cells into roughly five clusters.

(B) A PCA scree plot. Cluster #1 (PC1) contributes about 86% variance; cluster # 2 (PC2), 20%; and cluster # 3 (PC3), about 8%.

(C) t-SNE (distributed stochastic neighbor embedding analysis [van der Maaten and Hinton, 2008]) scatterplot reiterates the separation of the cell populations seen in the hierarchical clustering (Figure 2B).

(D and E) (D) Heatmap and (E) violin plot (partial) for cluster #1 (n = 58, 42 cortical, 10 equatorial, and 6 nuclear). Out of 11 crystallins (D, lower right), only Crygd, Crybb2, and Cryba4 are expressed in 44 of 58 fibers. A few of these fibers show expression of Crygf in cells, which express Hspb2 and Hspb6 (broken arrows), suggesting expression in discrete cell populations. Note also specific patterns of Hspb1, Igf2, Id2, and Itgb1 expression (top right broken arrow). These fibers do not express Cryaa and Cryba1 (* on violin plot). Note bimodal distribution in cortical (green) fibers (arrowheads).

(F and G) (F) Heatmap and (G) corresponding violin plots for cluster # 2 (n = 28). The majority of cells are from the nuclear (blue) region; all 17 crystallins are expressed, mostly at high levels.

(H and I) (H) Heatmap and (I) the corresponding violin plots for cluster # 3 (n = 48, 21 equatorial, 21 nuclear, 6 cortical). Note high expression of most crystallins. Examples of specific cellular distributions in these violin plots are indicated (▿) (complete violin and boxplots are presented in Figures S6–S8). Note that Hsf4 expression is lowest in cluster #3 and highest in cluster #1(•). Clusters #4 (n = 161) and #5 (n = 140) are both constituted with equal representations from equatorial, cortical, and nuclear regions. Both these clusters show higher expression of Cd84, Rab27b, Sparc, and Uba52 not shown). See also Figures S6–S8.

The cluster #2 is predominantly composed of nuclear fibers, blue (n = 28), which show high crystallin expression (Figures 4F and 4G, arrowheads); it contributes 20% variance (complete violin and boxplots for clusters #1–3 are presented in Figures S6–S8; the data for other two clusters [#4 and #5] are not presented). Interestingly cluster #3 has almost equal representation from the equatorial and nuclear fibers in addition to six fiber cells that it contains from the cortical region (green). The fiber cells in this cluster show robust crystallin expression with similar violin plots, but the cluster also has cells with certain specific distributions (see open arrowheads in Figure 4I), which indicates the presence of additional cell-type variability.

The Equatorial and Nuclear Fibers Present Similar Transcriptional Profiles

A comparison of the equatorial and cortical (Figure 5A) and the cortical and nuclear fiber cells (Figure 5C) reveals a significant number of differentially expressed genes. Noticeably, we do not pick any differential activity between the equatorial and nuclear fibers (Figure 5B). It is important to recognize that cortical fiber cell gene activity follows the equatorial gene activity, both physically (side by side) as well as temporally, and nuclear activity follows the cortical activity. The equatorial and the nuclear fibers are separated by the cortical fibers physically, yet the similarity of gene activity stands out. Interestingly, it is the non-crystallin gene expression (including the expression of crystallins, like Crygs, Cryba2, and Crygn, possibly because of their catalytic roles [Bhat, 2004]) that predominates in the cortical region fibers (Figures 5A and 5C).

Figure 5

Differential Expression of Genes in Fiber Cells Isolated from Different Regions

(A) One-way ANOVA analysis of differential gene activity. (y axis = −log10, p values; x axis = log2 fold change). Seventy-five genes are differentially expressed between the equatorial and cortical fibers; crystallins are highly expressed in the equatorial fibers. Note that crystallins, including Crygs, Crygn, and Cryge are part of the non-crystallin gene activities ( not marked).

(B and C) (B) Demonstrates that equatorial and nuclear fiber cells show no significant differences in expression. Cortical versus nuclear fibers reveal that 82 genes are differentially expressed (C). Note Crygs, Cryba2, and Crygn in the cortex, along with other non-crystallin genes.  = Crystallins,  = non-crystallins.

The data presented in Figures 6A–6C compare the average expression of each gene in the fiber cells of the three regions (the comparison is made between two regions at one time). Again, the gene activity in two separate regions, namely, the equatorial and nuclear fibers, is strikingly similar; no gene activity stands out (Figure 6B). Although the cortical gene activity is higher, the crystallins stand out in comparisons between the cortical and the equatorial fibers (Figure 6A, asterisks) as well as between the nuclear and cortical fibers (Figure 6C).

Figure 6

Equatorial and Nuclear Fiber Cells Have Similar Expression Profiles

(A–C) Comparative analyses of the expression of individual genes in fiber cells from different regions. y axis, average gene expression (log2); x axis, genes arranged based on probability (p < 0.05, significant to p > 0.05, not significant). The comparison between the equatorial and the cortical fibers (A) and the cortical and nuclear fibers (C) shows that it is the crystallins that stand out (open asterisks), although the non-crystallin activity is higher in both comparisons. (B) The expression levels in the equatorial and nuclear fibers are very similar; no gene activity stands out, including the crystallins (•).

Crystallin Expression in the Nuclear Fiber Cells Is Deterministic

Hierarchical clustering (unsupervised) of all the fiber cells revealed that out of 17 crystallins, 13 crystallins made a cluster (Figure 2B, open arrow #1). The other four crystallins, namely, Crygn, Crybg3, Crygs, and Cryba2, associate with different genes (Figure 2B). Interestingly, hierarchical clustering of fiber cells in each region (Figures S2A–S2C) indicated similar patterns of clustering, which was largely maintained in the equatorial and cortical heatmaps; however, in the fiber cells from the nuclear region, the cluster of crystallins was reduced to nine (Figure S2C, dotted rectangle). Figures 7A–7C indicate that the expression of these nine crystallin genes, namely, Cryaa, Cryba1, Crybb1, Crybb3, Cryga, Crygb, Crygc, Cryge, and Crygf, is highly correlated (Pearson correlations ≥0.8–1.0). These genes first appear in the equatorial region (Figure 7A); their expression is sustained through the cortical region (Figure 7B), and they finally become part of the nuclear/central region of the lens (Figure 7C). It is noteworthy that at the single-cell level, Cryab expression is not a part of this group of nine crystallins, which further corroborates the status of the two α-crystallins (Cryaa and Cryab) as separate, independent proteins and not as a single monolithic structural protein (Gangalum et al., 2012).

Line plots of the log2 expression levels of crystallin genes and non-crystallins, across all three regions reveal that the expression of nine crystallin genes goes up in the fiber cells of the nuclear region (Figure 7D, dotted rectangle) as if on cue, whereas the expression of the eight other crystallins declines along with the other 77 non-crystallins (Figure S9). It is important to recognize that the gene activity does not follow the size (length) pattern of the fiber cells (Figure 1).

Discussion

In this article, we present a snapshot of the development in the mouse lens at PND02. This snapshot is not a static picture but assembled of gene activities in single fiber cells, which represent three different states/regions of differentiation. A number of significant studies including intestinal crypts (Dalerba et al., 2011), lung epithelium (Treutlein et al., 2014), and tissues like liver (Halpern et al., 2017) and heart (Skelly et al., 2018) have all revealed developmental emergence of heterogeneity of gene expression associated with specific cell types. The data presented in this investigation, however, almost literally puts variability of gene expression, at the single-cell level, stark in the middle of the molecular progression that culminates in the functional phenotype of the tissue (Figure 3). The fiber cell morphogenesis that commences in the equatorial region ends at the center of the lens (which makes the visual axis). The data clearly suggest that the equatorial fiber cells and terminally differentiated fiber cells at the center of the lens have very similar transcriptional profiles, whereas the cortical fiber cells that separate the two show highly variable gene activity (Figures 3, 6, 7, and S9). The molecular commonality in the processes occurring at these two morphologically separate places in the lens warrants caution because we are only looking at about 100 genes; it is possible that we do not have those genes represented here that may vary in the fiber cells of the two regions.

The very nature of the fiber cell morphogenesis in the post-embryonic lens anatomically positions the oldest differentiated fiber cells at the center and the youngest at the periphery of the lens. Thus, at one time point (namely, the PND02 lens), we can weave a pathway of molecular progression that has a temporal as well as a morphological dimension to it. As the lens keeps growing (adding fiber cells to an already existing fiber mass), a nascent fiber cell early in the development may become part of the lens nucleus (the future visual axis) in a matter of days or weeks, although this process slows considerably with age (Augusteyn, 2007). It is obvious therefore that what shapes the lens nucleus/visual axis happens very early in development; importantly, these fiber cells and their constituents must last a lifetime, and remain transparent (Andley, 2007, Bassnett et al., 2011, Bhat, 2003, Costello et al., 2016).

What is quite interesting is that we see almost all known crystallins expressed here, but they are not uniformly expressed at the same level in all fiber cells. It is important to recognize that in the analysis of whole tissues or populations of cells, although there is no direct one-to-one correspondence seen between protein and RNA, it is becoming increasingly clearer that at the single-cell level, there is appreciable degree of correspondence between transcription and translation (Li and Biggin, 2015). If we, therefore, take the liberty of comparing the expression of crystallin genes here, with a 2D gel pattern of a newborn mouse lens (Ueda et al., 2002), an appreciation of working with the single cells versus with a population of cells becomes apparent. For example, the two proteins (Crybb2 and Crygs), which are least represented in the newborn gel pattern of the pooled lens homogenate, are highly expressed in a group of cortical fiber cells (Figure 2B, open arrow #2).

There are earlier studies wherein whole lens epithelia and whole fiber masses have been used for microarray analyses. In some studies, whole eyes were used for studying γ-crystallins (Goring et al., 1992). Studies involving in situ hybridizations with crystallin probes during early mouse lens development did not find many γ-crystallins in the equatorial region on postnatal day 1 (Treton et al., 1991); on the other hand, we find differential expression of γ-crystallins in individual fibers in the PND02 lens. These early pioneering investigations not only highlight the limitations of the technologies used but also point to the incongruences between population (whole tissue) studies with single-cell studies. Thus determinations of which crystallin is expressed when, within a population of cells, will only present an average picture, which may not allow meaningful functional correlations.

Previous work (Lieska et al., 1992) studiously dissected different regions of the adult bovine lens. This work demonstrated that that the in vitro translation of RNA isolated from the lens nucleus generates a profile similar to the protein composition of the adult lens nucleus. This clearly suggests that not only is the nuclear RNA intact and possibly functional, but more importantly, it represents the phenotype of the nuclear fiber cells. By extrapolation, therefore, the RNA transcripts we have detected in the nuclear fiber cells represent the phenotype (proteins) on these cells. Our work now provides a window into the dynamics of the differentiation process, which is harder to probe in an adult lens. The transcriptional profiles discovered in single fiber cells (Figures 5 and 6) reveal that the cellular progression from the nascent fibers to terminally differentiated fiber cells is not attended by a gradual ascension or decline of specific gene activities. Based on only 94 genes studied here, this process is book-ended by similar transcriptional profiles in the equatorial (nascent) and terminally differentiated fiber cells. The transition between the two (between the equatorial and the nuclear) is a phase of highly variable transcriptional activity in cortical fiber cells (Figure 3), suggesting that heterogeneity may be a critical intermediate in this process.

The high variability of gene expression in the cortex is not the result of gene activity variations in the same population of cells. Different populations of cells as indicated by violin plot distributions create this variability (Figures 3 and 4). It is interesting to recognize that the first cluster of cells in the PCA analysis (cluster #1, Figure 4) is predominantly composed of cortical fiber cells, which show very poor crystallin expression, and the second (cluster #2, Figure 4) is predominantly composed of nuclear fiber cells, which show high crystallin expression.

The significance of the group of fiber cells in the cortex, which do not synthesize crystallins appreciably (Figure 2B, arrow # 2, and Figure 4, cluster #1), is unknown. These cells do not express Cryaa, the chief crystallin of the vertebrate ocular lens (Horwitz, 2003). Proteins, other than crystallins, are known to have lower RI increments (Mahendiran et al., 2014, Zhao et al., 2011b); low crystallin expression therefore may suggest the presence of fiber cells with a slightly lower RI, which may have specific locations within the lens (Bassnett et al., 2011).

Non-crystallin gene activities, almost all of them, with very few exceptions show an increase in the cortical fiber cells (see Figure S9). These non-crystallin activities contribute to specific cellular identities (for example, see heat shock proteins Hspb 1, 2 and 6, and Igf2, Id2, Rab27b, and Cd84 in Figure 2B). This suggests a multitude of physiologies attended by molecular transformations, which prepare the fiber cells for passage to terminal differentiation. These physiologies, at one end, encompass the arrival of differentiating fibers from the equatorial region, and at the other, progression of these cells into terminal differentiation in the nucleus/center of the lens (Bassnett et al., 2011, Costello et al., 2016, Rowan et al., 2017, Subczynski et al., 2017). A critical role of cell types in region-specific appearance of pathologies (cataracts) is thus a possibility. For example, in MIP −/− (major intrinsic protein-null) mice, one of the important gene products that is downregulated is Hspb1 (Bennett et al., 2016).

We see that Lim2 (lens intrinsic membrane protein) is a highly expressed gene in the cortical fiber cells, followed by Cd24a, Bfsp1 (filensin), vimentin, Uba52 (a hybrid ubiquitin-ribosomal protein Rpl40 gene), and Capn3 (calpain 3). These gene products point to a large number of physiological activities, which encompass rearrangement of the membrane, cytoskeleton, proteolysis, and intercellular signaling. More work is required to separate and identify these activities for a specific understanding of each of these two progressions. The high non-crystallin gene expression in the cortex supports the argument that heterogeneity of gene activity may suggest multiple routes to the final functional phenotype, which in the present case is the synthesis of specific crystallins that make the nuclear region of the lens transparent and of high RI.

An important observation is that the correlated crystallin expression starts in the equatorial cells, persists through the highly robust gene activity in the cortical fiber cells (Figure 7B), and finally emerges at the center of the lens in the future visual axis where all other gene activities decline (Figures 7D and S9). We interpret these data to indicate that the highly correlated expression of nine crystallins (out of 17) suggests their refractive role in the future visual axis, whereas the other eight crystallins may have non-crystallin (non-refractive) functions (Bhat, 2004). The expression of crystallins also indicates that gene products of terminal differentiation need not appear at the end but may start at the beginning of this process and may, in fact, have a directive influence on the progression.

Notably, out of the nine crystallins (Figures S2C, 7C, and 7D) expressed in the nuclear fiber cells, five are γ-crystallins, which are known to have high RI increments (0.203 mL/g) (Zhao et al., 2011a) (Mahendiran et al., 2014, Pierscionek et al., 1987). These data provide a tangible and phenotypically rational purpose for the higher expression of these genes in the fiber cells of the lens nucleus/visual axis. The data suggest a deterministic program that leads to the expression of specific crystallin genes in the nuclear region of the developing lens. This relationship between the specific gene activity and the relevant phenotype will need further elucidation.

Although it is rather obvious to focus on the refractive increments of the γ-crystallins, which contribute to transparency, the physiology and structure of the fiber cells as modulated by gene expression in the early lens may have important functions, which contribute to the development of emmetropia (which is critical for the realization of focused functional vision in infancy) (Iribarren, 2015).

Finally, our data point to the existence of a causal continuum that commences with gene activity, which contains differentiation-specific genes, whose expression is sustained from the beginning until the end. The molecular progression detailed here in single fiber cells suggests that the process of differentiation in the ocular lens starts and ends with similar transcriptional profiles; it is highly significant that the two are connected via highly heterogeneous gene activity, which may appear “stochastic” and yet culminate in phenotypically relevant expression. Because we know the spatial location of this heterogeneity, a more incisive examination may explain how highly variable gene activity may facilitate deterministic emergence of a phenotype in terminal differentiation (Figure 7E).

Limitations of the Study

Procedures used in this investigation, allow us to place an isolated fiber cell within a broad but contiguous region(s) of morphological and molecular activity within the developing lens; the resolution is thus coarse. An analysis of fiber cells with spatial information that places them next to each other would ultimately be required to reveal the finer nuances of the gene activity and the phenotypic organization that makes this tissue functional (transparent). The data presented in this manuscript, however, clearly suggest that the ocular lens makes a meaningful paradigm for the study of the relevance of single cells to an understanding of the emergence of the tissue phenotype; nonetheless, there are vital limitations in these investigations. For example, first, the procedure we have employed to isolate cells is manual and time consuming; there is a need for a method that isolates more cells quicker, in less time, and in a way that maintains the physiological and morphological integrity of individual fibers while keeping their spatial information intact. Second, the use of proteases and/or storage in less-than-optimal buffers (physiological or otherwise) can introduce both visible (morphological) and invisible (molecular) changes, which may not reflect the native state of the cell within the tissue. Third, technical bottlenecks including the size and shape of cells within a tissue at different ages, their RNA content, and library preparations across hundreds and thousands of single cells needs to be standardized for each tissue (Oldham and Kreitzer, 2018) and here in the ocular lens for each developmental stage/age. Fourth, single-cell transcriptomics has a critical limitation in detecting molecular indicators of cellular interactions, which may only materialize when the cells are in contact with each other, in a fashion that is constrained by the tissue morphology and physiology. This information may only come from the deep sequencing of the whole-tissue RNA.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.