Transcriptomic landscape of quiescent and proliferating human corneal stromal fibroblasts

Abstract

This study analyzed the transcriptional changes in primary human corneal stromal fibroblasts (hCSFs) grown under quiescent (serum-free) and proliferating (serum-supplemented) culture conditions to identify genes, pathways, and protein‒protein interaction networks influencing corneal repair and regeneration. Primary hCSFs were isolated from donor human corneas and maintained in serum-free or serum-laden conditions. RNA was extracted from confluent cultures using Qiagen kit and subjected to RNA sequencing (RNAseq) analysis. Differential gene expression (DGE) and pathway enrichment analyses were conducted using DESeq2 and Gene Set Enrichment Analysis (GSEA), respectively. Protein‒protein interaction (PPI) networks were created exploiting the STRING database and analyzed with Cytoscape and the cytoHubba plugin. RNA-seq revealed 5,181 genes that were significantly differentially expressed/changed among the 18,812 annotated genes (p value ˂0.05). A cutoff value of a log2-fold change of ±1.5 or greater was used to identify 674 significantly upregulated and 771 downregulated genes between quiescent and proliferating hCSFs. Pathway enrichment analysis revealed significant changes in genes linked to cell cycle regulation, inflammatory, and oxidative stress response pathways, such as E2F Targets, G2M Checkpoint, and MYC Targets, TNFA signaling via NF-kB, and oxidative phosphorylation. Protein-protein interaction network analysis highlighted critical hub genes. The FGF22, CD34, ASPN, DPT, LUM, FGF10, PDGFRB, ECM2, DCN, VEGFD, OMD, OGN, ANGPT1, CDH5, and PRELP were upregulated, whereas genes linked to cell cycle regulation and mitotic progression, such as BUB1, TTK, KIF23, KIF11, BUB1B, DLGAP5, NUSAP1, CCNA2, CCNB1, BIRC5, CDK1, KIF20A, AURKB, KIF2C, and CDCA8, were downregulated. The RNA sequences and gene count files have been submitted to the Gene Expression Omnibus (accession # GSE260476). Our study provides a comprehensive information on the transcriptional and molecular changes in hCSFs under quiescent and proliferative conditions and highlights key pathways and hub genes.

Keywords: Cornea; Fibroblasts; Proliferative; Quiescent; RNA sequencing; Stroma.

Fig. 1.

Experimental design and workflow used to conduct transcriptomic and RNA-seq analyses of quiescent and proliferating human corneal stromal fibroblasts treated with or without 10% fetal bovine serum.

Fig. 2.

Funnel chart illustrating the results of the differential gene expression analysis. The funnel chart depicts all annotated genes (18,812), significantly changed genes (P-adj < 0.05) (5,181), and significantly altered genes with a log2(FC) cutoff of ±1.5 (1,391 genes; 674 upregulated genes with a log2(FC) ≥ 1.5 and 717 downregulated genes with a log2(FC) ≤ −1.5).

Fig. 3.

Volcano plot of differential gene expression in quiescent and proliferating human corneal stromal fibroblasts. The x-axis represents the log2-fold change in gene expression, and the y-axis represents the −log10 of the p value. Each point represents a gene, color-coded based on significance and magnitude of change: red dots for significant genes with log2 (FC) > 1.5, blue dots for significant genes with log2 (FC) < 1.5, green dots for nonsignificant genes with log2 (FC) > 1.5, and gray dots for nonsignificant genes with log2 (FC)< 1.5. The dashed lines indicate thresholds for statistical significance [p value < 0.05] and log2 (FC) = 1.5]. Notable differentially expressed genes are labeled.

Fig. 4.

Heatmap of gene expression and enriched pathways in quiescent and proliferating human corneal stromal fibroblasts. This heatmap displays the differential expression of genes between proliferating (P-hCSF-A1 to P-hCSF-A4) and quiescent (Q-hCSF-C1 to Q-hCSF-C4) human corneal stromal fibroblasts. The rows represent individual genes, while the columns represent different samples. Gene expression levels are color-coded, with blue indicating high expression and red indicating low expression, as shown by the color scale on the right. The left side of the heatmap displays the most enriched pathways identified through gene set enrichment analysis. Hierarchical clustering dendrograms above and to the left highlight distinct expression patterns and biological differences between Q-hCSF and P-hCSF.

Fig. 5.

Pathway enrichment analysis and heatmaps of enriched genes in the Q-hCSF and P-hCSF samples. (A) The bar chart shows the results of the pathway enrichment analysis performed via gene set enrichment analysis (GSEA), which ranks pathways based on their normalized enrichment score (NES), which indicates the degree to which genes in each pathway are overrepresented at the top or bottom of the ranked list of genes. The highest ranked pathways included E2F targets, the G2M checkpoint, MYC targets V1 and V2, TNFA signaling via NFKB, and Oxidative phosphorylation. (B-G) Heatmaps depict the expression profiles of genes specifically enriched in significant pathways: (B) E2F targets, (C) G2M checkpoint, (D) MYC targets V2, (E) MYC targets V1, (F) TNFA signaling via NFKB, and (G) Oxidative phosphorylation. The columns represent individual samples of Q-hCSF and P-hCSF, and the rows represent genes, with color scales indicating expression levels from high (blue) to low (red).

Fig. 5.

Pathway enrichment analysis and heatmaps of enriched genes in the Q-hCSF and P-hCSF samples. (A) The bar chart shows the results of the pathway enrichment analysis performed via gene set enrichment analysis (GSEA), which ranks pathways based on their normalized enrichment score (NES), which indicates the degree to which genes in each pathway are overrepresented at the top or bottom of the ranked list of genes. The highest ranked pathways included E2F targets, the G2M checkpoint, MYC targets V1 and V2, TNFA signaling via NFKB, and Oxidative phosphorylation. (B-G) Heatmaps depict the expression profiles of genes specifically enriched in significant pathways: (B) E2F targets, (C) G2M checkpoint, (D) MYC targets V2, (E) MYC targets V1, (F) TNFA signaling via NFKB, and (G) Oxidative phosphorylation. The columns represent individual samples of Q-hCSF and P-hCSF, and the rows represent genes, with color scales indicating expression levels from high (blue) to low (red).

Fig. 6.

Gene set enrichment analysis of key pathways in the Q-hCSF and P-hCSF samples. The enrichment plots illustrate the results of GSEA for six hallmark pathways that are significantly enriched in the comparison between the Q-hCSF and P-hCSF samples. Each panel displays the enrichment score curve, where the peak of the green line represents the enrichment score for the gene set. The plots also show the location of genes in the pathway within the ranked list of genes, with red indicating positive correlation (enrichment in Q-hCSF) and blue indicating negative correlation (enrichment in P-hCSF). The normalized enrichment score (NES) and the false discovery rate (FDR) q value are annotated in each plot.

Fig. 7.

Gene ontology analysis of DEGs and identification of biological processes and comodulated genes in Q-hCSF compared to P-hCSF. (A) Dot plot showing the enriched biological processes identified through Gene Ontology analysis of DEGs in the Q-hCSF group compared with the P-hCSF group. The x-axis represents the gene ratio, with the dot size indicating the count of genes involved in each process and the color representing the adjusted P value. (B) Cluster dendrogram displaying the hierarchical clustering of the enriched biological processes, with branches colored according to specific categories such as spindle checkpoint signaling and mitotic nuclear division. (C) Network diagram illustrating the comodulated genes and their associated biological processes, highlighting key genes. The node size represents the degree of connectivity, and the color indicates the log2-fold change. (D) Circle plot showing the interaction network of DEGs, with node size indicating the number of genes and edge thickness representing the strength of interactions among biological processes.

Fig. 7.

Gene ontology analysis of DEGs and identification of biological processes and comodulated genes in Q-hCSF compared to P-hCSF. (A) Dot plot showing the enriched biological processes identified through Gene Ontology analysis of DEGs in the Q-hCSF group compared with the P-hCSF group. The x-axis represents the gene ratio, with the dot size indicating the count of genes involved in each process and the color representing the adjusted P value. (B) Cluster dendrogram displaying the hierarchical clustering of the enriched biological processes, with branches colored according to specific categories such as spindle checkpoint signaling and mitotic nuclear division. (C) Network diagram illustrating the comodulated genes and their associated biological processes, highlighting key genes. The node size represents the degree of connectivity, and the color indicates the log2-fold change. (D) Circle plot showing the interaction network of DEGs, with node size indicating the number of genes and edge thickness representing the strength of interactions among biological processes.

Fig. 8.

Protein–protein interaction network of upregulated genes in quiescent vs. proliferating primary human corneal stromal fibroblasts.

Fig. 9.

Protein–protein interaction network of downregulated genes in primary quiescent vs proliferating-hCSF.