Identification and characterization of glycosylation-related CAFs subgroups in CRC using scRNA-seq data

In this study, we included 10 untreated CRC samples from GSE205506. First, the data underwent quality control to remove low-quality and abnormal cells. We utilized the “Harmony” R package to eliminate batch effects between samples (Fig. 2A). Subsequently, we performed dimensionality reduction and clustering on the data, followed by cell annotation based on known marker genes (Li et al. 2023; Zhou et al. 2024; Chu et al. 2024) (Fig. 2B & S1A). We also used the “inferCNV” R package with B cells as the reference to identify malignant cells (Fig. S1B). Ultimately, we identified a total of 50,015 cells, including 26,166 malignant cells, 2,825 epithelial cells, 6,299 B cells, 1,529 cancer-associated fibroblasts (CAFs), 8,566 T/I/NK cells, and 4 proliferating cells (Fig. 2C & D).

Fig. 2

Analysis of scRNA-seq data and identification of cell types. A t-SNE plot demonstrating the clustering results subsequent to batch effect correction, with distinct color-coding for individual samples. B Dot plot depicting the expression levels of marker genes across the identified cell types. C t-SNE plot visualizing the clustering of the various cell types. D Bar plot presenting the distribution of the identified cell types within each individual sample. t-SNE, T-distributed Stochastic Neighbor Embedding

To further compare glycosylation level differences among various cell types in CRC, we conducted an in-depth analysis of scRNA-seq data. Initially, a univariate Cox regression analysis was performed on glycosylation-related genes from previous studies (Table S2) (Qiu et al. 2023; Wu et al. 2023; Liang et al. 2022), which were then intersected with differentially expressed genes from the TCGA-CRC dataset. As a result, 15 glycosylation-related genes associated with CRC prognosis were identified. Our analysis indicated that these genes are predominantly expressed in malignant cells and CAFs (Fig. S2A). Using these 15 genes, we constructed a CRC-specific glycosylation-related gene set. Multiple scoring algorithms, including AUCell, UCell, singcore, ssGSEA, and AddModuleScore, were applied to evaluate glycosylation levels across different cell types in CRC. The results demonstrated that glycosylation levels were highest in CAFs (Fig. 3A).

Fig. 3

Identification and characterization of glycosylation-related CAFs subgroups in CRC using scRNA-seq data. A Dot plot illustrating the glycosylation gene set scores obtained from various scoring methods across different cell types. B UMAP plot showing the clustering of CAFs based on the NMF algorithm, highlighting distinct subgroups with varying glycosylation profiles. C Volcano plot depicting differentially expressed genes between CAFs subgroups with distinct glycosylation profiles. D Bubble plot presenting results from the GO enrichment analysis for CAFs subgroups characterized by glycosylation profiles. E Pseudotime trajectory analysis of CAFs subgroups, displaying variations in glycosylation subgroups (left) and pseudotime coloring of these trajectories (right). F Circular plot showing the number of cell–cell interactions among the subgroups. G Heatmap illustrating the contribution of outgoing and incoming pathways in cell communication. H Heatmap demonstrating metabolic state differences among various CAFs subgroups. UMAP, Uniform Manifold Approximation and Projection. GO, Gene Ontology

Subsequently, we explored the differences in glycosylation profiles within CAFs in CRC. The NMF algorithm was employed for dimensionality reduction and clustering of CAFs, leading to the identification of 16 distinct CAFs clusters. Differential expression analyses of these clusters revealed genes uniquely expressed in each cluster. Based on the biological functions of these genes in glycosylation, we merged the clusters, leading to a refined classification of CAFs according to their glycosylation profiles. Ultimately, 6 glycosylation-related CAFs subgroups were identified: Glycosaminoglycan-C1, Glycosphingolipid-C2, N-glycosylation-C3, O-GalNAc-C4, Other-C5, and None-C6 (Fig. 3B & C).

We further analyzed these 6 glycosylation-related CAFs subgroups to investigate their intrinsic heterogeneity. GO Enrichment analysis (Fig. 3D & S2B) revealed that Glycosaminoglycan-C1 is associated with the glycosaminoglycan biosynthetic process, while N-glycosylation-C3 is linked to protein N-linked glycosylation. These results further validate the accuracy of the NMF algorithm in identifying glycosylation-related CAFs subgroups. Notably, Glycosphingolipid-C2 was found to be related to antigen processing and presentation, as well as the positive regulation of T cell-mediated immunity.

The different glycosylation-related CAFs subgroups also exhibit significant heterogeneity in their developmental trajectories (Fig. 3E). Glycosaminoglycan-C1 and Glycosphingolipid-C2 are distributed across various developmental stages, whereas O-GalNAc-C4 is primarily involved in intermediate and late developmental stages.

CellChat (Jin et al. 2024) analysis further revealed the interactions among various cell types and the signaling patterns of incoming and outgoing signals in different CAFs subgroups (Fig. 3F & G). In terms of outgoing signaling pathways, O-GalNAc-C4 serves as the primary secretory subgroup among CAFs, transmitting signals through VEGF, PERIOSTIN, ANGPTL, FGF, CALCR, and COMPLEMENT signaling pathways. When CAFs subgroups act as the target cells in incoming signaling pathways, they are primarily driven by PDGF, PERIOSTIN, and MK signaling pathways. We further explored the potential ligand-receptor communication between CAFs subgroups and other cell types. The analysis indicated that MAMPT-(ITGA5 + ITGB1) signaling is particularly active between CAFs and other cells (Fig. S2C).

Additionally, we compared the relative activity of various metabolic pathways across the fibroblast subgroups using “scMetabolism” R package (Wu et al. 2022) (Fig. 3H). Glycosaminoglycan-C1 demonstrated heightened activity in sulfur metabolism, while Glycosphingolipid-C2 was notably active in steroid hormone biosynthesis and the pentose phosphate pathway. N-glycosylation-C3 showed increased activity in linoleic acid metabolism, and O-GalNAc-C4 exhibited elevated activity in fatty acid metabolism-related processes. This highlights the diverse metabolic landscape of glycosylation-related fibroblast subgroups, offering insights into their distinct functional roles and potential metabolic vulnerabilities, which could serve as therapeutic targets.

In conclusion, the glycosylation-related CAFs subgroups display substantial heterogeneity in terms of biological functions, intercellular communication, metabolic profiles, and developmental trajectories. These differences likely reflect their distinct roles within the TME, further emphasizing their potential as therapeutic targets.

A prognostic model based on glycosylation-related genes effectively predicts the outcomes in CRC patients

Given the intrinsic heterogeneity exhibited among different CAFs subgroups, it is important to investigate whether significant prognostic differences exist in CRC patients with different CAFs subgroups. Specific gene signatures were constructed based on the differential genes of each subgroup, and their prognostic significance was compared in the training and validation sets (Fig. S3A). It was observed that patients with higher scores in the C2, C4, and C5 subgroups in both the training and validation sets had worse prognoses (Fig. 4A-F & S3B-D). Next, we aimed to more precisely predict the prognosis of CRC patients by analyzing the differential genes associated with these three CAFs subgroups. Therefore, the differential genes from these subgroups were intersected with the differentially expressed genes in tumor and normal samples from TCGA-CRC, and univariate Cox regression analysis was performed, identifying 19 prognostically relevant differential genes. Subsequently, 101 predictive models were developed using 10 different machine learning algorithms, including CoxBoost, Enet, GBM, Lasso, plsRcox, Ridge, RSF, StepCox, SuperPC, and survival-SVM (Fig. 4G). Ten-fold cross-validation was employed, and the robustness of these models was evaluated based on the average C-index in the training and validation sets and the number of genes included in the final models (Fig. S3E, Table S3).

Fig. 4

A prognostic model based on glycosylation-related genes effectively predicts the outcomes in CRC patients. A-F Kaplan–Meier plots illustrating prognostic differences among differential genes identified within subgroups C2, C4, and C5 across the TCGA (overall survival) and Meta-GEO cohorts (relapse-free survival). G Heatmap presenting the performance of 101 machine learning algorithms employed to construct prognostic risk models, ranked by the average C-index calculated across the TCGA and Meta-GEO cohorts. H–K Analysis of prognostic differences between high-risk and low-risk patient groups in the TCGA and Meta-GEO cohorts, accompanied by ROC curves for 1, 2, and 3-year prognostic projections

Notably, the StepCox[both] and RSF algorithms combination demonstrated the second-highest average C-index (0.757) and included only 5 genes (CCDC85B, MYL9, PLOD3, PSMA7, and SLC22A17), compared to the top-ranked RSF model. Therefore, it was selected as the best model. To comprehensively evaluate the robustness of the model, the risk scores for each patient in the training and validation sets were calculated and divided into low-risk and high-risk groups based on the median risk score. Kaplan–Meier analysis revealed that patients with lower risk scores had better overall survival (OS) and relapse-free survival (RFS) in the training set, and better RFS in the validation set. Consistently, the area under the ROC curve (AUCs) for 1-year, 2-year, and 3-year OS in the training set were 0.98, 0.99, and 0.99, respectively; the AUCs for 1-year, 2-year, and 3-year RFS in the training set were 0.78, 0.73, and 0.80 (Figure. S3F&G); and the AUCs for 1-year, 2-year, and 3-year RFS in the validation set were 0.61, 0.57, and 0.57, respectively (Fig. 4H–K). These results indicate that the model can accurately predict the prognosis of CRC patients. Overall, these findings suggest that the glycosylation-related gene model exhibits stable and robust performance in a large cohort of CRC patients.

The signature scores of glycosylation-related genes reshape the TME in CRC patients

Recent studies have suggested that aberrant glycosylation modifications may alter the interactions between tumor cells and immune or stromal cells, which, in turn, impacts intercellular communication and TME homeostasis. Therefore, a comparison of TME differences was further conducted between patients in the high-risk and low-risk groups. Multiple immune infiltration algorithms were employed to assess the infiltration levels of various immune cell types (Fig. 5A & S4A). CIBERSORT analysis indicated increased infiltration of M2 macrophages in high-risk patients, whereas low-risk patients demonstrated higher levels of T_cells_CD4_memory_activated and plasma cell infiltration (Fig. 5B). ESTIMATE analysis indicated that high-risk patients had significantly higher ESTIMATE, immune, and stromal scores compared to low-risk patients (p < 0.001), whereas low-risk patients demonstrated a higher tumor purity score (Fig. 5C-F). ssGSEA analysis revealed higher infiltration levels of gamma delta T cells, regulatory T cells, and MDSCs in high-risk patients, suggesting an immunosuppressive TME (Fig. 5G). Furthermore, the ssGSEA scores of various immune pathways were calculated for each patient, enabling a comparison of immune status between the different risk groups. The results demonstrated that high-risk patients were more active across various immune pathways (Fig. 5H). Expression levels of immune checkpoint molecules were also compared between the high- and low-risk groups, revealing higher expression levels in high-risk patients (Fig. 5I & S4B), suggesting potential immune evasion in these patients when undergoing immunotherapy.

Fig. 5

The signature scores of glycosylation-related genes reshape the TME in CRC patients. A Evaluation of immune infiltration differences between patients of varying risk levels using four distinct algorithms (TCGA cohort). B Violin plot depicting the proportions of 22 immune cell types infiltrating high- and low-risk patients, as determined by CIBERSORT. C-F Violin plots illustrating differences in stromal score, immune score, ESTIMATE score, and tumor purity among risk groups. G Boxplot showing ssGSEA scores for 22 immune cell types in high- and low-risk patients. H Heatmap displaying ssGSEA scores of immune-related pathways in high- and low-risk patients. I Boxplot demonstrating the expression levels of immune checkpoint molecules in high- and low-risk patient groups(*p < 0.05; **p < 0.01; ***p < 0.001)

Predictive value of glycosylation-related gene signature for immunotherapy

Next, the predictive efficacy of the model on patient responses to immunotherapy was investigated. High-risk patients had higher TIDE scores (Fig. 6A), suggesting that this group is more likely to experience immune evasion, leading to poorer immunotherapy outcomes. The effectiveness of the risk model was further validated using immunotherapy datasets (GSE78220 and IMvigor210), where patients in the PD/SD group had higher risk scores (Fig. 6B, Table S4). The bar chart clearly displayed the distribution of treatment responses between high- and low-risk patients, with a significantly lower proportion of CR/PR in high-risk patients compared to low-risk ones (Fig. 6D & H). Subclass Mapping (Submap) analysis suggested greater immunotherapy sensitivity in low-risk patients. It was also observed that high-risk patients had poorer prognoses (Fig. 6C & E–G). Notably, in the IMvigor210 cohort, high-risk patients had poorer OS in early-stage (stage I and II) cases (p = 0.04) rather than in late-stage (stage III and IV) cases, highlighting the model's strong predictive potential in early-stage patients. This comprehensive analysis emphasizes the capability of the prognostic model to predict immunotherapy outcomes in CRC patients.

Fig. 6

Predictive value of glycosylation-related gene signature for immunotherapy. A Correlation analysis between risk scores and TIDE scores in the TCGA cohort. B Correlation between risk scores and immunotherapy response in the GSE78220 cohort. C Kaplan–Meier curve depicting overall survival of high- and low-risk patients in the GSE78220 cohort. D Bar chart illustrating the proportion of immunotherapy responses in high- and low-risk patients within the GSE78220 cohort. E–G Kaplan–Meier curves showing overall survival for all patients (E), early-stage patients (F), and late-stage patients (G) in the IMvigor210 cohort. H Bar chart representing the proportion of immunotherapy responses in high- and low-risk patients in the IMvigor210 cohort. I Submap analysis assessing the immunotherapy response of high- and low-risk patients in the TCGA cohort. (*p < 0.05; **p < 0.01; ***p < 0.001)

A nomogram integrating glycosylation-related gene signature scores and clinical parameters for precise prediction

In seeking other clinical applicability of the glycosylation-related gene signature, univariate and multivariate Cox regression analyses were conducted to further explore potential factors influencing the prognosis of CRC patients (Fig. 7A & B). The results indicated that the risk score is an independent risk factor for CRC. Furthermore, the correlation between the risk score and various clinicopathological features was compared (Fig. 7C). In the low-risk group, there were more early-stage (Stage I and Stage II) patients, whereas the high-risk group demonstrated a higher proportion of late-stage (Stage III and Stage IV) patients (p < 0.001). Additionally, the low-risk group had a higher proportion of patients in the early T stages (T1 and T2) (p < 0.001). There were also significant differences between the high- and low-risk groups in terms of lymph node metastasis and distant metastasis (p < 0.001), with the high-risk group exhibiting a higher proportion of patients with lymph node and distant metastases. However, MSI status and gender were not significantly correlated with the risk score. Next, a predictive nomogram was developed that combines the risk score with key clinicopathological parameters (N stage) to predict the 1-year, 3-year, and 5-year OS (Fig. 7D). Calibration curves confirmed that the nomogram's predicted 1-year, 3-year, and 5-year survival rates were highly consistent with the actual survival outcomes (Fig. 7E). Decision curve analysis (DCA) showed that the nomogram and risk score provided the greatest standardized net benefit at most high-risk thresholds (Fig. 7F). Additionally, compared to other clinical indicators, the nomogram and risk score maintained the highest concordance index over a 10-year period, demonstrating their superiority and stability in long-term prediction (Fig. 7G).

Fig. 7

A nomogram integrating glycosylation-related gene signature scores and clinical parameters for precise prediction. A Univariate independent prognostic analysis in the TCGA cohort. B Multivariate independent prognostic analysis in the TCGA cohort. C Association between various clinicopathological characteristics and risk scores. D Nomogram integrating clinical features, specifically N staging, with risk scores. E Calibration plot assessing the agreement between actual OS rates and predicted survival rates, with the 45-degree line representing an ideal prediction. (F) DCA demonstrating that the nomogram offers significantly greater clinical benefit compared to other clinical features. (G) AUC curve showcasing the predictive performance of various clinical characteristics, nomogram scores, and risk scores

Prediction of biological mechanisms associated with glycosylation-related gene signature

To investigate the potential biological mechanisms associated with the prognosis of patients in the high-risk and low-risk groups, differences in somatic mutation patterns within the TCGA cohort and functional pathway enrichment analysis were examined. The top 5 differential mutations were compared in patients with different risk levels (Fig. 8A). Notably, patients in the high-risk group exhibited a higher incidence of insertion and deletion mutations in APC and deletion mutations in TP53, whereas patients in the low-risk group showed more alternative splicing in TNN. Moreover, among the top 20 genes by mutation frequency in high-risk and low-risk patients, it was found that genes from the Ryanodine Receptor family (RYR1, RYR2 and RYR3) showed a higher mutation frequency (Fig. 8B & C). Studies indicate that RYR1 and RYR2 mutations may enhance the release of intracellular calcium ions, thereby affecting cancer cell growth and drug resistance. Consequently, high-risk patients may not respond well to therapies that modulate calcium ion signaling pathways. Significantly different somatic mutation gene interaction patterns were discovered between high-risk and low-risk patients (Fig. 8D & E). In the high-risk group, RYR1 and TP53, as well as RYR3 and KRAS, were mutually exclusive genes, a pattern not observed in the low-risk group. Additionally, there were significant differences in co-occurrence gene interaction patterns between the two groups. A series of enrichment analyses were also conducted for the high-risk and low-risk groups to explore the underlying biological mechanism differences that might exist. The results showed that pathways such as HALLMARK_APICAL_JUNCTION, HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_INFLAMMATORY_RESPONSE, and HALLMARK_KRAS_SIGNALING_DN are enriched in the high-risk group (Fig. 8F). On the other hand, pathways such as HALLMARK_G2M_CHECKPOINT, HALLMARK_SPERMATOGENESIS, HALLMARK_E2F_TARGETS, and HALLMARK_MYC_TARGETS_V1 are enriched in the low-risk group (Fig. 8G). Additionally, a correlation analysis was conducted between the risk score and the GSVA scores of each pathway, revealing that the risk score is positively correlated with the HALLMARK_KRAS_SIGNALING_DN pathway score (Fig. 8H).

Fig. 8

Prediction of Biological Mechanisms associated with glycosylation-related gene signature. A Waterfall plot illustrating the top 5 somatic mutation differences between high-risk and low-risk patients. B-C Waterfall plots showing the top 20 somatic mutations in high-risk (B) versus low-risk (C) patients. D-E Heatmaps depicting co-mutations and mutually exclusive mutations among high-risk (D) and low-risk (E) patients. F GSEA analysis indicating enriched pathways among high-risk patients. G GSVA analysis revealing differences in enriched pathways between high- and low-risk patients. H Heatmap illustrating the correlation between risk scores and various enriched pathways

Exploring potential therapeutic options for the high-risk patients

We further evaluated the response of CRC patients in the high- and low-risk groups to candidate drugs using the “pRRophetic” R package. We assessed each patient’s sensitivity scores to different drugs and compared their correlation with the risk scores (Table S5). A total of 103 compounds were identified that showed significant differences between the high- and low-risk groups. Erlotinib, VX-680, Lapatinib, A-443654, and Sorafenib were associated with higher sensitivity scores in high-risk patients, suggesting potential resistance to these compounds (Fig. 9A-E). Conversely, increased sensitivity to PLX4720, SL-0101–1, WO2009093972, CCT007093, and AMG-706 was observed in these patients (Fig. 9F-J). Among these, PLX4720 is identified as a selective BRAF^V600E inhibitor, suggesting that high-risk patients may benefit from combination therapy with BRAF inhibitors to enhance their prognosis.

Fig. 9

Exploring potential therapeutic options for the high-risk patients. A-E Analysis of differences in drug sensitivity for the top 5 drugs associated with resistance in high-risk patients, comparing high and low-risk groups (left panels), and their correlation with risk scores (right panels). F-J Assessment of differences in drug sensitivity for the top 5 drugs associated with sensitivity in high-risk patients. (*p < 0.05; **p < 0.01; ***p < 0.001)

Validating the expression and function of CCDC85B in the glycosylation-related gene signature

To advance understanding of the biological roles of the five genes highlighted in the glycosylation-related gene signature and elucidate their molecular mechanisms, a comprehensive investigation was undertaken. Initially, the relationships between each gene and the infiltration levels of various immune cells were assessed (Fig. S5A). These findings indicated that MYL9 showed a significant positive correlation with fibroblast infiltration levels and StromalScore, suggesting its potential role in stromal remodeling within the TME. CCDC85B was negatively correlated with the infiltration levels of plasmacytoid dendritic cells (pDCs) and CD4 + effector memory T cells, indicating its possible involvement in regulating the immunosuppressive microenvironment of CRC. Additionally, SLC22A17 exhibited significant positive correlations with the infiltration levels of multiple immune cell types. The ssGSEA analysis demonstrated that CCDC85B was positively correlated with KRAS_SIGNALING_DN, and PLOD3 with WNT_BETA_CATENIN_SIGNALING (Fig. S5B). To identify the most functionally central genes among the five-gene signature, we used the “GOSemSim” R package to calculate the semantic similarity scores for GO biological processes among the five genes. Genes with higher average semantic similarity scores were considered to play more crucial roles, as they are functionally more connected with other genes in the signature. This analysis revealed that PSMA7 and CCDC85B exhibited the highest scores, suggesting their central roles in glycosylation-related pathways in CRC (Fig. 10A). Further analyses showed that only CCDC85B displayed significantly higher expression in CRC tissues and correlated with worse patient prognosis, leading us to select CCDC85B for functional validation (Fig. 10B & C, S5C). Subsequently, CCDC85B was further investigated regarding its biological functions in CRC. The results demonstrated that CCDC85B expression was significantly higher in CRC cell lines (SW620, SW1116, and DLD1) compared to normal colonic epithelial cells (NCM460) (Fig. 10D & E). Immunohistochemical (IHC) analysis of surgically resected CRC tissues from FUSCC confirmed elevated unregulated CCDC85B expression in tumor tissues compared to adjacent normal issues (Fig. 10F & G). Furthermore, CCK8 and colony formation assays showed that CCDC85B knockdown inhibited CRC cell proliferation (Fig. 10H & I). The transwell assay demonstrated that CCDC85B knockdown reduced invasion capacity in CRC cells (Fig. 10J & K). Given that CCDC85B is a differential gene among various glycosylation CAFs subtypes, it was hypothesized that it might influence glycosylation levels in tumor cells. Western blot analysis showed a significant decrease in O-glycosylation levels in CRC cells with CCDC85B knockdown (Fig. 10L). Additionally, multiplex immunohistochemical (mIHC) staining of FUSCC CRC tissue microarrays revealed that higher CCDC85B expression was related to reduced CD8 + T cell infiltration (Fig. 10M). These results suggest that CCDC85B promotes the proliferation and migration of CRC cells in vitro and may serve as a potential therapeutic target for CRC.

Fig. 10

Validating the expression and function of CCDC85B in the glycosylation-related gene signature. A Comparative analysis of the importance of five genes (CCDC85B, MYL9, PLOD3, PSMA7, SLC22A17) using the “GOSemSim” R package. The box plot displays the relative importance of each gene in the model based on semantic similarity analysis. B Box plots showing the differential expression of CCDC85B in READ and COAD patients compared to normal tissue. C Kaplan–Meier survival analysis showing overall survival of CRC patients stratified by high and low CCDC85B expression levels. Patients with high CCDC85B expression show significantly lower survival rates. D Relative mRNA expression levels of CCDC85B in various CRC cell lines (HCT8, HCT15, HCT116, HT29, SW620, SW1116, RKO, DLD1,Caco2) and normal colon epithelial cell (NCM460). E Western blot analysis of CCDC85B protein levels in the CRC cell lines and normal colon epithelial cell. GAPDH was used as a loading control. F Immunohistochemical staining of CCDC85B in human colorectal cancer tissue (upper panel) and normal colon tissue (lower panel). G Violin plot showing the difference in immunohistochemical scores of CCDC85B between tumor and adjacent normal tissues in the FUSCC cohort. H Quantitative RT-PCR analysis showing the knockdown efficiency of CCDC85B in SW620 and HCT116 cells transfected with siRNA. Data are presented as mean ± SD. I Cell proliferation assay (OD values) of SW620 and HCT116 cells with and without CCDC85B knockdown. J Colony formation assay showing the colony-forming ability of SW620 and HCT116 cells following CCDC85B knockdown (Left). Quantification of cell colony numbers from colony formation assays in SW620 and HCT116 cells. Data are presented as mean ± SD (Right). K Transwell migration assay illustrating the migration ability of SW620 and HCT116 cells with and without CCDC85B knockdown (Left). Quantification of migrated cells from transwell assays for SW620 and HCT116. Data are presented as mean ± SD (Right). L Western blot analysis of O-glycosylation levels and CCDC85B protein in SW620 and HCT116 cells following CCDC85B knockdown. M Immunofluorescence staining of CCDC85B and CD8 in CRC tissues with low (upper panels) and high (lower panels) expression levels. DAPI was used for nuclear staining. Scale bars represent 50 μm. (*p < 0.05; **p < 0.01; ***p < 0.001)