Single Cell Analysis Dissects the Effects of Vitamin D on Genetic Senescent Signatures Across Murine Tissues

1. Introduction

Vitamin D is a steroid hormone obtained by ingestion or photoreaction in the skin cells [1]. This hormone is necessary in mineral and bone homeostasis, although new insights are proving vital in our understanding of the non-canonical effects of the vitamin, for it seems to exert action in maladies like cancer and age-related diseases [2].

Remarkably, vitamin D deficiency is very prevalent worldwide, for approximately a billion of the world’s population is at risk of presenting inadequate serum vitamin D concentrations [3]. This has caused massive public health efforts to establish interventions, such as enhancing the employment of fortified food with vitamin D; although VD deficiency is still present in countries which have taken said measure [4]. Consequently, physicians and researchers now have associated many diseases with these VD deficits depending on the magnitude of patient deficiency.

An important functional axis in aging alterations is cellular senescence. This phenomenon refers to replicative cell halt and continuous growth under determined stress or physiological cues like oncogenic stress, deoxyribonucleic acid (DNA) damage, or early embryonic development [5]. This process represents an alternative cell fate to, for instance, apoptosis, and has been viewed as a culprit of aging, mainly in non-communicable diseases such as diabetes mellitus, hypertension, and cancer [6]. Current research regarding this issue focuses on eliminating or taking advantage of the characteristics of these cells for treating diseases. As such, vitamin D points to antagonizing senescent cells in its anti-apoptotic, energy-demanding, and pro-inflammatory profile while it promotes its anti-proliferative and prodifferentiative characteristics. Thus, vitamin D may be a fit candidate for cellular senescence-focused therapeutics, also known as senotherapeutics [7].

Cellular senescence has proven to be a challenging cellular state given its heterogeneity across cell types, and trigger mechanisms, so there is no definite biomarker nor genetic signature that certainly predicts if a cell is properly senescent [8]. Furthermore, vitamin D shows pleiotropism between tissues, and even cell types, making it difficult to target its non-canonical effects [9]. To tackle these complex interactions, a more comprehensive approach is needed to fully elucidate vitamin D role within senescence.

Along these lines, single cell approaches may be an exhaustive tool to thoroughly integrate the big picture of biological phenomena. In this particular case, single cell RNA sequencing can be a valuable experimental tool to suggest cellular function [10]. This is to characterize cell populations, cell states in interventions, key biological pathways and even differentiation pathways and cellular communication [11].

Here we explore vitamin D and cellular senescence phenomena from mouse skin, prostate and bone tissue from a single cell transcriptomic perspective. This with the purpose to unveil how this hormone affects the aging process in a non-canonical view which may, in turn, lead to future research lines to validate and suggest interventions against aging related diseases harnessing the effects of vitamin D.

2. Materials and Methods

2.1. Data Acquisition

A research in literature was made with the keywords vitamin D and single cell RNA seq in PUBMED and GoogleScholar. Six datasets were originally recovered. Four datasets were selected for further data exploration [12,13,14,15]. Nevertheless, further analysis was performed only with the databases of mouse bone, skin and prostate [12,13,14] given quality control results. It is worth noticing that the bone tissue datasets are from 8 healthy mice, the prostate dataset is from a PTEN- mouse and the skin tissue datasets were collected with UV-irradiated mice.

[12]

2.2. Data structure

Raw prostate dataset included: a vehicle count matrix, with its features and barcodes with 5163 cells and 17744 genes; and a VD count matrix, features and barcodes with 7912 cells and 17934 genes (GSE164858). Raw skin dataset included: an annotated vehicle count matrix with 17666 genes across 9081 cells and an annotated VD count matrix with 17257 genes across 6200 cells (GSE173385). Raw bone dataset included: a vehicle count matrix, with its features and barcodes with 18208 genes in 4629 cells; and a VD count matrix, features and barcodes with 19008 genes across 3914 cells (GSM6820990). Integrated prostate dataset included 13075 cells and 18535 genes in the expression matrix; bone dataset 8543 cells and 19516 genes; and skin dataset 15281 cells and 18545 genes.

2.3. Senescent Associated Gene Sets

Three senescence associated gene sets were retrieved for enrichment: SENmayo with 119 genes mainly associated with the senescent exosome [16], Geneage with 177 genes associated with overall aging [17] and Cellage with 949 genes associated with cellular senescence specifically [18].

2.4. Preprocessing

Statistical analyses were performed using R Statistical Software (v4.3.2; R Core Team 2023) and all plots were made with the ggplot2 package [19]. Seurat package was used for the management of the scRNAseq data [20]. First, all datasets were subjected to quality control evaluating the number, count, amount of genes and the percentage of mitochondrial genes. Threshold values were set as the number of genes to be greater than 200 and mitochondrial RNA transcripts of less than 7.5%. Subsequently, data was normalized with lognormalize, scaled by a factor of 10000 and further scaled for identifying relevant variable genes. Then, dimensionality reduction using principal component analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) using the first 10 first dimensions was performed. Finally, data was clustered using the KNN algorithm with PCA euclidean distances.

2.4.1. Annotation

Cellular type annotation was carried out with the SingleR package [21] using mouse prostate [22], skin [23] and bone [24] atlases as reference.

2.4.2. Integration

For batch effect and condition vs control (VD vs vehicle) integration, the Harmony algorithm included in the Seurat package (Hao et al. 2024) was applied. Finally, UMAP plots were obtained for cluster and cell type visualization. Annotation and integration raised questions regarding a cell type identified as immune in hanai (bone) dataset. Under further investigation, they were identified as myeloid precursors, consistent with what was reported in the original article of the dataset. Given the interest of studying cortical bone only and not bone marrow cells, they were removed from the dataset.

2.5. Functional Enrichment

Biological function enrichment of the processed datasets was performed using the AUCell algorithm (Aibar et al. 2017) via the Escape package (Borcherding et al. 2021) with the three senescent associated gene sets. Here, for gene conversion from human to mouse, the biomart package was used (Durinck et al. 2009). To see if there was a difference between enrichment scores of vitamin D and vehicle groups, a pairwise Wilcoxon test for each cell type and for overall tissue was made (Patil 2021) with a significance of p < 0.05.

2.6. Gene Regulatory Network Analysis

2.6.1. Network Processing

Gene regulatory networks (GRN) were made using the hdWGCNA package (Morabito et al. 2023). GRNs were acquired for the 3 overall tissues using a pseudobulk approach. Briefly, all the UMI were added according to the cell type and the intervention (vitamin D vs control) in a gene expression matrix. The matrix was normalized and subsequently, the soft powers and the co expression matrix were obtained. Finally, network topological variables like module eigengene connectivity were calculated for analysis. Worth noticing that, given the huge dropout in sc RNA seq data, the fraction of cells that a gene needs to be expressed in order to be included is set to 5%.

2.6.2. Network Analysis

Upfront analysis involved identifying hub genes and module eigengenes via topological analysis of the networks. Afterwards, differential module eigengenes expression (DME) between vitamin D and vehicle was assessed in every tissue in every cell type with more than 70 cells for significant statistical testing. Subsequently, to evaluate which genes were involved in the senescent phenomenon, an overlap analysis was performed between module genes and genetic senescent signatures via the ComplexUpset package (Krassowski 2020). Then, to assess the relevance of the overlapped genes in the modules, module eigen-centrality was evaluated for the overlapped genes. Later, GO enrichment of the modules and an intersection of the overlapped genes with said enrichment was carried out to evaluate in which biological processes of the module the senescent genes participated. This was done using the Enrichr package (Chen et al. 2013). Finally, DME of the overall tissue was obtained for each module. Moreover, to properly address the cellular response to VD in each module, the cell types with the minimal and maximal log fold change for each module were assessed.

2.7. Cellular Communication

Data analysis was performed using two complementary approaches, by looking at the interactions between senescent cells exclusively and the interactions between senescent cells and other cell types in the tissue. For this reason, cells of each tissue and intervention (VD vs vehicle) with a SENmayo enrichment score in the 90th quantile were clustered in a cell group named senescent_cells.

First, for each analysis, differential cellular interaction number and strength was investigated. Then, differential pathway analysis was performed. This was done in three ways: by manifold learning with functional analysis, by pattern recognition of dominant communication, and via information flow differences. All these analyses were performed with the cellchat package (Jin et al. 2023).

3. Results

3.1. Cell Annotation UMAP Between VD and Vehicle

UMAP plots of cell annotation (Figure 1) showed qualitative information regarding the cell type quantity and distribution. Comparison of the VD and vehicle plots show variances in the quantity of certain cell types. In the skin UMAP, we can visualize that the SB2 (suprabasal 2) cell type is more dense toward the end of the cluster shape in the VD plot, while in the vehicle plot is more dense the site of the intersection between SB2, SB1 and B cell types (SupraBasal 1 and Basal, respectively) (Figure 2B), suggesting that vitamin D promotes differentiation of the cell epidermis.

3.2. Functional Enrichment

SENmayo, geneage and cellage enrichment scores (ES) showed different trends in mouse organs (Supp. table 1). In the prostate and skin, SENmayo was significantly reduced in the vitamin D group (Figure 2A, 2C), while in bone, there was a slight increase in SENmayo score in the VD group (Fig 2B). Cellage difference between groups was significantly decreased in the prostate VD group, but in skin and bone, it was significantly increased in the VD group. Finally, Geneage score difference was significantly decreased in prostate and skin VD groups, whereas in bone it was decreased (Supp. table 1).

Regarding cell type analysis in each tissue, skin suprabasal 2,1 and hair follicle cells ES (with SENmayo) were significantly reduced by VD, while macrophages ES were increased. Subsequently, prostate T cell, B cell, macrophages and fibroblasts ES were significantly reduced with VD compared to the vehicle. Finally, bone chondrocytes, osteocytes and osteoblasts ES increased, while adipocytes were decreased (Supp. table 2).

3.3. Gene Regulatory Network Analysis

For the pseudobulk analysis, between 32-42 modules were obtained for each tissue (Fig 3A-3C). Remarkably, in prostate tissue, fibroblasts tend to have downregulated modules, while immune cells, particularly macrophages, tend to have upregulated modules (Fig 3A). In the skin, we can see a clear pattern of upregulation across all cell types except in basal and cycling keratinocytes, where modules are heavily downregulated. (Figure 3C). Finally, non parenchymal cells of the cortical bone (like tenocytes, myocytes and chondrocytes) had their modules upregulated, while parenchymal cells (like osteoblasts, osteocytes and osteoclasts) were downregulated. Notably, only green and pink modules are statistically significant among parenchymal cells, and they are upregulated in almost every other cell type (Figure 3B).

Regarding overlap analysis, in the prostate (Figure 3D), the senescent genes connectivity median is 0.40 (0.24-0.55). As for the senescent genes overlap DME (barplot in Figure 3D), two trends appear, 1) fibroblast downregulation/macrophage upregulation and 2) overall downregulation of the modules, with dark red and blue modules being non significant in the overall DME. Interestingly, 1) blue is downregulated in myeloid leukocytes and promoted in lymphoid leukocytes; 2) blue senescent genes downregulate myeloid immunity, degranulation and dendritic cell chemotaxis while they upregulate B cell differentiation, leukocyte aggregation, Natural killer T cell chemotaxis and NF-kB; 3) this may suggest that VD upregulates macrophage activity while it diminishes its aggregation and reduce lymphoid degranulation but enhance their chemotaxis and differentiation. For the other biological phenomena, vitamin D appears to upregulate DNA repair, replication, and transcriptional regulation; ECM organization and assembly; protein ubiquitination and neddylation; Insulin growth factor pathway and development; mitophagy, apoptosis and estrogen response in macrophages mainly, while reducing it in fibroblasts.

For bone, senescent genes median connectivity is 0.36 (0.27-0.47), lower than prostate. In the barplot we can appreciate that three modules are upregulated, colored magenta, salmon and blue, while the rest are downregulated (Figure 3E). Regarding biological processes, for the upregulated modules: magenta senescent genes seem to play a role in chondrocytes with the unfolded protein response (UPR) and gene expression pathways; salmon is involved in the interferon pathway although no cell type significantly changed this module; and blue mainly revolves around rRNA regulation and macroautophagy, which appears to be upregulated in chondrocytes and downregulated in tenocytes. For the downregulated modules, there is an opposing trend between tenocytes and myocytes, since in the former is inhibited the inhibition of transcription in response to stress and protein degradation, while in the latter is upregulated. Also, myocytes seem to be inhibited by vitamin D to undergo vascular reorganization.

For skin, senescent genes connectivity median is 0.29 (0.24-0.42). In the barplot, there is a clear contrary trend between suprabasal cells and basal cells, being the former upregulated and the latter downregulated, except for the green module which is overall downregulated (Figure 3F). For biological phenomena, the only upregulated module (blue) is in charge of activation of angiotensin-signaling, establishment of skin barrier and epithelium development. As for the rest of the modules, cytoplasmic translation, ECM reorganization, oxidative phosphorylation and autophagy are downregulated overall, although these processes seem to be cell dependent. For example, tan, darkgrey and brown are downregulated in all significant cell types (involved in translation, stress mediated response and peptidyl-serine auto-phosporylation processes) while yellow, lightcyan, blue and cyan are upregulated (involved in autophagy, oxidative phosphorylation and skin barrier function processes) in all significant cell types. The rest of the modules follow the basal-suprabasal contrary trend.

3.4. Cell communication analysis

Information flow between cells is vital to coordinate and maintain biological processes necessary for homeostasis of the organ and, ultimately, the individual. Given the characteristic exosome of the senescent cells, it is of interest to see how they are affecting each other and the other cells. For this, in Figure 4 and Figure 5 the main interactions, pathways and differential relations are shown between VD enriched and control tissues.

In regard to prostate senescent cell communication, senescent source fibroblast_Fbn1 and receptor luminal Clu and B cell interactions were upregulated while source luminal and fibroblast with receiving macrophages were downregulated (Figure 4A). For the functional analysis, 4 clusters based on similarity were obtained. Notice that the pathways with most distance (JAM and BMP) are in a completely different cluster (VD in violet and vehicle in blue) (Figure 4D). In Figure 4G Vehicle-enriched pathways emphasize inflammation (CD34, IL4, TNF, COMPLEMENT) and migration (SELPLG, ICAM, FN1), while vitamin D-enriched pathways focus on a more complex immune modulation (IL6, IL10, CD200), enhanced adhesion (CLDN, OCLN, NECTIN), and growth (PDGF, VEGF, IGF).

For bone senescent cell communication, only receiving pericyte interactions appear to be downregulated, while the rest senescent cell types interactions are upregulated, being the senescent chondrocytes the most upregulated. (Figure 4B). In the functional analysis, the pathways with the most distance are SEMA6, JAM and VEGF (Figure 4E). The pathways upregulated by VD are involved in vascularization (EGF), structural integrity (Collagen, fibronectin, tenascin), inflammatory response (TNF, CD80, SLURP) neural regulation (UNC5, SEMA6, SEMA7) and metabolic regulation (ENHO, Apelin, PTH and melanocortin). The ones downregulated were involved in bone formation and maintenance (IGF, SPP1, Periostin, BSP and RANKL), immune chemotaxis (CCL, CXCL, MIF), neural development (NGF, NTS, NRXN and SEMA4), cell adhesion and integrity (NCAM, CEACAM, CDH, NECTIN, and CLDN) and vascular regulation (ANGPTL and EDN) (Figure 4H).

Finally, in respect of skin senescent cell communication, the major upregulated senders were the fibroblasts (particularly dermal papilar fibroblast), while the most downregulated were the germinative layers of the hair follicle. For the receivers, notably downregulated cells were myocytes and macrophages, whereas the upregulated cells were mainly basal cells (Figure 4C). Functional analysis reveals that the most altered pathways were ADGRG and NOTCH (Figure 4F). Touching the overall signaling patterns, the flow information diagram shows that VD enriched pathways are involved in cell adhesion (DESMOSOME, LAMININ, CDH1, CLDN) and differentiation and development (EPGN, WNT, ADGRA) while it inhibits inflammation with pathways enriched in vehicle like COMPLEMENT, GAS and IL1, adhesion and angiogenesis (PECAM, VCAM). These factors suggest a tissue repair environment promoted by vitamin D.

Prostate analysis reveals that VD enhances certain senders like fibroblast cxcl9 and fbn1 and mesothelial cells, while it down-regulates mainly the rest of the fibroblasts and senescent cells (Figure 5A). For functional analysis, it can be appreciated that SPP1 and HSPG1 are the most differentiated pathways between conditions (Figure 5D)., The pathways enriched in VD were more related to extracellular matrix remodelling (NECTIN, ANGPTL, GAP, TENASCIN, FN1) differentiation, chemotaxis (VCAM, ncWNT, CX3C) and inflammation (CD40, TNF, PTN) while the control group is more associated with antigen presentation (MHCI). (Figure 5G).

For sender information, bone evaluation shows that cortical parenchymal bone cell (periosteum, osteocytes, osteoblasts) interactions were upregulated by VD, with fibroblasts, senescent cells and chondrocytes interactions being the most downregulated. Receiver information was markedly upregulated for periosteum and osteocytes and heavily downregulated for senescent cells and tenocytes (Figure 5B). In functional analysis, the top 2 pathways with major differences in distances between conditions were EPHB and ADGRL (Figure 5E). Conversely, the information flow shows that the majority of pathways are enriched in the vehicle and very few in VD. Notably, resorption pathways are upregulated by VD (RANKL, TNF, MMP, IL6) and downregulated mainly for chemotaxis (CXCL, PTN, CD34) (Figure 5H).

Finally, skin senders upregulated by vitamin D were senescent cells, germinal layer 2,3 cells and melanocytes while germinal layer 1 and cuticle cells were downregulated. However, for receivers, it was remarkably upregulated in the T cell population by vitamin D (Figure 5C), In the functional analysis, the pathways with the most distance or difference between conditions were CDH1, NECTIN and JAM (Figure 5F). Finally, we can see in Figure 5I that there are a lot of enriched pathways for VD, although with a minimal information flow, these correspond in grand majority to the T cell enrichment seen in Figure 5C. For the pathways with the most information flow, VD seems to upregulate antigen presentation (MHCI) chemotaxis and inflammation (CXCL, TNF, TWEAK, KLK) and matrix remodelling (TENASCIN, ANGPTL, CDH1).

4. Discussion

4.1. Preprocessing

Cell distribution and quantity across different lineages in different tissues can bring important clues to start single cell analysis. For instance, the differentiation pattern seen in the cell type UMAP plot of skin tissue (Figure 1B) is concordant with the current literature that vitamin D promotes keratinocyte differentiation (Hawker et al. 2007; Hu, Bikle, and Oda 2014). This has relevant implications in the aging process, since an impaired differentiation triggers a decreased permeability barrier function, which leads to more cutaneous inflammation and ultimately, aged related disorders (Bocheva, Slominski, and Slominski 2021; Wang et al. 2020).

Interestingly, Hu et al. shows that VD promotes keratinocyte differentiation, in part at least, by inhibiting the $\beta$-catenin pathway, as well as other authors (Ferrer-Mayorga et al. 2019; Jeong et al. 2015; Larriba et al. 2013). Accordingly, Figure 4F demonstrates that skin barrier establishment module is upregulated by vitamin D.

4.2. Functional enrichment

Firstly, to see if vitamin D affected cellular senescence across different tissues, the approach was to evaluate how much of the senescent gene sets were present with or without vitamin D. The results varied according to the senescent gene set used (Sup. table 1). However, senmayo has shown to be an optimal candidate for senescent cell identification using gene signatures (Saul et al. 2022; Tao, Yu, and Han 2024). Knowing this, the general trend was a downregulation of cellular senescence by vitamin D in prostate and skin, and an increase in bone. Current literature suggests that the role of vitamin D is to prevent and rescue senescence phenotype in fibroblasts, bone mesenchymal stem cells, smooth muscle, endothelium and chondrocytes (Chen et al. 2023, 2024; Lu et al. 2023; Marampon et al. 2016; Sayegh et al. 2024; Valcheva et al. 2014). Contrastingly, in the supplementary table 2, in bone, chondrocytes were enriched for senescence with VD. Delving deeper, communication between senescent chondrocytes and the other cell types involved mainly extracellular matrix remodeling (collagen deposition). To further corroborate these results, most perturbed modules in the GRN were reviewed for chondrocytes (Figure 3B). Interestingly, most downregulated ones mainly involved matrix remodeling, chondrocyte differentiation, oxidative phosphorylation and protein endoplasmic reticulum retention. Upregulated ones were about apoptotic immune clearance, glutathione synthesis and negative regulation of apoptosis. Thus, the apoptotic nature of the chondrocytes promoted by VD may contribute to the senescent profile of these chondrocytes.

Cell types that had a reduced senescent score with VD were prostate immune cells, bone adipocytes and suprabasal 1,2 keratinocytes. Vitamin D immunomodulatory effects are well described in literature (Ghaseminejad-Raeini et al. 2023; Mora, Iwata, and von Andrian 2008; Yeh et al. 2024) which concordates with these findings since adipocytes also play an important role in systemic inflammation via adipokines. Accordingly, in skin, as discussed previously, VD upregulates keratinocyte differentiation, which may intervene with the senescent phenotype and ultimately inhibit it.

4.3. Gene Regulatory Network Analysis

Secondly, to see how vitamin D affected cellular senescence across different tissues, the angle we took was to make GRN of the datasets and a differential module eigengene analysis to see how vitamin D affected the networks. We then identified the genes in each module that were involved in senescence with an overlap analysis between the module genes and the senescent gene signatures. Afterwards, we investigated what was the relevance of these overlapped genes in the GRN. Finally, we inquired about the biological processes in which the modules were involved, focusing on the ones that included the overlapped genes.

The first thing we noticed was that almost all modules with the most senescent genes in all the networks were overall downregulated by VD (Figure 3). This has interesting implications since vitamin D has been associated with anti-aging effects (Fantini et al. 2023; Martinelli et al. 2023; Zhu et al. 2019).

In the prostate network, the main findings were that Vitamin D upregulated the majority of the modules in myeloid leukocytes (macrophages and monocytes), it downregulated them in fibroblasts and diminished their overall expression. Current literature suggests that Vitamin D has anti-inflammatory properties in the macrophages (Song et al. 2016; Zhu et al. 2019), which agrees with our findings in the way that the blue module (involved with inflammation) is downregulated in almost every cell type. However, questions arose as one of the main processes in the blue module involved negative regulation of myeloid leukocytes, meaning that macrophage activity was upregulated by its inhibition. On further inspection, gene differential expression of the genes of this process revealed that they were actually upregulated by vitamin D, particularly Spi1, involved in macrophage lineage commitment and survival (Eichbaum et al. 1997; Zasłona et al. 2016) and suggesting an overall anti inflammatory effect of vitamin D. Moreover, as chronic inflammation is considered a hallmark of aging (López-Otín et al. 2023), here VD may counter this pathological state.

Interestingly, the most downregulated module by VD was red, which is involved in matrix reorganization in fibroblasts. This results agrees with the seen inhibitory effect of Vitamin D in cancer associated fibroblasts profile (Ferrer-Mayorga et al. 2017; Shany, Sigal-Batikoff, and Lamprecht 2016). Finally, turquoise module was the module with the most senescent genes. This module was involved in mitophagy, growth hormone (GH) receptor and splicing. VD upregulated this module in macrophages and downregulated in fibroblasts. This could be due to the fact that VD upregulates an anti-inflammatory profile via bioenergetic alterations, mainly involved with the control of the GH and mitophagy of worn out mitochondria in fibroblasts (Soler Palacios et al. 2023) and promote an immunomodulatory profile in macrophages. This process may also alleviate another hallmark of aging which is mitochondrial dysfunction.

Notably, prostate tissue had the senescent genes with the higher connectivity, even considering that this is from a PTEN knockout mouse tissue, or in other words, a precancerous model. This highlights the crossroads between aging and cancer, since they share many common mechanisms like genomic instability, telomere attrition, epigenetic alterations, immune dysregulation, mitochondrial dysfunction and cellular senescence (Havas, Yin, and Adams 2022). This also points out that vitamin D does influence the gene regulatory networks of both mechanisms significantly.

In the bone network, there are some senescent modules upregulated by VD. For instance, salmon involved IFN response. Delving deeper in this biological process genes, Bst2 gene turned out to be the most differentially expressed between conditions, and it was downregulated by VD in adipocytes mainly. This may explain how the adipocytes were downregulated by VD in the senescent gene signatures (Sup. Figure 2). For the overall trend of upregulation, this may be explained by the hypothesis developed by Newmark (Newmark, Dantoft, and Ghazal 2017) in which they explain the link with VD and IFN. Briefly, Vitamin D receptor (VDR) evolved to be a key regulator of immunity via sterol network control, another key player in this network evolved to be IFN in vertebrates, so the upregulation of one promotes the upregulation of the other, agreeing with our current results.

For the downregulated modules, interestingly, the module with the most senescent genes (turquoise) was involved in negative transcription under stress response, implying that VD upregulates stress transcription via c-JUN activation, which is key for the anti proliferative effects of vitamin D (Bi et al. 2016; Li et al. 2007). Another interesting module, although not included in the overlap analysis, is the green module, which remains as a very significant module in the network.

Interestingly, this module is involved in Signal Recognition Particle (SRP) pathways and is heavily upregulated in all cells, suggesting a strong interaction of VD with this SRP pathway. The SRP pathway is very important in the control of misfolded proteins, so is very involved in diseases associated with these pathophysiology, like neurodegenerative diseases, autoimmune myositis and cancer (Kellogg, Tikhonova, and Karamyshev 2022). Remarkably, in the work of Kimmel et al. of single cell transcriptomic analysis of aging in murine tissues, they mention that a common pathway involved in aging across tissues was the SRP pathway, underscoring the importance of this pathway in aging (Kimmel et al. 2019). This concordates with the current research that VD plays a major role in bone aging (Qiao et al. 2020). However, its definitive function remains elusive since the recent VITAL clinical trial has shown no significant reduction of fractures in non deficient nor osteoporotic patients after VD supplementation (LeBoff Meryl S. et al. 2022).

In the skin network, interestingly, work from Segaert suggests that VDR is actively expressed in cycling cells, and when arrested, VDR is decreased (Segaert, Degreef, and Bouillon 2000). This suggests that the downregulation of gene modules may be actively linked with VDR expression, particularly in senescent cells genes, which may hint the downregulation of senescent cells by VD to promote his effect on cycling cells. While for upper level keratinocytes, literature supports the VD role in the differentiation of these cells (Chaiprasongsuk et al. 2019), concordant with the results (Figure 3C). Another interesting result is that the downregulation of the green module (which has genes involved in inflammation pathways) is preserved across almost every cell type (Figure 3C), agreeing with VD immunomodulatory effects. All the biological phenomena regulated by VD agree with current literature. For example, the downregulation of oxidative phosphorylation in yellow module except in the epidermis (basal and suprabasal 2 cells) to reckon the stress induced by UV (Ambagaspitiya, Appuhamillage, and Dassanayake 2024; Pena et al. 2019), inhibition of protein and mRNA synthesis in brown and dark grey modules, suggesting the energy conserving profile of vitamin D, which is particularly important in the context of an anti-inflammatory profile.

Interestingly, upregulation of skin barrier function in blue module and ECM disassembly in purple suprabasal 1 cells appears, this supports the barrier function of VD in skin and the effect of VD on the hair follicle (Demay et al. 2007), since the most upregulated cell in the blue module was germ layer 2. Another remarkable finding is the inhibition of autophagy in lightcyan module except in suprabasal cells of epidermis and hair follicle, this may be because VD protects the cells in charge of maintaining the skin barrier (the suprabasal cells) but do not protect the rest of the tissue to avoid excessive energy consumption in other non vital sites for barrier function (and thus avoid energy demanding processes like inflammation). Altogether, vitamin D proves useful to maintain skin barrier function, provide anti inflammatory effects, and overall conserve energy for vital skin processes only. These main effects underscore the importance of vitamin D in cellular senescence and aging, since photoaging and inflammaging are crucial to develop a pathological aging; and to cellular senescence, since VD modulates SASP profile via anti inflammation, ECM regulation and energy conserving actions.

4.4. Cell Communication

For wrapping up our analysis, we performed a cellular communication network to see how senescent cells affected each other and other cell types in each tissue.

In the prostate network, the main findings indicated a reduced senescent fibroblast source- senescent macrophage receiver crosstalk in Figure 4A. Under further inspection, the main pathways involved in this communication involved collagen, agreeing with the information flow of this pathway seen Figure 4G and the reduced module expression of the fibroblasts of Figure 3A,D. However, a remarkable pathway was the APP one, since in the vehicle, the fibroblast express APP and the macrophages receive them via Cd74, but the APP pathway in the vitamin D group was nowhere to be found. While there is no evidence about the role of vitamin D on APP expression in prostate, APP has been shown to promote androgen-dependent growth of cancer (Takayama et al. 2008) and metalloproteinase expression, promoting tumor metastasis (Miyazaki et al. 2014). So this may suggest that VD attenuates prostate cancer through androgen dependent inhibition and cancer cell migration via APP. Other interesting finding between prostate senescent cells is that IL-6 and IL-10 are enriched in VD senescent cells. This suggests that these cells promote cellular senescence via SASP induction in the prostate and promote an antiinflammatory shift and matrix reorganization. This may be a protective effect since the activation of IL-6 can be useful to prevent cancer progression via senescence induction (Culig et al. 2016).

For the communication between senescent cells and the other cell types, we noted a significant downregulation of the sender communication of the top 10 percentile of senescent cells in the VD group. This suggests that VD downregulates the SASP, which agrees with current literature (Chen et al. 2024; Sayegh et al. 2024), however, under further inspection and in accordance with our previous findings with IL-6, TGFb and TNF were upregulated by VD in the senescent cells sender profile. This triggered investigation of these pathways, revealing that TGFb was shifted toward fibroblasts with high CXCL9, while TNF remained targeted mainly to macrophages and senescent cells in both groups. This was interesting, since CXCL9 fibroblasts communication was overall upregulated by vitamin D as seen in Fig 5A. These findings suggest that VD modulates SASP to improve immune function and recruitment via TNF, MIF and CSF pathways; and CXCL9 fibroblasts/macrophage interactions, but overall reduces SASP which is reflected in reduced collagen pathways as seen in FIg 5G.

Another interesting finding was the one with MHCI underrepresented in VD. We investigated the pathway and discovered that proliferating cells CD8 positive were the only receivers in the control group, and they were absent in the VD one. This could be explained by the VD immunosuppressive and anti proliferative characteristics (Cantorna et al. 2015). So, in conclusion, VD seem to modulate SASP to promote cellular senescence, immune recruitment and function of myeloid leukocytes, downregulate fibroblast activity (although not CXCL9 immune functions) thus diminishing ECM deposition, and enhance senescence profiles in existing senescent cells. This could be useful in early prostate cancer stages, however in more advanced cases it could prove deleterious.

For the bone network, in the senescent cell communication we found a marked increase in pericyte receiver information and an increase in the rest of the communication of the senescent cell types, particularly in chondrocyte (Figure 4B). On further dwelling, we found that senescent chondrocytes exhibited a curious communication profile, having some pathways like pre angiotensinogen, melanocortin and apelin with some information flow (Figure 4H). While there is some evidence of adropin’s role in alleviating aging in neurodegeneration (Banerjee et al. 2021) there is a lack of evidence regarding the role of adropin in cellular senescence. With apelin, there is information on alleviating angiotensin-induced senescence in endothelium via AMPK and SIRT1 activation (Yang et al. 2018). This agrees with our communication network, since apelin from chondrocytes is sent to the endothelium. This suggests that vitamin D attenuates cellular senescence, although the concentrations were high, which could raise the senescent signatures as seen in Figure 2B.

Another interesting finding was the downregulation of NOTCH in the vitamin D group as seen in Figure 4H. This pathway has been shown to induce secondary senescence in a juxtacrine way (Teo et al. 2019). Moreover, NOTCH inhibition in mice via myosin light chain 3 maintenance has been shown to inhibit chondrocyte senescence and osteoarthritis (Cao et al. 2023). Thus, vitamin D inhibition of NOTCH seems useful to delay senescence in cartilage.

In the communication between senescent cells and the other cell types, we noticed that parenchymal-cell interactions (mainly periosteum and osteocytes) were upregulated, while fibroblasts, chondrocytes and senescent cells were downregulated (Figure 5B). To further dwell in this characterization, we reviewed the interactions between vitamin D group and we noticed an increase in SASP expression, like IL-6, MMP and TGF-b (Figure 5H). Other pathways like SPP1, RANKL and netrin suggest that VD at high concentrations promotes SASP production to osteocytes, osteoblasts and periosteum to promote the resorption process mainly, while the SASP to other cell types was attenuated since they do not participate in bone resorption, although some pathways remained active like GAS and APP. This agrees with studies which mention that high doses of VD can promote loss of bone mass density and resorption (Burt et al. 2019; Nakamichi et al. 2018), and this may suggest that VD modulates SASP to play a role in resorption. GAS and APP maintenance in non parenchymal cell types from senescent cells could play a role in bone homeostasis promoted by vitamin D, promoting bone formation (Liu et al. 2023; Pan et al. 2018).

Finally, for the skin network, we discovered that communication between senescent cells was hampered in the Vitamin D group as seen in Figure 4C and in the downregulation of main SASP pathways like MMP, IL1, TGFb, GAS and CCL in Figure 4I, concordant with current literature (Sayegh et al. 2024). To further delve on these interactions, we studied in Figure 4I the main pathways altered by VD. Integrating this information, the first thing that we noted was the upregulation of fibroblasts communication and the shift from collagen, tenascin and gap ECM components to desmosomes and cadherins. This suggests that VD modulates SASP to a skin repair via epithelium regeneration rather than a fibrotic wound repair, which agrees with other authors results (Ge et al. 2022).

Interestingly, important contributors to the communication network were the germinal layers and the dermal papilla fibroblast cells, which concordates with current literature of the importance of VD in hair growth (Joko et al. 2023). Remarkably, this study by Joko et al. corroborated the importance of VD via VDR K.O. in the removal of senescence, but in this instance, in hair follicle generative layers specifically. However, senescence may have another role in this context of UV radiation, since germinative layer 2 senescent cells with vitamin D were overly expressing ANGPTL4, which has been associated with skin regeneration in wound process (Yang et al. 2023). Another interesting finding in Figure 4I was the upregulation of the galectin pathway in the VD group. Investigating deeper this pathway, we found out that the main cells receiving information were the macrophages. In literature, galectin is known to enhance macrophage phagocytic function (Karlsson et al. 2008) and polarize macrophages toward an immunomodulatory profile (M2) (Correa 2003) which could potentiate skin regeneration via vitamin D SASP modulation.

Subsequently, we investigated the senescent cell communication with other cell types and the overall changes of vitamin D in the radiated skin. We noticed that the senescent cells were the most altered source and T cells the most changed receiver in the network (Fig 5C). Further investigating both phenomena, main pathways sent from senescent cells involved matrix remodelling as the one seen in Figure 4I, concordant with what is seen in the overall pathway differential analysis in Figure 5I. Interestingly, all the fibroblasts were included in the senescent cells group, reflecting the fact that all fibroblasts went in the 10% of skin cells with highest enrichment SENmayo score. This may highlight the importance of cellular senescence in fibroblast-mediated skin repair (Demaria et al. 2014) and the way in which Vitamin D modulates its secretome to profile a dermal regeneration phenotype rather than a fibrotic process, which agrees with other works in the subject (Ge et al. 2022).

Another interesting approach was the notable upregulation of ANGPTL seen in Figure 5I. We observed that the only sources for this pathway in the VD group were the sebaceous glands, the upper hair follicle basal and the germinative layer cells. Concurrently, according to the work from Dahlhoff et al., ANGPTL4 proves to be important in reducing lipid droplet size and mass, suggesting ANGPTL4 mimics as therapeutic agents for acne (Dahlhoff et al. 2014), thus, this could provide some evidence for VD as a possible relieving agent in acne patients via ANGPTL activation (Wang, Zhou, and Yan 2021). For the role of ANGPTL and the hair follicle, various studies have shown that this pathway promotes follicle regeneration and hair growth via angiogenesis induction, thus, this could prove useful to promote VD as a possible treatment for hair follicle diseases.

5. Conclusions

Some perspectives and future directions can also be drawn form this work. Among these, we can mention the following:

• Vitamin D seems to play a role in senescent cell expression, particularly in the SASP profile, given the most significant change in SENmayo scores and cellular communication relevance in senescent cells across prostate, bone and skin.
• Gene regulatory networks of the 3 tissues show that almost all the modules with the most senescent genes were downregulated by VD, suggesting the importance of VD in senescent profiles and in aging.
• The role of VD in the SRP pathway in bone reveals trailblazing insight in aging bone, pointing out the relevance of future investigations to further characterize these findings.
• The maintenance of cellular senescence by vitamin D via IL-6, APP reduction, matrix reorganization and macrophage immunomodulation seems useful for early stages of prostate cancer, but inconclusive for later stages.
• Senescent chondrocyte profiles triggered by VD could provide clues to the role of these cells in energy homeostasis in bone.
• Paracrine senescence seems to be inhibited in bone by vitamin D via NOTCH downregulation.
• ANGPTL and galectin promotion by Vitamin D in senescent cells may play a crucial role in the radiated skin repair process via lipid metabolism and immune clearance.
• Senescent fibroblast SASP profile in irradiated skin appears to be relevant in extracellular matrix remodeling, and vitamin D seems to modulate it towards an anti fibrotic type.
• Further validation of these findings with other experimental methods are crucial for further understanding the role of vitamin D in cellular senescence.