Single-Cell RNA Sequencing of Rabbit Sclera at Different Developmental Stages: Unveiling Scleral Cells Atlas and the Heterogeneity of Fibroblasts

Abstract

Purpose: This study aims to construct a single-cell transcriptomic atlas of the developing rabbit sclera to elucidate fibroblast heterogeneity, differentiation trajectories, matrisome expression patterns, and intercellular communication, while revealing conserved molecular features of scleral cells through cross-species analysis.

Methods: Single-cell RNA sequencing (scRNA-seq) was performed on scleral tissues from New Zealand rabbits at embryonic day 25 (E25) and postnatal days 7 (P7), 21 (P21), and 180 (P180). Libraries were prepared using the DNBelab C Series Kit and sequenced on the BGISEQ-2000 platform. Sequencing reads were aligned to the OryCun2.0 genome using STAR, and unique molecular identifier (UMI) count matrices were generated with PISA. Data preprocessing was conducted using Seurat. Fibroblast lineage differentiation was analyzed via VIA, intercellular communication via CellChat, matrisome expression patterns via AUCell, and cross-species analyses via CACIMAR and hdWGCNA.

Results: We identified 7 major cell types and 15 subpopulations, with fibroblasts dominating the cellular landscape. Distinct fibroblast subtypes exhibited varied expression profiles and functions: KERAlow SPARCL1⁺ fibroblasts showed stem/progenitor-like features, while KERAhigh myocilin (MYOC)⁺ fibroblasts displayed senescence-associated phenotypes. Matrisome analysis revealed dynamic alterations in collagen and extracellular matrix (ECM)-related genes, and intercellular communication analysis highlighted complex signaling networks, particularly the MDK/PTN pathway. Cross-species comparisons demonstrated high conservation of fibroblasts between rabbit and human sclera, identifying four conserved co-expression modules.

Conclusions: This study presents the first single-cell atlas of rabbit scleral development, unveiling fibroblast heterogeneity, ECM remodeling mechanisms, and cross-species conserved features. These findings enhance our understanding of scleral biology and provide valuable insights for future research on ocular development and associated diseases, including myopia.

Figure 1.

Comprehensive single-cell landscape atlas of the sclera of New Zealand rabbits at different developmental stages. (A) Schematic of single‐cell sequencing in rabbit sclera. (B, C) The 133,673 scleral single cells were visualized using uniform manifold approximation and projection (UMAP), colored according to Seurat unsupervised clustering and identified cell types. (D) Bar graph showing the number of each cell type. (E) Dot plot identifies marker genes for each cell type, and subpopulations of cells with similar transcriptional states are clustered. (F) Scale diagram showing the distribution of each cell type at different developmental stages and among samples, left: are the different developmental stages, and right: are the different samples. (G) UMAP visualizing the marker genes identified by FindALLMarkers for each cell type of the sclera. SWCs, Schwann cells; ECs, endothelial cells; IMCs, immune cells; cycling_FBs, cycling_fibroblasts.

Figure 2.

Subcluster analysis of rabbit scleral fibroblasts. (A, B) UMAP plots of fibroblasts colored by cell subtype and different developmental stages. (C) Scale plots showing the distribution of each fibroblast subtype at different developmental stages and among samples, left: are the different developmental stages, right: are the different samples. (D) UMAP visualization of marker genes used to distinguish fibroblast subtypes. (E) Heatmap showing differential genes across fibroblast subtypes. (F) Box line plot showing ROGUE scores for individual cell types. (G) Heatmap showing the Spearman correlation between each fibroblast subtype and mural cells. (H) Dot plot showing the GO enrichment results for each fibroblast subtype. (I) Violin-box line plots showing the cellular senescence GSVA scores for each fibroblast subtype. * Indicates significance, *: P < 0.05, **: P < 0.01, ***: P < 0.001.

Figure 3.

Pseudotime analysis of fibroblasts. (A) Cytotrace stemness scoring of fibroblasts at stage E25. (B) VIA-inferred PARC clustering plot with arrows indicating direction of differentiation, colored by fibroblast subtype, developmental stage, and proposed time. (C) Vector flow map embedded according to VIA results, colored by fibroblast subtype. (D) Atlas View shows more complex differentiation trajectories of fibroblasts. (E) Proposed temporal heatmap of differential transcription factors at different developmental stages. (F) Expression of HIF1A, NFE2L1, and NFE2L2 with the proposed chronology, with increased expression at the end of differentiation. (G) The Venn diagram shows the number of specific transcription factors of KERA^low FB at each developmental stage. (H) The Venn diagram shows the number of specific transcription factors of KERA^high FB at each developmental stage.

Figure 4.

Matrisome expression patterns in fibroblasts. (A, B) UMAP density plots showing the activity of each matrisome gene set in fibroblasts. (C, D) Violin box line plots showing the AUCell scores of each matrisome gene set for the two fibroblast subtypes at each developmental stage. (E, F) Violin box line plots showing the AUCell scores of the matrisome gene sets for each developmental stage of fibroblasts. (G) Violin-box line plots showing the expression of major collagen genes at different developmental stages of fibroblasts. * Indicates significance, *: P < 0.05, **: P < 0.01, ***: P < 0.001; NS, no significance.

Figure 5.

Differential matrisome genes of KERA^high FB and KERA^low FB at each developmental stage. (A–-D) Volcano plots showing the differential matrisome genes of KERA^high FB and KERA^low FB at each developmental stage, with the 20 most significant genes labeled and ECM-AFFILIATED genes common to each developmental stage highlighted.

Figure 6.

Sclera cellular communication network. (A) Number of cellular interactions and strength of interactions at different developmental stages. Left: The number of interactions, and right: the strength of interactions. (B) Scatter plots showing the signal strengths sent and received by each cell type at different developmental stages. The abbreviations of cell types are the same as in Figure 1. (C) Bar graph showing the information flow of signaling pathways at different developmental stages. (D) Patterns of signals sent by each cell type at each developmental stage. (E) Signal-receiving patterns of each cell type at each developmental stage. (F) Heat map showing the expression pattern of MDK signaling at each developmental stage. (G) Ligand-receptor pair communication probability of KERA^high FB with the MDK/PTN signaling pathway for each cell type. (H) Ligand-receptor pair communication probability of KERA^low FB with the MDK/PTN signaling pathway for each cell type.

Figure 7.

Cross-species analysis of scleral cells of humans and rabbits. (A) Heatmap showing the degree of conservation between human and rabbit scleral cells. (B) Heatmaps showing the conserved genes of human and rabbit scleral cells. (C) Bar plots showing the top 10 characteristic genes of each module. (D) Bar plots showing the enrichment results of genes in each module, with the top 10 results of enrichment displayed. (E) Venn diagram showing the number of myopia genes expressed in humans and rabbits. (F) Violin-box line plots showing the myopia gene set score of human and rabbit scleral cells. (G) Venn diagram showing the highly expressed myopia genes in human and rabbit scleral fibroblasts. * Indicates significance, ***: P < 0.001.

Figure 8.

PinpoRNA in situ hybridization for the detection of KERA and SPARCL1 mRNA expression in the sclera at different developmental stages.