Transcriptome Meta-Analysis Uncovers Cell-Specific Regulatory Relationships in Embryonic, Juvenile, Adult, and Aged Mouse Lens Epithelium and Fibers

Abstract

Purpose: The lens transcriptome has been examined using microarrays and RNA-sequencing (RNA-seq). These omics data are the basis of the bioinformatics web-resource iSyTE that has identified new genes involved in lens development and cataract. The lens predominantly contains epithelial and fiber cells, and yet, presently, iSyTE is based on whole lens data. To gain cell-specific regulatory insights, we meta-analyzed isolated epithelium and fiber transcriptomes from embryonic/postnatal, adult and aged lenses.

Methods: Mouse lens epithelium and fiber transcriptome public datasets at embryonic (E) and postnatal (P) stages E12.5, E14.5, E16.5, E18.5, P0.5, P0, P5, P13, and age one month, three months, six months, and two years were analyzed. Microarray or RNA-seq data were analyzed by appropriate methods and compared to other resources (e.g., Cat-Map, CompBio).

Results: Across all RNA-seq datasets examined, 2466 genes are differentially expressed between epithelium and fibers, of which 106 are cataract-linked. Gene ontology enrichment validates epithelial and fiber expression, corroborating the meta-analysis. Whole embryonic-body–in silico subtraction and other analyses identify several new high-priority epithelial- and/or fiber-enriched genes (e.g., Casz1, Ell2). Furthermore, new insights into cell-specific regulatory processes at distinct stages are identified (e.g., ribonucleoprotein regulation in E12.5 epithelium). Finally, this data is made accessible at iSyTE (https://research.bioinformatics.udel.edu/iSyTE/).

Conclusions: This spatiotemporal transcriptome meta-analysis comprehensively informs on epithelium- and fiber-specific regulatory processes in developing, adult and aged lenses. Notably, it includes the first description of an embryonic stage (i.e., E12.5) representing early primary fiber differentiation, thus informing on the initial transcriptome changes as lens cell-types are readily distinguishable.

Over the past two decades, many studies have focused on the analysis of the lens transcriptome by microarrays^1–18 or high-throughput RNA-sequencing (RNA-seq).^{1,12,13,16,17} This rich genome-wide transcript expression data on various stages of the lens has been used in developing the bioinformatics web resource iSyTE (integrated Systems Tool for Eye gene discovery),^11,19,20 which has facilitated identification and characterization of new genes and pathways linked to lens development and cataract pathology.^{10,14,16,21–36} The lens is made of two principal cell types, namely, epithelial cells in the anterior region, and fiber cells that span the entire anteroposterior space to majorly contribute to the mass of the lens.^37,38 In early development, as the lens pit progresses to form the lens vesicle, cells located in the anterior region maintain epithelial fate whereas cells located in the posterior region differentiate into primary fiber cells.^37,38 Cells in the epithelium maintain their potential to proliferate, but only those located in the “germinative zone” proliferate, and cells located near the equatorial region termed the “transition zone,” exit the cell cycle and initiate differentiation into secondary fiber cells. The fiber differentiation program involves downregulation of “epithelial” genes (e.g., Cdh1, Foxe3, etc.) and upregulation of “fiber” genes (e.g., several crystallins, gap junction proteins, Aqp0 (Mip), Caprin2, Tdrd7, etc.). 

Indeed, the epithelial and fiber cell transcriptomes differ significantly in order to serve their specific cellular functions, as the epithelium is the “repertoire” of proliferating cells (and, hence, can be considered analogous to “stem-like cells” in the lens), whereas differentiated fiber cells express high levels of select proteins and undergo dramatic cell shape changes, as well as organelle degradation, all of which are necessary for lens transparency. Thus mutation or mis-regulation of genes expressed in either of these cell types are known to cause lens defects and cataract.^39,40 

Therefore, to gain insights into lens epithelial and fiber biology, it is important to examine cell-type specific transcriptomes. However, the lens transcriptome data that is presently used in iSyTE is generated on whole lenses and therefore is not optimal for understanding the cell-type–specific impact on lens biology or pathology, representing a significant knowledge gap. Incidentally, although human data is limited,⁶ several datasets on isolated mouse lens epithelium and fiber cells have been generated. Thus, to address the current knowledge-gap in iSyTE and to gain insights into stage-specific changes in these cell types, we meta-analyzed publicly available transcriptome datasets on isolated mouse lens epithelium and fiber cells. These data are generated using microarrays or high-throughput RNA-sequencing (RNA-seq) on embryonic, early postnatal, adult, and aged mice. The stages covered are embryonic day (E) 12.5 (unpublished), E14.5,⁴⁰ E16.5,⁴⁰ E18.5,⁴⁰ postnatal (P) 0 (newborn, referred as P0.5 by Cvekl group and P0 by Robinson group [for simplicity, we refer to these data as P0a and P0b, respectively]),^39,40 P5,⁷ P13,⁴¹ and age one month,¹⁸ three months,⁴² six months,⁴³ and 24 months.⁴² Although from diverse sources and based on different platforms (Affymetrix microarray chip, Illumina microarray chip, custom microarray, RNA-sequencing), meta-analysis of these datasets effectively identified epithelial and fiber expression of established lens genes. 

Additionally, we also analyzed other datasets on isolated epithelium, and outer cortical and inner cortical fiber cells generated on Affymetrix microarray,¹⁸ as well as datasets on isolated young and mature fiber generated on a custom-ordered 20K Codelink Bioarray.⁷ This analysis is the first to comprehensively inform on cell-type–specific lens transcriptome changes from early embryonic and postnatal stages through adulthood and aging. Particularly, this is the first report that provides detailed analysis of isolated epithelium and fibers in an early stage of lens development, namely at E12.5, when these cells are readily distinguishable. Finally, to allow user-friendly visualization of these data, genes exhibiting expression or enriched-expression in the isolated lens epithelium or fiber cells are made publicly accessible at iSyTE (https://research.bioinformatics.udel.edu/iSyTE/), to facilitate studies in lens biology and cataract gene discovery. 

Methods

Mouse Studies

All mouse assays in this work conform to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research. All experimental protocols involving animals were reviewed and approved by the Animal Care and Use Committee of the University of Delaware. Wild type ICR mice (Taconic, Albany, NY, USA) were used for in situ hybridization analyses, which was performed as described previously.¹¹ Mice housed in a 14-hour light to 10-hour dark cycle were bred for timed pregnancies, and the detection of vaginal plug in the morning of the day was considered as embryonic day E0.5. 

Meta-Analysis of RNA-Sequencing (RNA-seq) Datasets on Isolated Lens Epithelium and Fibers

Multiple studies from different laboratories have generated RNA-seq data on isolated lens epithelium and fibers at different developmental stages (Table, Supplementary Table S1). To perform meta-analysis of these datasets, raw RNA-seq data from the SRA database was collected so that they can be analyzed using a similar pipeline (Fig. 1). The raw RNA-seq data were obtained through the SRA Toolkit (v2.10.9) and then mapped onto the mouse genome (GRCm38.p6) using the STAR software (v2.7.8a). Only the uniquely mapped reads were retained for further analysis. Differentially expressed genes (DEGs) were identified in R (v4.0.3) using the edgeR package (v3.32.1). Normalization of the samples was based on genes with expression above 0.2 counts per million (CPM). To control inter-sample relationships, principal component analysis and multidimensional scaling (MDS) plots were generated for samples from individual studies. DEGs were identified with the following thresholds: |log₂FC| > 1, false discovery rate (FDR) < 0.05, and an average expression in log₂ CPM > 1.0; where FC refers to fold-change and CPM refers to counts per million. Fragments per kilobase of transcript per million mapped reads (FPKM) scores were calculated using the CPM divided by gene length. Gene lengths were defined as the sum of exonic sequences for each gene, based on the Ensembl annotation (Mus_musculus.GRCm38.101.gtf). 

Table.

 

Isolated Lens Epithelium and Fiber Omics Datasets

 

Figure 1.

 

Experimental design for meta-analysis of isolated lens epithelium and fibers datasets. The flowchart outlines the experimental approach and pipeline for analysis of RNA-seq data from isolated mouse lens epithelium and fibers.

Analysis of Microarray Datasets on Isolated Lens Epithelium and Fibers

Two analyses have been conducted on the transcriptomic differences between lens epithelium and fibers—isolated by laser capture microdissection (LCM)—using microarrays at embryonic day (E) 12.5 and at postnatal day (P) 13.⁴¹ The E12.5 data was generated in Dr. David C. Beebe's laboratory as previously described.^14,44 Briefly, wild-type C57BL/6J mice were bred for timed matings for collection of E12.5 embryos. The embryonic tissue was frozen in OCT and sectioned. The sections were subjected to LCM to collect lens epithelium and fibers regions. Leica LMD 6000 laser microdissection system was used. Total RNA was isolated from the captured epithelium and fiber samples, reverse transcribed and amplified using the NuGEN Pico system as described.^14,44 The amplified cDNAs were analyzed on Illumina MouseWG-6 v2.0 Bead Chip whole-genome microarrays. With regard to the P13 data, a quality filter was applied to remove genes with probe tags labeled as NEG or MSR (negative or NA values). For the E12.5 data, only probes with a detection P value < 0.05 in at least one cell-type (either lens epithelium or fibers) were included in the analysis. Nonsignificant values were replaced with the minimum significant value for each cell-type to enable the calculation of fold-change between lens epithelium and fibers. For both datasets, gene expression was calculated as the logarithmic mean of all probes associated with each gene. Median signal normalization was performed to adjust the intensity values of the microarray data, using the median intensity of each array to reduce systematic variations between samples. Due to the lack of replicates, statistical significance calculations could not be performed. To highlight cell-type specific genes, only genes with |log₂FC| > 1 and those with average expression in the top 80% were retained. 

Analysis of Microarray Datasets on Isolated Lens Epithelium and Cortical and Inner Fibers

In addition to isolated lens epithelium and fibers, microarray datasets have been generated on the subregions of the mouse lens, namely, the epithelium (manually isolated), the inner cortex (isolated by LCM), and the cortical cortex (isolated by LCM).¹⁸ To compare these different regions, only probes with a detection P value < 0.05 in at least one sample were retained. Gene expression was calculated as the log mean of all probes associated with each gene. The R package edgeR was used to normalize the different samples. Principal component analysis and MDS plots were generated to control inter-sample relationships. The expression values were converted to log expression signals, and the thresholds used to define DEGs were |log₂FC| > 1, FDR < 0.05, and an average expression in log₂ signal > 0. 

In Silico Subtraction Analysis of Isolated Lens Epithelium and Fibers Datasets

To establish an enrichment score for each gene in either lens epithelium or fibers, the transcriptomic data from whole embryonic body tissue (WB) in iSyTE was used. We used RNA-seq and microarray data generated on the different platforms.^19,20 WB scores were normalized to the lens epithelium and fiber data using median signal normalization. The enrichment score was calculated as the log (expression score of a specific cell type/expression score of the WB). A gene was considered enriched if the log enrichment score was greater than 1, indicating it was expressed at least twice as much in the specific cell type compared to WB. 

Analysis Using the Cat-Map Database

To identify genes already associated with cataract due to mutation or mis-regulation, genes in the Cat-Map database⁴⁵ were used in comparative analysis. Human genes (n = 454) listed in the CatMap database (vOct 21) were compared with the mouse gene identified in the present analysis. 

Gene Ontology Term Enrichment and Pathway Analysis

The R package ClusterProfiler (v3.18.0) was used to identify gene ontology (GO) terms enriched in the cell-type-specific DEGs through GO enrichment analysis. This analysis included GO biological process (BP), GO cellular component (CC), and GO molecular functions (MF). 

Analysis Using Comprehensive Multi-Omics Platform for Biological InterpretatiOn

To identify relevant biological themes in the meta-analyzed transcriptome data from isolated lens epithelium and fibers across the different developmental and aged stages (E12.5, E14.5, E16.5, E18.5, P0, P13, three months, six months, and two years), we used the web-based tool called Comprehensive Multi-omics Platform for Biological InterpretatiOn (CompBio) developed at Washington University (CompBio platform (GTAC@MGI, WashU School of Medicine, https://gtac-compbio-ex.wustl.edu—Academic/Non-Profit; PercayAI Inc, www.percayai.com/compbio—Commercial).^46,47 For a specific list of genes (e.g., a “test” gene list of differentially expressed mRNAs (i.e., “entities,” numbering between 5 and 2500) in a particular cell/tissue), CompBio applies contextual language processing and biological language dictionary to identify relevant concepts from the published literature (e.g., >30 million Pubmed abstracts and >3 million full-text articles). This allows CompBio to perform de novo analysis independent of KEGG or Gene Ontology.^48,49 CompBio informs on the biological “concepts” that are significantly enriched in the test gene list, compared to a random list of genes. Enriched concepts that are commonly identified together are grouped into “themes.” Themes that have entities in common are further connected by edges to form co-occurring networks. These networks in turn inform on “central ideas.” Enrichment of concepts in a given theme (over randomly generated themes) is denoted by an “enrichment score.” This enrichment score can be used to further calculate a normalized enrichment score (NEScore) for a given theme. P values provided confidence for the significant enrichment of a given theme over random noise. In each dataset (e.g., E12.5, E14.5, etc.), genes that had higher expression in either epithelium or fiber cells (log₂ fold-change ≥ 2) were input into CompBio. An NEScore ≥ 1.3 and an associated P value < 0.1 were considered significant for theme identification. The identified themes were annotated to highlight key pathways and enriched entities. CompBio also provides auto-annotation of themes that were adopted for some themes. 

Results

Meta-Analysis of Isolated Lens Epithelium and Fiber RNA-seq Transcriptome Data

We sought to examine the transcriptomes of isolated lens epithelium and fibers generated in various studies performed on different stages, spanning embryonic and early postnatal development, through adulthood and aging. We identified several publicly available transcriptome datasets generated on microarrays or by RNA-seq that represent stages E12.5, E14.5, E16.5, E18.5, P0a (Cvekl data), P0b (Robinson data), P13, age three months, six months, and two years that could be used in this meta-analysis (Table, Supplementary Table S1). Of these, two datasets were on microarray platforms (E12.5, P13) whereas eight were generated by high-throughput RNA-seq. Meta-analysis was performed on these eight RNA-seq lens expression datasets (E14.5, E16.5, P18.5, P0a, P0b, age three months, six months, and two years). We first focused on meta-analysis of the isolated lens epithelium and fiber cell RNA-seq data by following the outlined workflow for each of these datasets (Fig. 1). This is important because the application of the same workflow to reanalyze these datasets serves to minimize variables that arise from the different analytical methods used in the original reports. This now allows for comparable downstream studies between these datasets. This RNA-seq meta-analysis identified 15,411 genes to be expressed in the lens epithelium or fiber cells in at least one stage (log₂CPM >0) (Fig. 2A). We next identified 2841 genes with robust expression (log₂CPM >1) that were also significantly differentially expressed (FDR <0.05) in the lens epithelium or fiber cells across all developmental stages (E14.5, E16.5, P18.5, P0a, P0b, age three months, six months, and two years) (Fig. 2A, Supplementary Table S2). A heat-map representation of the differential expression of these 2841 candidates shows separation of epithelium and fiber expressed genes (Fig. 2B). Furthermore, these differentially expressed genes also appropriately segregate according to the age of the lens, regardless of the origin of the sample (Fig. 2B). For example, the three-month data (generated in the Duncan laboratory) and six-month data (generated in the Fan laboratory) cluster closer compared to the three-month and two-year data, both simultaneously generated in the same laboratory (Duncan laboratory). Furthermore, the embryonic/early postnatal data (E14.5 through P0; generated by two different laboratories, namely, Cvekl and Robinson) cluster together while segregating separately from the adult/aged data (three months through two years). We examined the 2841 significantly differentially expressed genes across all the stages (fold-change (|log₂FC|>1)) where candidates showed distinct differential expression (e.g., higher in epithelium compared to fiber and vice versa) or dynamic differential expression (e.g., higher in epithelium compared to fiber and vice versa at certain stages, but exhibiting an opposite pattern in other stages). This identifies 1145 genes as higher in epithelium and 1427 genes as higher in fiber (Fig. 2A). This also identifies 172 genes that do not pass the stringent (|log₂FC|>1) fold-change threshold. These genes are differentially expressed to a lower extent; they are still significant in terms of statistical power for differential expression at all the stages. Finally, 97 genes are found to exhibit dynamic differential expression in epithelium and fiber cells at different stages. A few representative genes that exhibit these trends are shown (Fig. 2C). 

Figure 2.

 

Meta-analysis of lens-isolated epithelium and fiber RNA-seq datasets. (A) RNA-seq meta-analysis identifies 2841 robustly expressed genes (log₂CPM > 1, FDR < 0.05), which are differentially expressed between lens epithelium and fibers. Expression of significant genes (n = 2841) between lens epithelium and fibers across all the stages. This figure illustrates the distribution of the 2841 significantly DEGs between mouse lens epithelium and fiber cells across all developmental stages: (i) 1145 genes are consistently upregulated in epithelium, (ii) 1427 genes are consistently upregulated in fibers, (iii) 97 genes show dynamic differential expression between the two cell-types at different stages, and (iv) 172 genes did not meet the stringent fold-change threshold we set (i.e., they were differentially expressed to a lower extent), but these genes are still significant in terms of statistical power for differential expression. (B) Heat-map based on significant DEGs (n = 2841) in the lens epithelium or fibers across all developmental stages. P0a and P0b represent data from Cvekl and Robinson laboratories, respectively. (C) Examples of expression—across the different stages—of genes that are abundant in epithelium compared to fibers (e.g., Aqp1, Gja1), abundant in fibers compared to epithelium (e.g., Dnase2b, Gja3, Mip), initially abundant in epithelium and later in fibers (e.g., Gcnt2), initially abundant in fibers and later in epithelium (e.g., Gja8), and similarly abundant in epithelium and fibers (e.g., Yap1).

GO-Based Validation of Epithelial and Fiber Genes in Meta-Analyzed Data

We next performed GO analysis to validate the nature of the cell-type-specific expression represented in these datasets. We examined 1145 genes that are expressed higher in the epithelium and 1427 genes expressed higher in fibers for this analysis. In the epithelial cell data, the following GO terms were enriched: “morphogenesis of an epithelium,” “tissue morphogenesis,” “epithelial tube morphogenesis,” “regulation of cellular response to growth factor stimulus,” and more (Fig. 3A, Supplementary Table S3). In the fiber data, the following GO terms were enriched: “lens development in camera-type eye,” “lens fiber cell differentiation,” “mitochondrion organization,” “regulation of extrinsic apoptotic signaling pathway via death domain receptors,” and more (Fig. 3B). Together, these findings suggest that the meta-analysis effectively identifies dynamic and differential gene expression patterns in epithelial and fiber cells across the various embryonic, early post-natal, adult, and aged lens stages. 

Figure 3.

 

GO analysis of significant differentially expressed genes in isolated lens epithelium and fibers. (A) Top enriched GO terms for genes with higher expression (n = 1145) in lens epithelium. (B) Top enriched GO terms for 1427 genes with higher expression (n = 1427) in lens fibers.

Meta-Analyzed Epithelial and Fiber Datasets Identify Cataract-Associated Genes in Cat-Map

Next, we examined whether the meta-analyzed datasets could readily identify genes that are linked to cataract. Therefore we examined genes in the Cat-Map database, which is a resource for inherited and age-related cataract in human and animal models.⁴⁵ A majority of the cataract-linked genes (∼78%, n = 355 out of total 454) were identified in the meta-analyzed datasets (Supplementary Table S4). Of these, 37% were higher in epithelium whereas 34% were higher in fibers in at least one stage. Further, 6% were higher in epithelial or fibers depending on the specific stage. Moreover, 23% genes did not pass the fold-change threshold (|log₂FC|>1) in any of the different stages. This indicates that the meta-analyzed datasets are a useful resource for gaining insights into the specific cell-types that express cataract-linked genes. 

Identification of New Candidates in Lens Epithelium and Fibers

We next sought to identify new candidates and validate their expression in epithelium and fibers. Our analysis identified Btg1 (BTF anti-proliferation factor 1) and Igfbp5 (insulin-like growth factor binding protein 5) among the genes that exhibit higher expression in the epithelium compared to fibers (Supplementary Table S2). Furthermore, this analysis also identified Cpeb2 (cytoplasmic polyadenylation element binding protein 2), Id2 (inhibitor of DNA binding 2), and Pgam2 (phosphoglycerate mutase 2) among the genes that exhibit higher expression in fibers compared to the epithelium. RNA in situ hybridization (ISH) was used to validate these expression patterns. Furthermore, ISH provided insights into the subtleties of expression in epithelium and fibers. For example, Igfbp5 expression is restricted to anterior-most region of the epithelium, whereas Btg1 exhibits a more expanded expression pattern in the entire epithelium, as well as in the transition zone (Figs. 4A, 4B). Conversely, Cpeb2, Id2, and Pgam2 are all expressed in fibers, and although some are also expressed in the epithelium (e.g., Pgam2), they are still found to be higher in fibers (Figs. 4A, 4B). Thus ISH analysis independently validates localized expression of new candidate genes in lens biology. 

Figure 4.

 

Validation of meta-analysis by in situ hybridization. The gene expression predictions derived from the meta-analysis were independently validated through in situ hybridization. (A) In situ hybridization confirms the spatiotemporal expression of candidate genes as predicted by the meta-analysis. For instance, Igfp5 is expressed in the lens epithelium, Btg1 shows higher expression levels in the lens epithelium compared to lens fiber cells, Cpeb2 is more abundant in the transition zone and early differentiating fiber cells, Id2 is prominently expressed in the early differentiating fiber cells and Pgam2 is expressed in both lens epithelium and lens fiber cells. (B) Expression of these genes in epithelium or fibers based on RNA-seq meta-analysis at embryonic, adult and aged postnatal stages is shown. Intensity of heat-map indicates gene expression levels at different stages. e, epithelium; fc, fiber cells; tz, transition zone.

Epithelial and Fiber Expression Dynamics Across Age

Next, from the 15,411 genes that are found to be expressed in the lens epithelium or fiber cells in at least one stage at a cutoff of log₂CPM >0, we sought to identify genes that exhibit dynamic changes in expression pattern with age. Therefore, first, a total of 12,955 genes that exhibited higher expression in either epithelium or fibers in at least one of the eight time-points, by meeting the combined criteria of |log₂FC|>1 and log₂CPM>1 and FDR <0.05, were selected for further examination. This analysis identified 119 genes that were significantly expressed higher in the epithelium in embryonic development but later were significantly higher in fibers in adult stages. Interestingly, these genes are enriched in the GO term “biological process: hydrogen peroxide catabolic process,” which could suggest their importance for the resistance of oxidative stress in aged lenses. Conversely, 118 genes were found to be expressed higher in fibers in embryonic stages, but were found to be higher in the epithelium in adults (Supplementary Table S5). These genes exhibit interesting GO term enrichment, such as “CC: synaptic membrane,” “CC: contractile fiber,” and “CC: extracellular matrix,” which may reflect the cellular functions such as elongation and capsule formation. 

WB In Silico Subtraction Identifies Genes With Enriched Expression in Epithelium or Fibers

Previously, we showed that “WB-in silico subtraction” effectively identifies genes that have enriched expression in a specific cell/tissue-type compared to the whole embryonic body tissue.¹¹ In the past, WB-in silico subtraction has been applied to whole lens datasets to identify new genes and pathways linked to lens biology and cataract.^1,50 Although this has been effective, it does not inform on the genes that exhibit enriched expression specifically in the epithelium or fibers. Therefore we performed WB-in silico subtraction on the meta-analyzed data, which leads to identification of genes with elevated expression in the specific cell/tissue type compared to WB (log₂ fold-change (expression in lens cell-type/expression in WB) >1). This identified 9312 genes to exhibit enriched expression in the epithelium (in at least one stage across all the ages) and 9580 genes to exhibit enriched expression in the fibers (in at least one stage across all the ages). The heatmap figure was generated by retaining the expression values of genes across all time points (n = 11,256) (Fig. 5A). Unsurprisingly, a large number of genes are enriched in both fibers and epithelium (compared to WB), with the highest enrichment scores being exhibited by genes encoding crystallins (e.g., Cryaa, Crybb, etc.) (Supplementary Table S6). Notably, although the enrichment levels of crystallins are lower in the epithelium compared to fibers, their overall higher expression compared to WB results in these genes being called as exhibiting the high enrichment in the epithelium. To effectively identify cell-type-specific enriched genes, a cell- type enrichment difference (Delta) was calculated (Delta = enrichment score in fibers minus enrichment score in epithelium). The Delta value highlights genes preferentially enriched in specific cell-types (Figs. 5B, 5C). For example, Birc7, which is known to be highly expressed in fibers, shows the largest enrichment difference score (Fig. 5C). Other genes with known high expression in fibers, such as Dnase2b and various crystallins, are also identified among genes with the highest enrichment difference scores. Interestingly, there are several genes that have not yet been studied in the lens context that also show high enrichment difference scores, suggesting they may play critical roles in lens fiber differentiation and homeostasis. Conversely, genes with the negative delta enrichment scores suggest significantly enriched expression in the epithelium. Among the top genes with enriched expression in the lens epithelium is Gja1 that is known to be expressed in the epithelium, along with several other genes that have not yet been studied in the context of the lens (Fig. 5B, Supplementary Table S6). Also among the enriched expressed genes are many other new promising candidates for future studies in the lens. These include Ell2, and Prdm16, among others. 

Figure 5.

 

In silico-subtraction identifies genes with enriched-expression in isolated lens epithelium and fibers. In silico-subtraction allowed analysis of genes with respect to their enriched-expression compared to WB in lens-isolated epithelium and fibers (log₂ fold-change (expression in lens cell-type/expression in WB) >1) across the different stages. (A) Heatmap representing enriched-expression of genes in both lens-isolated epithelium and fibers compared to WB across all the stages. (B) Isolated lens epithelium-specific enriched expression of genes calculated by negative delta values and indicated by heatmap color gradients (left column: maximum enrichment score in fibers, right column: maximum enrichment score in lens-isolated epithelium, delta values indicated by separate heat-map). (C) Isolated lens fiber-specific enriched expression of genes as denoted by positive delta values.

Transcriptome Analysis of Isolated Lens Epithelium and Fibers at the Early Embryonic Stage of E12.5

In mouse development, soon after the lens placode (formed at E9.5) invaginates into a lens pit (E10.5), this structure transitions and progressively develops into a lens vesicle (at E11.5-E12.5). At this stage, morphological differences between the anteriorly and posteriorly localized cells can be readily distinguished. We next sought to examine the differences in the transcriptome between the anteriorly and posteriorly localized cells at this early stage (at E12.5), when primary fiber differentiation is commenced, which has never been addressed. We identified a publicly available microarray expression dataset for LCM mouse lens epithelium and fibers at E12.5 that was generated by Dr. David Beebe (Fig. 6A, Table, Supplementary Table S7). Analysis of this dataset identified 1618 genes with elevated expression in the epithelium and 1548 genes with elevated expression in fibers (Fig. 6B). GO term analysis shows that E12.5 epithelium is enriched for “DNA replication” and other cell proliferation related terms, while E12.5 fibers are enriched for genes involved in the biological processes for “structural constituent of eye lens,” “lens development in camera-type eye” and “camera-type eye development” (Figs. 6C, 6D). Indeed, several key genes with well-established expression in the epithelium (e.g., Cdh1, Pax6) or fibers (e.g., Cryaa, Crybb1, Crybb3, Crygc, Tdrd7) exhibit the expected cell-type expression (Supplementary Fig. S1). These analyses independently validate the effective isolation of the two different lens cell types at this technically challenging early stage. Furthermore, this data shows that even at this early stage, the posteriorly localized cells exhibit a transcriptome switch that will lead to primary fiber differentiation (Supplementary Table S7). 

Figure 6.

 

Analysis of microarray transcriptomic dataset of isolated-lens epithelium and fibers at E12.5 and P13. (A) Schematic of LCM of mouse lens epithelium and fibers at E12.5. Area marked by broken lines denote the tissue region that was isolated as epithelium (Epi) or fibers. (B) MA plot demonstrates relative expression of genes in lens-isolated epithelium (1618 higher in epi) versus fibers (1548 higher in fibers). (C) Gene ontology (GO) analysis for the 1618 genes elevated in E12.5 lens-isolated epithelium. (D) GO term analysis for the 1548 genes elevated in E12.5 lens-isolated fibers. (E) GO analysis for the 1532 genes elevated in P13 lens-isolated epithelium. (F) GO term analysis for the 1793 genes elevated in P13 lens-isolated fibers.

Transcriptome Analysis of Isolated Lens Epithelium and Fibers at Mid-Postnatal Stage P13

To identify the cell type–specific changes in mid-postnatal lens development, we analyzed a publicly available custom-microarray dataset on isolated epithelium and fibers at P13. Analysis of these data identified 1532 genes with elevated expression in the epithelium and 1793 genes with elevated expression in fibers. GO term analysis at stage P13 shows that the epithelium is enriched for “DNA replication,” while P13 fibers are enriched for genes involved in the biological processes for “Structural constituent of eye lens,” “lens development in camera-type eye” and “cell-cell junction” (Figs. 6E, 6F, Supplementary Table S8). This independently validates the cell-type specific expression in these datasets at P13. Furthermore, this data describes the transcriptome differences between epithelium and fibers at mid-postnatal stage (Supplementary Table S8). 

Analysis of Isolated Lens Epithelium and Outer Cortical and Inner Cortical Fiber Transcriptomes

All the datasets discussed above focus on isolated epithelium and fibers in bulk. However, fiber cells exhibit changes in gene expression based on their differentiation stage.^51,52 To gain insights into these nuanced gene expression changes from distinct fiber cell populations, we analyzed the transcriptome datasets from fibers isolated at two different “depths” in one-month-old mouse lens—from the outer and inner cortex—in addition to analyzing isolated lens epithelium datasets that are generated by Dr. Stephen Bassnett's laboratory and are publicly available¹⁸ (Fig. 7A, Table). This analysis identifies 1422 DEGs between one-on-one comparisons of isolated epithelium, outer cortex, and inner cortex (Fig. 7B, Supplementary Table S9). Only 284 DEGs are found to be common between all comparisons (Fig. 7C). Not surprisingly, comparison between the isolated epithelium and the inner cortical fibers (representing the two most different cell-types) accounted for the vast majority (94%) of DEGs (1349 out of total 1422 DEGs). Comparison of the outer cortical fibers with either the isolated epithelium or the inner fibers identified a relatively lower number of DEGs, 56% (796 of 1422) and 56% (795 of 1422), respectively. This is not a surprising finding because it suggests that the outer cortical fiber dataset is representative of the transcriptome progressing from an epithelial to fiber cell types. 

Figure 7.

 

Comparative transcriptomic analysis of isolated epithelium, outer cortical fibers and inner cortical fibers from 1-month old mouse lens. (A) Schematic representation of lens outer cortical fibers (green) and inner cortical fibers (red) cell populations captured by of LCM from one-month-old mouse lenses in a study published from the Bassnett laboratory. The epithelium (asterisk) from the one-month old mouse lenses was isolated by manual dissection. (B) Heat-map of significant DEGs (n = 1422) between isolated epithelium, outer cortical fibers, and inner cortical fiber datasets. (C) Venn diagram shows 248 common genes between the three datasets. (D) Comparison of DEGs from the one-month-old lens isolated epithelium, outer cortical fibers, and inner cortical fibers datasets and the meta-analyzed RNA-seq data on isolated epithelium and fibers. The pie-charts inform on the respective proportions that data from the Bassnett study (e.g., the epithelium (i.e., “Bassnett Epithelium”), outer cortical fibers (i.e., Outer Fibers), inner cortical fibers (i.e., Inner Fibers) falls into, as per the meta-analyzed RNA-seq data on isolated epithelium and fibers (given in the inset key).

Next, we compared the E12.5 and 1 month old datasets generated by the groups of Drs. Beebe and Bassnett to examine whether these data commonly identified a subset of genes expressed in the specific lens cell-type. Although these data are generated on different microarray platforms and at vastly differences stages, 562 genes were found to have similar differential expression in the different cell-types (epithelium vs. fibers [bulk or outer or inner cortical fibers]) (Supplementary Table S10). Out of these 562 genes, 155 were commonly identified in the epithelium, 215 were commonly identified in the fibers in the E12.5 and age one month datasets. Furthermore, 43 genes were identified in epithelium at E12.5 and then were identified in fibers at age one month, whereas 105 genes were identified in fibers at E12.5 and then were identified in the epithelium at age one month. Notably, 47 genes were specific to the outer cortical fibers, and among these, three were epithelial at E12.5, whereas the remaining 44 were identified in E12.5 fibers. This suggests that majority of the commonly identified genes in these two datasets retain their cell-type-specific expression, while a minority of genes show dynamic cell-type expression changes at these different stages. 

Next, we compared the 1422 DEGs in Bassnett datasets with the genes identified in these cell-types by meta-analyzed RNA-seq datasets. Expectedly, the majority of the 670 genes that are expressed in the epithelium in the Bassnett dataset (i.e., higher expressed in the epithelium compared to inner cortical fibers and outer cortical fibers) are also found to be expressed in the epithelium in the meta-analyzed RNA-seq datasets (Fig. 7D). A subset of genes was found to show differential expression with time or was recognized to be expressed in fibers or was not specific to any cell-type. Similarly, the vast majority of the 703 genes that are expressed in the inner cortical fibers in the Bassnett dataset (i.e., higher in the inner cortical fibers compared to outer cortical fibers and the epithelium) are also recognized as fiber expressed genes in the meta-analyzed RNA-seq datasets (Fig. 7D). However, smaller subsets of genes were found to show differential expression with time or were recognized to be expressed in the epithelium or were not specific to any cell-type. Because outer cortical fibers were between the epithelium and inner cortical fibers, we were interested in finding genes that are specifically expressed in this region. In the Bassnett dataset, 37 genes were found to have specifically higher expression in the outer cortical fibers compared to both epithelium and inner cortical fibers, suggesting that these genes were dynamically expressed in a spatial manner in the lens (Fig. 7D). Most of these 37 genes were recognized as fiber-expressed genes when compared to the meta-analyzed RNA-seq datasets. We were also interested in identifying genes that show a transient “dip” in expression between the epithelium and the inner cortical fibers. Our analysis found seven genes to have higher expression in the epithelium and inner cortical fibers but not outer cortical fibers, suggesting their dynamic control in the lens (Fig. 7D). Together, these findings suggest that cell-type–specific data in the Bassnett study largely agrees with the Beebe microarrays and RNA-seq meta-analysis data. 

Pathway Analysis of Epithelium and Fibers From Embryonic Through Aging

Next, we sought to gain biological insights from genes and pathways that are preferentially expressed in the lens epithelium or fibers at different stages. Toward this goal, we used CompBio (Comprehensive Multi-omics Platform for Biological InterpretatiOn), a web-based tool that generates maps of core biological concepts that are associated with omics-dataset inputs (e.g., a list of genes (termed “entities” in this analysis) that are expressed in a particular cell-type or condition as revealed by RNA-seq). In an ontology-independent manner, CompBio conducts de novo analysis of >30 million PubMed abstracts and >3 million full-text articles to identify themes that are associated with a specific input, as outlined below. A CompBio search for a given list of entities (genes) lead to the identification of individual “concepts.” Concepts that are frequently identified concomitantly are grouped into “themes.” Themes that share entities are connected by edges to represent networks that can be used to derive “central ideas.” We applied CompBio to analyze genes that were preferentially expressed in either the lens epithelium or fibers (log₂ fold-change ≥ 2) at different stages (e.g., E12.5 through 2 years). 

Distinct themes with significant enrichment were identified in the epithelium and fibers. Not surprisingly, several themes corresponding to cell proliferation (e.g., “DNA replication,” “cytokinesis, cell division,” “cell cycle regulation,” etc.) were identified in the epithelium. Other themes identified in the epithelium include “ECM” (extracellular matrix), “Ephrin signaling,” “Wnt signaling,” and several themes related to SLIT/ROBO signaling and axon pathfinding/guidance. (Fig. 8, Supplementary Fig. S2). Additionally, other notable enriched themes identified in epithelium were “Ectopia lentis” and themes pertaining to cell adhesion. In fibers, the themes identified include “Zonular cataract,” “Hsp27 protein binding,” and “Molecular chaperones and chaperonins.” Multiple other themes corresponding to protein synthesis, processing and stability were also identified in fibers (e.g., “Protein synthesis,” “Ribosomal proteins,” “Protein processing,” etc.). In some fiber stages, “Autophagy” as a theme was also identified. Notably, in fibers, themes related to RNA processing (“mRNA processing,” “mRNA processing and stability, nuclear import,” “AU-rich element binding; mRNA stability,” “Nonsense-mediated decay,” “RNA modification and regulation”) emerged as enriched themes specifically in all embryonic stages and at P0. This emphasizes the importance of post-transcriptional regulation and RNA-binding proteins (RBPs) in embryonic and early postnatal lens development. Indeed, fiber differentiation necessitates that several cellular challenges are met, wherein RBPs may have a functional role.^50,53 For example, select mRNAs (e.g., crystallins) are translated into unusually elevated levels of proteins. Furthermore, fiber cells need to optimize gene expression control taking into account the spatial context as they undergo ∼1000-fold length-wise increase while undergoing migration from outer to inner regions of the lens and prior to organelle degradation. Indeed, recent studies have defined the role of different RBPs (e.g., Caprin2, Celf1, Rbm24 and Tdrd7) in post-transcriptionally regulating several important targets (e.g., Cdkn1b (p27^Kip1), several crystallins, Dnase2b, Hspb1, Pax6, Prox1) involved in fiber differentiation.^{10,22,27–30,33,54–57} Additionally, in both epithelium and fibers, at least one theme related to lipid metabolism and signaling (e.g., “lysophosphatidic acid signaling,” “lipid signaling,” “phospholipid metabolism,” “lipid metabolism,” phosphoinositide metabolism,” “S1P signaling,” etc.) was identified across all developmental stages, highlighting the importance of lipids in the lens. Similar themes were highlighted when ranked on the basis of Normalized Enrichment Scores (NEScore) (Fig. 9, Supplementary Fig. S3). Together, these data uncover key themes associated with gene expression in lens epithelium and fiber cells across different stages of development and aging. 

Figure 8.

 

The bioinformatics tool CompBio identifies biological themes enriched in isolated lens epithelium and fibers. CompBio-derived map of distinct biological themes enriched for genes with higher expression in isolated lens epithelium or fibers at different stages, namely, embryonic day (E) 14.5, postnatal day (P) 0, three months, and two years. Note that P0 represents data from the Cvekl laboratory. Similar colors for related themes are shown across the different stages.

Figure 9.

 

CompBio-based analysis of isolated lens epithelium and fibers using Normalized Enrichment Score (NEScore). The top ten enriched biological themes based on NEScore rank in isolated lens epithelium (Epi) and fibers at different stages, namely, embryonic day (E) 14.5, postnatal day (P) 0, three months, and two years. Note that P0 represents data from the Cvekl laboratory.

Furthermore, the entities (in this case, genes preferentially expressed in specific cell-types) that contributed majorly to the identification of the different concepts—and thereby the themes—across multiple stages in the epithelium were Cdh1, Cdk1, Gja1, Pole, Trpv4, and more, whereas those in the fibers were Aqp5, Bfsp2, Cryaa, Hspb1, Prkaa2, Sptan1, Tjp1, and more (Fig. 10A, Supplementary Table S11). 

Figure 10.

 

CompBio-based identification of key genes and central biological ideas in isolated lens epithelium and fibers across embryonic, postnatal, adult and aged stages. (A) The top 10 genes based on their contribution to the prediction of distinct themes in isolated lens epithelium and fibers at different stages, namely, embryonic day (E) 12.5, E14.5, E16.5, E18.5, postnatal day (P) 0a, P0b, P13, three months, six months, and two years. P0a and P0b represent data from Cvekl and Robinson laboratories, respectively. (B) Central ideas highlighting key biological processes in isolated lens epithelium and fibers at the different stages stated above.

Not surprisingly, the central ideas (i.e., biological insights) in the epithelium were associated with cell proliferation and extracellular matrix (Fig. 10B, Supplementary Fig. S4). Some central ideas observed in later stages in the epithelium were those associated with channel proteins and retinal photoreceptors (three months through two years). In fibers, central ideas that are associated with lens development, cytoskeleton, fiber cell apoptotic process, and DAG signaling were identified in several stages (Fig. 10B, Supplementary Fig. S4). Interestingly, in later stages, fibers showed central ideas associated with chaperones/heat shock proteins and ribosomal subunits. Thus CompBio analysis provides key biological insights into the dynamics of epithelial and fiber-enriched gene expression in the lens at embryonic, postnatal, adult and aged stages. 

Accessing the Isolated Lens Epithelium and Fiber Data in iSyTE

To make the meta-analyzed data on isolated lens epithelium and fibers readily accessible, we developed a new user-friendly web portal on the iSyTE webpage that can be navigated at http://research.bioinformatics.udel.edu/iSyTE. Examples are provided below to describe the utility of the new web features in iSyTE for fast and effective visualization of the meta-analyzed data. One or more genes can be examined through the new web portal at iSyTE. For example, one can navigate to the iSyTE webpage, select “Gene Expression,” “Epi Vs Fib,” Select dataset: “FPKM Table (RNA_Seq)” (e.g., for RNA-seq data), Select category: “All” or “Epi” or “Fib,” and search by gene symbols (Fig. 11A). Examples are shown for genes with higher expression in lens epithelium or lens fibers at the different stages. Namely, Cdh1 and Foxe3 exhibit higher expression in epithelium compared to fibers (Fig. 11A), while Bfsp2 and Mip exhibit higher expression in fibers compared to epithelium as expected based on their established expression pattern (Fig. 11B). Additionally, the new iSyTE portal also allows examination of epithelium or fiber gene expression based on (1) “Expression (Microarray)” (at E12.5, P13), (2) “Enrichment vs WB (Microarray + RNA_Seq),” (3) “Fold Change Epi vs Fib (Microarray + RNA_Seq), (4) “Enrichment vs WB (RNA_Seq),” and (5) “log₂FC Epi vs Fib (RNA_Seq).” The iSyTE 2.0 even captures a well-established expression switch in the TFs Sox2 and Sox1 (Fig. 11C), indicating the sensitivity of the meta-analyzed data. Together, these findings demonstrate that the meta-analysis and the new web-portal-enabled data visualization in iSyTE allows effective examination of gene expression in normal lens epithelium or fibers in embryonic, adult and aged stages. 

Figure 11.

 

Access of isolated lens epithelium and fiber gene expression in iSyTE. (A) The meta-analyzed data on isolated lens epithelium and fibers is made freely accessible in a new user-friendly web portal on the iSyTE webpage at http://research.bioinformatics.udel.edu/iSyTE. This provides effective visualization of single or multiple genes at the different stages by following Steps 1 through 5. Examples are provided of genes with higher expression as per the meta-analyzed RNA-seq-based FPKM data in the lens epithelium (Cdh1, Foxe3), (B) lens fibers (Bfsp2, Mip), (C) transcription factors with role in the lens (Sox2, Sox1), and (D) candidates linked to cataract from genetic association studies (Dpy19l3, Herc4).

The iSyTE Expression Analysis of Candidate Genes Identified by GWAS and Transcriptome-Wide Association Studies for Cataract

In the past, information on expression in the lens has been helpful in independently supporting gene candidates identified by genome-wide association studies (GWAS)^31,58 or by multi-tissue transcriptome-wide association studies.⁵⁹ However, this relied on whole lens transcriptome data. The meta-analysis in this study now allows cell-type expression data to be analyzed for these candidates as shown by a few examples (Fig. 11D). These data will help to develop hypotheses regarding the specific function of these genes in the lens. 

Discussion

Although effective in discovery of genes linked to lens development and cataract, the earlier version of iSyTE—which predominantly currently represents data generated on whole lenses—does not inform on gene expression in lens epithelium or fibers. The present meta-analysis study was undertaken to address this gap. These efforts identified 2841 significantly differentially expressed genes across all the stages, namely, those higher in epithelium (n = 1145 genes), higher in fibers (n = 1427 genes), or dynamically expressed in epithelium and fibers (n = 97 genes). Furthermore, we identified genes (n = 172) that, although being significantly differentially expressed, exhibited this pattern to such a low extent that did not meet the stringent criteria which was set to identify only the top preferentially expressed genes in the epithelium and fibers. These lower differentially expressed genes could still be of interest in the lens. Furthermore, in early lens development stages, posteriorly located cells in the lens vesicle that are derived from the lens ectoderm go on to differentiate into primary fiber cells. There has been no omics level analysis of isolated epithelium and isolated fibers at an early stage such as E12.5 when the two cell-types can be isolated with confidence—which the present meta-analysis study also addresses. 

In addition to informing on epithelium and fiber gene expression and its utility as a tool in iSyTE for discovery of cataract-linked genes, the meta-analysis serves to provide new insights into the biology of the lens across multiple developmental and aging stages. The tool CompBio provides new hypothesis-generating insights in the lens. For example, CompBio analysis shows that aging lens epithelium (stages three months, six months, two years) exhibits expression of genes normally found to be expressed in retinal photoreceptors. This suggests that the developing/early lens likely has mechanisms to repress retina/photoreceptor fate—mechanisms that are likely compromised with age. Indeed, previous studies have suggested that the calcium-binding protein S100A4 (gene deletion of which causes misexpression of retina genes in the lens⁶⁰), as well as the Polycomb repressive complex 2 (PRC2) epigenetic regulators Ezh2 and Jarid2 (which are localized to retina transcription factors gene loci that are normally repressed in the lens³⁵), are likely involved in repressing a retina differentiation program in the lens. Indeed, our meta-analysis shows that expression of Ezh2 and Jarid2 is reduced in the epithelium as the lens ages. These analyses lead to an exciting hypothesis that can be tested in the future, namely, whether lens-specific conditional knockout of Ezh2 and Jarid2 results in misexpression—in early stages of the lens—of genes involved in retina differentiation. 

Pathway analysis by CompBio also gave insights with regards to lipid signaling and lipid metabolism in lens epithelium and fibers. There is evidence that lipid and sterol composition are critical for lens transparency.^61,62 Indeed, lanosterol synthase (LSS) gene mutations are linked to congenital cataract in human⁶³ and rat,⁶⁴ mutation in the human SC4MOL gene, which encodes a methyl oxidase necessary for cholesterol synthesis, are linked to syndromic cataract,⁶⁵ and several sterol pathway genes are misexpressed in Mafg^−/−:Mafk^± compound mice that develop cataract.²³ The present meta-analysis shows sterol pathway genes (e.g., Fdft1, Fdps, Hmgcr, Idi1, Lss, Mvd, Mvk, Pmvk, Sqle, etc.) to be expressed higher in fibers compared to epithelium. Additionally, LPA (Lysophosphatidic acid signaling) and S1P (Sphingosine-1-phosphate signaling) are important for cell proliferation and cell migration. CompBio shows that the receptors for these pathways (e.g., Lpar1, S1pr3) are more highly expressed in the epithelium compared to fibers, suggesting their role in the cellular properties (e.g., proliferation) of the epithelium. 

Additionally, CompBio prioritized several interesting “entities” (genes) for further investigation in the lens. For example, CompBio prioritized the cytoskeletal protein Sptan1 (Spectrin alpha chain, non-erythrocytic 1) in fibers. Sptan1 protein was also detected in the top 150 lens-expressed proteins and top 150 lens-enriched proteins in the E14.5 mouse lens.⁶⁶ SPTAN1 was identified among the genes associated with human congenital cataract by exome-sequencing.⁶⁷ Furthermore, a SPTAN1 single nucleotide polymorphism variant has been associated with age-related cataract in humans.⁶⁸ Interestingly, initial findings from International Mouse Phenotyping Consortium (IMPC, www.mousephenotype.org), also suggest that Sptan1-deficient mice exhibit ocular defects including microphthalmia and cataract. Among other genes prioritized by CompBio were Gfap (Glial fibrillary acidic protein), Kdr (Kinase insert domain receptor) for epithelium, and Avp (Arginine vasopressin), Dnm3 (Dynamin 3), Prkaa2 (Protein kinase AMP-activated catalytic subunit alpha 2) for fibers. 

Furthermore, the meta-analysis also prioritized signaling pathways and their components for further investigation. For example, the Eph-ephrin signaling pathway—previously known to be important for lens development and transparency^69–75—was also found to be enriched as a theme in the epithelium at different stages. In addition to Epha2 (receptor), Efna5 (ligand), and the recently predicted³² Epha5 (receptor), the present analysis also highlighted Efna1 (ligand), Efna4 (ligand), Efnb3 (ligand) as candidates that may be relevant to lens biology. Additionally, components of SLIT/ROBO signaling and axon pathfinding/guidance were highlighted in the lens epithelium datasets. Interestingly, SLIT/ROBO signaling is shown to be involved in nerve repulsion during cornea innervation⁷⁶ and therefore its role in the lens may need to be further investigated. Furthermore, components of Wnt signaling (e.g., Fzd2, Fzd7, Wnt5a, etc.) were also identified among the themes enriched in the epithelium. In fibers, all embryonic and early postnatal stages showed enrichment of certain themes related to post-transcriptional gene expression control. These include mRNA processing, nuclear transport, stability, or AU-rich element binding. This is significant to lens biology, because since as long as 60 years ago, studies have postulated the involvement of post-transcriptional gene expression control in the lens, based on initial observations that for certain genes, the mRNA and the encoded protein levels did not correlate, or that epithelial and fiber cells differ in the stability of specific mRNAs.^50,77–85 Furthermore, fiber cells need to translate specific mRNAs into proteins at unusual abundance (e.g., crystallins) while also ensuring that other key regulators, the mRNA levels of which are comparatively low (e.g., c-Maf, Mafg, Mafk, Prox1, Sox1, etc.), are also made in optimal quantities. Recent studies have shown that post-transcriptional regulation mediated by RBPs (e.g., Caprin2, Celf1, Rbm24, Tdrd7, etc.) is necessary to achieve optimal spatiotemporal control over the lens proteome and that their loss/alteration causes lens defects and/or cataract.⁵⁰ Thus, in addition to these few established RBPs, it will be interesting to examine new RBPs (e.g., Cpeb3, Cpeb4, Msi1, Rbm38, Zfp36) that may be involved in these processes. 

This is the first report of a detailed analysis of an early stage (i.e., E12.5) when the lens epithelium and the fibers can be reliably isolated. CompBio analysis identifies themes that correlate with epithelium or fiber biology. For example, in the E12.5 epithelium, the following themes are enriched: proliferation, ECM, chromatin remodeling complex, and histone modification. In the E12.5 fibers, the following themes are enriched: autophagy, regulation of actin cytoskeleton, cell adhesion and junctional complexes, lipid metabolism, regulation of protein stability, ubiquitination and protein degradation, and zonular cataract. 

Furthermore, GO analysis of these meta-analyzed datasets independently highlights critical components related to cellular division, ECM organization, and signaling pathways in the lens epithelium, and those related to cellular stress responses and protein structure in lens fibers, which is consistent with the properties of these cells. For example, GO analysis reveals an interesting finding that several genes associated with mitochondrial function are enriched in fiber cells. It can be speculated that this may in turn reflect the physiological demand of the early differentiating cells, for example in building up the high amounts of lens transcripts and proteins (crystallins, membrane proteins, etc.) before reaching the terminal differentiation stage. We note that several genes involved in mitochondrial function (e.g., Gpx1, Fundc2, Ndufa2, Ndufb10, etc.) that were found to be upregulated in chicken lens fibers in a previous study⁸⁶ are also found to be expressed higher in mouse fibers in the present analysis. 

Additionally, WB in silico subtraction prioritizes many interesting genes for future studies in the lens. For example, this approach prioritizes factors with known involvement in transcription, such as, Ell2 (transcription elongation regulator) and Prdm16 (transcription factor). Interestingly, Ell2 enriched-expression is independently found to be conserved in fish lens (Zfin database, https://zfin.org/ZDB-GENE-030131-1717/expression) and Prdm16 expression is reduced in established gene perturbation models of cataract (e.g., Celf1^cKO lens), further supporting their examination in the lens. Although considerable progress has been made on transcriptional control in the lens, finding new transcriptional regulatory candidate genes in the lens will allow new insights in future studies, thus advancing our understanding of how specific genes (e.g., crystallins) are expressed at high levels in lens fiber cells. There are also many other promising candidates potentially involved in a variety of different functions in the lens. These include Casz1 (potential transcription regulator also found to be associated with human cataract in GWAS), Htra3 (HtrA serine peptidase 1), Msi1 (RNA-binding protein), Rmst (long non-coding RNA), Rnf123 (ring-finger protein), Scel (sciellin), Zbtb8b (potential transcription regulator also prioritized by previous versions of iSyTE), among many others. Below we discuss a few of these candidates in the context of lens biology. 

The present meta-analysis has also identified potential new genes with relevance to lens biology. For example, Rmst, which encodes a long non-coding RNA and exhibits higher expression in fibers, may potentially contribute to Sox2’s regulatory impact in the lens, because it is known to regulate the recruitment of this transcription factor to its genomic targets.⁸⁷ Control of p27^Kip1 levels in the lens is important as it determines lens epithelial cell-cycle exit and initiation of fiber differentiation⁸⁸ as well as nuclear degradation in fiber maturation.^27,89 To expand the list of known regulators (e.g., Celf1) of p27^Kip1,²⁷ the present study identifies Rnf123 as a potential new regulator of p27^Kip1 in the lens. This is because Rnf123 encodes a protein with E3 ubiquitin ligase activity that is known to negatively regulate p27^Kip1,⁹⁰ and based on our findings Rnf123 exhibits high expression in fibers and thus is an excellent candidate for functional examination in the lens. 

Additionally, our study identifies new candidates such as Htra1 and Htra3 that encode serine proteases and exhibit dynamic expression in epithelium (Htra3 high in earlier stages and Htra1 high in aged stages in lens epithelium) and may be involved in controlling the extracellular matrix. Further, this metanalysis also identified new candidates such as Scel, whose expression is progressively elevated with age in the lens epithelium. Scel is involved in epithelial to mesenchymal transition (EMT),⁹¹ and therefore may be relevant to studies examining the factors that contribute to posterior capsular opacification (PCO) as a complication of cataract surgery. This analysis also shows that eIF5B exhibits progressively high expression in fibers through the different stages. Recent studies have shown a role for eIF5B in negatively regulating apoptosis by promoting translation of anti-apoptotic factors encoded by IRES-containing transcripts (e.g., Bcl2l1, Cflar (c-FLIPS)),⁹² which as our study shows, also exhibit progressively high expression in fiber cells. This may be relevant to the biology of the aging lens because it is shown that reduced eIF5B results in elevated apoptosis by reactive oxygen species.⁹² Thus the present meta-analysis identifies many new candidates genes for further investigation in various aspects of lens biology including its development, homeostasis, aging, and pathology. 

In addition to the various above approaches identifying new genes, this study allows us to gain new comprehensive information on the spatiotemporal expression dynamics of key transcription factors in the lens from embryonic, early post-natal, adult and aging stages. For example, Pax6, Foxe3 and Tfap2a exhibit high expression in epithelium, whereas Maf, Mafg, and Sox1 exhibit high expression in fibers, at all these stages. In contrast, Sox2, which exhibits high expression in epithelium at embryonic and early postnatal stages, shows close to uniform expression in these cell types (with slightly higher expression in fibers) at adult and aging stages. These data give new insights to examine a potential new role of Sox2 in adult/aging stages in lens fibers. 

Among its many downstream applications, the present meta-analysis allows us to correlate the findings from GWAS for cataract by examining the expression profiles of GWAS hit loci/genes in lens cells at embryonic through aging stages, to gain more insight into their relevance to lens biology. For example, a recent such comparative study identified, among other targets, the genes DPY19L3 and HERC4, to be associated with cataract in humans.⁹³ Dpy19l3 encodes a c-mannosyltransferase that mannosylates a WNT regulator and is expressed higher in epithelium compared to fibers in aged lenses, suggesting its potential function (and potentially also that of other Wnt associated/target proteins, e.g., Lrp5, Lrp6, and Camk4, which also exhibit higher expression in lens epithelium compared to fibers) in the aging lens epithelium. Herc4 encodes a ubiquitin ligase that is abundant in fibers compared to the epithelium in all stages, suggesting a role in perhaps post-translational modification of proteins in these cells. Thus both genes are good candidates for future examinations in the lens. 

The primary goal of the present research was to meta-analyze and make readily available all the omics-level datasets on isolated epithelium and fibers that are currently deposited in GEO. We had no control over the choice of the mouse strains (e.g., C57BL/6J, CD-1, FVB/N) or the specific omics-platforms (i.e., RNA-seq vs. microarrays, etc.) used for generating these data, which were the decisions of the PIs who generated the data. Although different strains may likely harbor differences in expression of specific genes, these do not seem to have a large impact on the global gene expression profiles on these cell-types. This is because our analysis shows that regardless of the laboratory/strain, the overall expression of genes correlated closely with the stages/age and not with the specific laboratory or the strain that the data was generated on. For example, in Figure 2B, the expression datasets generated at stage P0 from the Cvekl laboratory (generated on CD-1 background) and Robinson laboratory (generated on FVB background) correlate closely, compared to the correlation of P0 data generated by the Cvekl laboratory with data at E14.5 and E16.5 generated on same strain (generated on CD-1 background by the Cvekl lab) —demonstrating that regardless of the lab or the strain, overall the data correlate with the specific age. Similarly, at later stages, data made on C57BL6 background at age three months (Duncan laboratory) and at age six months (Fan laboratory) correlates closely, compared to the correlation of the data at three months (made by the Duncan laboratory) with data at age two years (also made the Duncan laboratory)—demonstrating that regardless of the laboratory source, overall the data correlate closer with the stage. However, there is a possibility that expression differences in specific genes could arise from strain/background differences. For example, it is possible that changes observed between embryonic/early postnatal lens datasets (generated on CD-1, FVB mouse backgrounds) and aged (generated on C57BL/6J mouse background) lens datasets could partially arise from the differences between different strains. Thus it is recommended that in future studies, expression should be independently validated for candidate genes in the specific mouse strain model to be used. 

In sum, this study gives new insights into the biology of lens epithelium and fibers and provides a new web-based interface for the examination of gene expression in these cell-types. These new updates to iSyTE are expected to further aid the understanding of lens biology and the discovery of cataract-linked genes. Finally, although single-nucleus and single-cell transcriptome approaches are now being applied to lens studies, these do not retain spatial information, because they necessitate dissociation of individual cells or nuclei. The present study, based on meta-analysis of transcriptome data on isolated lens epithelium and fibers retain aspects of spatial information on the cell-types and therefore provide a good reference dataset for comparative analysis with scRNA-seq and snRNA-seq approaches. 

Acknowledgments

Supported by National Institutes of Health [R01 EY021505 and R01 EY036923] to S.L. Bioinformatics support from the University of Delaware was made possible through funding from the State of Delaware and National Institutes of Health / National Institute of General Medical Sciences INBRE Program Grant [P20 GM103446]. M.D. and L.P. were supported by a grant from Association Retina France, CNRS PICS. 

Disclosure: M. Duot, None; S.Y. Coomson, None; S.K. Shrestha, None; M.V.M.K. Nagulla, None; Y. Audic, None; R.A. Barve, None; H. Huang, None; C. Gautier-Courteille, None; L. Paillard, None; S.A. Lachke, None