MicroRNA Profiling of Highly Enriched Human Corneal Epithelial Stem Cells by Small RNA Sequencing

The objective of the study was to elucidate the microRNA (miRNA) profile of an enriched human corneal epithelial stem cell (CESC) population in comparison to differentiated central corneal epithelial cells (CCECs) by small RNA sequencing. The CESCs were enriched by differential enzymatic treatment to isolate the basal limbal epithelial cells followed by laser capture microdissection of cells with nucleus to cytoplasm ratio ≥0.7, from donor tissues. Small RNA sequencing was carried out using Illumina NextSeq. 500 platform and the validation of differentially expressed miRNAs by quantitative real-time PCR (qPCR) and locked nucleic acid miRNA in-situ hybridization (LNA-ISH). The sequencing identified 62 miRNAs in CESCs and 611 in CCECs. Six miRNAs: hsa-miR-21-5p, 3168, 143-3p, 10a-5p, 150-5p and 1910-5p were found to be significantly upregulated in enriched CESCs, which was further confirmed by qPCR and LNA-ISH. The expression of hsa-miR-143-3p was exclusive to clusters of limbal basal epithelial cells. The targets of the upregulated miRNAs were predicted to be associated with signaling pathways -Wnt, PI3K-AKT, MAPK and pathways that regulate pluripotency of stem cells, cell migration, growth and proliferation. Further studies are essential to elucidate their functional role in maintenance of stemness. The findings of the study also hypothesize the inherent potential of hsa-miR-143-3p to serve as a biomarker for identifying CESCs.

make use of a highly enriched CESCs in order to identify the profile of miRNA specific to adult tissue resident stem cells. Therefore, in this study we have made use of the strategy we have developed earlier 7,18 to enrich CESCs to the extent of 80% and this enriched population are known to express stem cell markers ΔNp63α and ABCG2 18 . It was thus possible to identify a specific profile of miRNA, significantly up-regulated in CESCs compared to differentiated CCECs 19 . Accordingly, the major findings of the study will serve as a platform to decipher the unresolved questions relating to the miRNA regulation of CESCs.

Methods
Samples. Human donor tissues were handled in accordance with the tenets of the Declaration of Helsinki and the study was approved by Institutional Ethics Committee Aravind Medical Research Foundation (RES2013038BAS). Informed consent was obtained for all donor eyes including the minors from the legally authorized representative -either the donor's parents or family through the Eye Banks of Aravind Eye Care System (Madurai, Coimbatore, Pondicherry and Tirunelveli). The donor globes not suitable for transplantation (procured within 8 hours of death and received within 24 hours for research) and limbal rims obtained after corneal transplantation (received within 10 days of storage in Cornisol media at 4 °C in sterile condition) were included in this study. Inclusion criteria were donor age less than 70 years (range: 11 to 70 years), non-vascularized, with no history of diabetes and ocular infection. The donor globes were observed under stereo binocular microscope and those with intact limbus with radial ridges of palisades of Vogt were used for the study.
Enrichment of corneal epithelial stem cells. The corneal region together with adjacent sclera was dissected out from the donor globes and the corneal endothelium was removed using a sterile cotton tip 20,21 . Central cornea was punched out using 8 mm trephine and the central corneal epithelial cells (CCECs) were isolated following the protocol of Arpitha et al. 20 . Briefly, the corneal epithelial sheet was separated by dispase II (2 mg/ml) treatment and the individual cells were obtained after 0.25% trypsin treatment. For enrichment of the CESCs, a two step-protocol was followed 18 . The limbal basal epithelial cells (LBECs) were isolated from limbal tissues by differential enzymatic treatment -0.25% trypsin followed by 2 mg/ml dispase II treatment 7 . The cytosmears of isolated LBECs on membrane slide were stained with Giemsa. The cells were focused with 20X magnification and the image was projected on the computer monitor connected to an automated microscope. The area of nucleus and cytoplasm of each cell and the nucleus to cytoplasmic (N/C) ratio was then calculated using PALM Robo Software 4.3 SP2. The corneal epithelial stem cells (CESCs) with N/C ratio >0.7 were selectively marked and cut using Positioning and Ablation Laser MicroBeam (PALM) and collected on the adhesive cap of the collection tube which positioned to overhang the area under focus. Catapulted cells collected in the adhesive cap can be checked by moving it directly over the objective lens at any point of time during the collection. The enriched CESCs in the adhesive cap of the collection tube were then incubated for 30 minutes at 42 °C in 15 µL of extraction buffer (ARCTURUS PicoPure RNA Isolation Kit) and collected to the bottom of the tube from the cap by a brief centrifugation. The collected cell lysate was stored at −80 °C until RNA isolation 18,22 . Small RNA sequencing and data analysis. The total RNA was extracted from pooled CESCs (2046 cells) and CCECs (8.55 × 10 5 cells) using ARCTURUS PicoPure RNA Isolation Kit from two donor tissues (Age: 60 and 62 years) for small RNA sequencing. Total RNA was used for small RNA library construction using TruSeq Small RNA Library Prep Kit (Illumina, California, USA) following manufacturer's protocol. Libraries were quantified using Qubit dsDNA HS Assay Kit (Invitrogen, California, USA) and validated in bioanalyzer using High Sensitivity DNA Kit (Agilent, California, USA). Based on the library QC report, the library generated was suitable for sequencing on Illumina. For CCECs, library distribution was in the range of 120bp-180bp and for CESCs the range was 120bp-160bp indicating the presence of small RNA inserts. The effective insert length was in the range of ~10bp-40 bp; with combined size of adapters flanking the library being ~120 bp. Cluster generation and sequencing was carried out using NextSeq. 500 High-Output v2 Kit (Illumina, Inc.) using 75 cycles chemistry at the Genotypic Technology Pvt Ltd (Bangalore, India).
First, the quality assessment of raw deep sequencing data was performed with the FastQC tool. Adapter and low-quality reads were discarded from raw sequencing data using our Perl script, allowing no mismatches for adapter identification. Further, the data were uploaded into Oasis 2.0 online software 23,24 for the miRNA analysis. Briefly, the data with the read size (15-32 nt) and low abundance reads (<5 reads) were discarded and remaining data were aligned to Homo sapiens hg19 genome reference. For differential expression (DE) analysis, NOISeq R software package were applied on the TMM normalized reads (R version 3.5.1). The variable expression of miR-NAs between two cell types was considered significant when the fold change was 2 or greater and the probability score was more than 0.9. The results were represented in the volcano plot constructed by R version 3.5.1 25 based on the NOISeq scores: the M-value and D-value.
Quantitative real time PCR. For confirmation of the sequencing data, quantitative real time PCR was carried out using (i) CESCs from 33 pairs of limbal rims (cells from 11 pairs were pooled as one sample to get around 3500 cells/pool, n = 3) and CCECs from three pairs of corneal button and (ii) LBECs and CCECs from limbal rims and corneal buttons respectively, were pooled to obtain a minimum of 1 × 10 6 cells each (pool 1: 7 pairs; pool 2: 7 pairs; pool 3: 6 pairs of tissues), but after ensuring that there was no bias of age. The RNA concentration was estimated using Qubit RNA HS Assay Kit (Invitrogen, California, USA) in Qubit2.0. fluorometer. RNA (15 ng) from CESCs (n = 3) and CCECs (n = 3) and 1 µg of RNA from LBECs (n = 3) and CCECs (n = 3) was reverse transcribed using miScript II RT Kit (Qiagen, Hilden, Germany) according to manufacturer's protocol. Due to the low quantity of RNA (15 ng) obtained, the reverse transcribed cDNA from CESCs was pre-amplified using miScript PreAMP PCR Kit (Qiagen, Hilden, Germany) for 12 cycles (Denaturation: 30 seconds at 94 °C and annealing/ extension: 3 minutes at 60 °C) preceded by initial activation: 15 minutes at 95 °C. Hence the same protocol was followed for CCECs for comparison. The qPCR amplification was performed for 40 cycles using 2X miScript SYBR Green PCR Master Mix, 10X miScript Universal Primer and miRNA specific miScript Primer Assays (Supplementary Table S1 for list of miRNA Primer Assays used and Table S2 for custom designed miRNA primer sequence). Experiment was done in triplicates and signals were normalized to small nucleolar RNA U6 (RNU6B) which was run in parallel as reference miRNA. Relative miRNA expression was calculated using comparative threshold cycle (Ct) method (2 -ΔΔCt ) and the result was represented as mean ± SEM. Statistical analysis of qPCR data was carried out using Mann-Whitney U test and miRNAs with fold change >2 and p value <0.05 were considered significant.
LnA in-situ hybridization. For validation of differential expression of miRNAs in its native environment, LNA in-situ hybridization was carried out following protocol of Obernosterer et al. 26 with some modifications. Briefly, The corneoscleral region from the donor eye ball was separated and fixed with 4% paraformaldehyde (Sigma-Aldrich, Missouri, United States) for 30 minutes at room temperature and then the tissue was incubated overnight at 4 °C in 30% sucrose solution prepared in 1X phosphate buffered saline (PBS) of pH 7.4 (Ambion, California, United States). The tissue was embedded in optimal cutting temperature (OCT) compound (Leica, Wetzlar, Germany), and cryosections (10 μm) of limbal epithelium and corneal epithelium were taken on SuperFrost Plus slides (Thermofisher Scientific, Massachusetts, United States). The cryosections were fixed with ice cold acetone for 10 minutes at room temperature and then hybridized with biotinylated LNA miRNA probes (Eurogentec, Liège, Belgium; Supplementary Table S3 for list of LNA miRNA probes used) at 50 °C overnight by placing coverslip over the sections to prevent drying. U6 small nucleolar RNA probe and LNA scrambled microRNA probe were used as positive and negative control respectively. After hybridization, the slides were soaked in prewarmed 5X SSC Buffer (Ambion, California, United States) to remove the coverslips. The sections were incubated at 60 °C for 1 hour in 0.2X SSC buffer followed by blocking with 10% fetal bovine serum (Invitrogen, California, USA) in 0.1 M Tris (pH 7.5) and 0.15 M sodium chloride solution at room temperature for 1 hour. To detect the hybridized probes, the sections were incubated with Streptavidin-Fluorescein Isothiocyanate (FITC) Conjugate (BD Biosciences, New Jersey, United States) in 1:1000 in 5% bovine serum albumin (Sigma-Aldrich, Missouri, United States) at 25 °C in dark for 2 hours. The sections were then washed twice with 1X PBS for 5 minutes and mounted using fluorescent mounting media with propidium iodide (Vector Laboratories, California, United States). The images were acquired using confocal laser scanning microscope (Leica SP8, Germany). miRNA Target prediction, Gene Ontology and Pathway Analysis. miRWalk (Version 3.0) and mirDIP (Version 4.1.1.6) were used to predict the target genes of miRNAs. Experimentally validated targets from miRWalk and which were targets with> 10 sources in mirDIP databases were filtered for the analysis. For pathway prediction, ClueGO package in Cytoscape 3.2.0 was used. KEGG pathway list was chosen with p-value ≤ 0.01. Further, the network of miRNAs, target genes and the associated pathways was constructed by Cytoscape 3.2.0 27 .

Statistical analysis.
Mann-Whitney U test was performed to determine the statistical significance between the two experimental groups and p < 0.05 was considered statistically significant. STATA 14.0 (Texas, USA) statistical software was used for analysis. Non-parametric test, Mann-Whitney U was chosen since the qPCR data followed non-gaussian distribution based on the Shapiro-Wilk normality test.

Results
Characteristics of small RNA sequencing and differentially expressed miRNAs between CESCs and CCECs. The entire quantity of total RNA isolated from enriched CESCs (pooled cells from two pairs of donor tissues) was used for library preparation as the concentration was lesser than the detection limit in the Qubit analysis. For CCECs, 10 ng of total RNA (RIN 2.8) was used. The sequencing generated more than 11 million reads in CESCs and CCECs. After quality assessment and refinement, the high quality reads were mapped to human hg 19 genome. A total of 38 miRNAs in CESCs and 301 miRNAs in CCECs were detected. Further, differential expression (DE) analysis using NOISeq after normalization were performed on all detected miR-NAs. A total of 127 DE miRNAs were identified in CESCs compared to CCECs with fold change >±2. Among them 29 miRNAs were upregulated and 64 miRNAs were downregulated. The fold change values of the miRNAs were summarized as volcano plot in Fig. 1. The top ten miRNAs that were highly expressed in CESCs (Fig. 2a) and CCECs (Fig. 2b) based on raw reads counts were represented graphically. Based on fold difference between CESCs and CCECs and raw reads counts, fourteen differentially expressed miRNAs and two miRNAs without significant difference in expression were selected for validation by qPCR expression analysis ( Table 1).
But only a few reports are available on association of miRNAs in the regulation of stemness. The major challenge of profiling the miRNAs in CESCs was the low stem cell content (3-5%) in the total limbal epithelium. Hence an attempt was made to profile the miRNAs specific to adult stem cells of the limbus using a small number of highly enriched CESCs as close as eighty percent 18 by small RNA sequencing platform. Since the yield of CESCs was minimum after enrichment, Arcturus picopure RNA isolation kit was used in this study for RNA isolation, to reduce the loss of short structured RNA with low GC during extraction from small number cells 42 .
Small RNA sequence data of CCECs was used as a quality control for the sequencing. Comparison of the data with published literature 29,43 using total cornea or corneal epithelium (by NGS or microarray) identified a similar profile of miRNAs, thereby confirming the validity of the data. However, only 38 miRNAs were identified in CESC population and this could possibly due to the following concerns (i) low input RNA, that could not be quantified before library preparation ii) usage of LCM for enrichment process and iii) difference in method of isolation of CESCs and CCECs. In spite of this drawback, miRNAs identified to be highly expressed in the CESCs were confirmed to be specific by in-situ hybridization and qPCR. www.nature.com/scientificreports www.nature.com/scientificreports/ Though the number of miRNAs detected in CESCs was only one tenth compared to CCECs, hsa-miR-3168, hsa-miR-21-5p, hsa-miR-143-3p, hsa-miR-150-5p, hsa-miR-1910-5p and hsa-miR-10a-5p had greater than seven folds increase in expression in the enriched CESCs indicating their specificity to stem cells. Among these six miRNAs, hsa-miR-3168 44 and hsa-miR-1910-5p 45 were reported as novel miRNAs in human embryonic stem cells. Understanding their functional role in regulating adult tissue resident stem cells will be of significance. The remaining four miRNAs were reported to have a role in regulating other stem cells. Hsa-miR-21-5p suppresses marrow derived endothelial progenitor cells proliferation by activating TGFβ signaling through down regulation of WW domain-containing protein 1 (WWP1) 46 . Similarly, miR-10a down regulates proliferation of human cardiomyocyte progenitor cells by targeting GATA6 47 . Hsa-miR-143-3p and hsa-miR-150-5p are known to be associated with self-renewal of mouse embryonic stem cells and human liver cancer stem cells respectively. MiR-143 promotes self-renewal by suppressing de novo methyltransferase gene DNMT3A 48 and miR-150 inhibits cell proliferation by targeting C-MYB 49 .
Validation of the sequencing data by qPCR analysis revealed that fourteen differentially expressed miRNAs to have similar expression pattern (Fig. 3). In addition, all the differentially expressed miRNAs had several magnitudes of fold difference in CESCs than LBECs when compared to CCECs independently (Fig. S2). In addition, the up regulation of hsa-miR-3168 as well as the down regulation of hsa-let7a-5p and hsa-miR-204-5p was identified only with the enriched CESCs but not with LBECs. Thus, these observations highlight the importance of enrichment of CESCs and confirming the specificity of the identified miRNAs.
LNA in-situ hybridization revealed the location of the miRNAs in corneal and limbal tissue. Among the miRNAs that were highly expressed in CESCs, hsa-miR-3168, hsa-miR-21-5p, hsa-miR-150-5p and hsa-miR-1910-5p had higher expression in limbal basal epithelial layer compared to that of the corneal basal epithelial layer. Further, a few clusters of small cells had strong positivity (probably stem cells) compared to that of the other limbal basal epithelial cells. Likewise, the expression of hsa-miR-143-3p was confined only to clusters in limbal basal epithelial layer, indicating its very high specificity to stem cells.
The predicted targets of these miRNAs that were highly expressed in CESCs were identified to be associated with pathways regulating pluripotency of stem cells and other pathways including MAPK signaling, Ras signaling, Hippo signaling, ErbB signaling, Wnt signaling and PI3K-Akt signaling. Reports are available on the significance of Wnt 50 , PI3K/Akt 51 , MAPK 52 signaling pathways in the regulation of stem cells. The direct targets of the miRNAs (highly expressed in CESCs) predicted to be involved in these signaling pathways are listed in Table 2. These targets are known to be associated with cell proliferation, differentiation, migration and apoptosis. Further studies are therefore essential to confirm the functional role of the identified miRNAs and their targets in the regulation of stemness in CESCs.
One of the limitations of this study is the use of donor tissues. The RNA expression could have been altered due to the associated ischemia 53 and the additional processing of the cells for enrichment of stem cells 54 . In spite of the low input of RNA for sequencing and low number of miRNAs identified in this study, confirmation of the high expression of miRNAs hsa-miR-3168, hsa-miR-21-5p, hsa-miR-143-3p, hsa-miR-150-5p, hsa-miR-1910-5p and hsa-miR-10a-5p in CESCs indicates their significant role in the regulation of stem cells. Transfection studies using hsa-miR-143-3p and hsa-miR-150-5p mimics in the primary limbal epithelial cells (n = 3) identified their increased colony forming potential, specifically their ability to form holoclones (based on colony morphology as defined by Barrandon and Green 55 ), thereby indicating a strong regulatory influence of these miRNAs on maintenance of stemness 56,57 . Thus, the data generated serves as a platform to study the miRNAs associated with the maintenance of stemness. Further studies are essential to elucidate the functional role of these miRNAs in regulating the stemness in CESCs and the associated signaling pathways.

miRNA Targets Functional role Reference
Hsa-miR-143-3p MAPK7 Inhibits adipocyte proliferation and enhance their differentiation 58 Inhibits proliferation, cell migration and invasion of breast cancer 59

KRAS
Suppresses colorectal cancer cell growth 60

IGF1R
Regulates cell proliferation and apoptosis 61 Regulation of myogenesis 62

CREB1
A novel Wnt effector that may enhance epithelial mesenchymal transition of colorectal cancer cells 67

CDK2
Inhibits proliferation and tumorigenicity via retarding G1/S phase transition 68