Measurable Metrics of Mesenchymal Stem Cell Aging

Introduction

Preparations of allogenic and autologous mesen­chymal stem cells (MSCs) are considered as a promising component of regenerative cell therapy [1]. In clinical trials, high doses reaching hundreds of millions of cells per one procedure are used. It often requires the expansion of the culture in vitro, which, together with the compromised functional state of the donor’s organism, may reduce the quality of cell preparations [2–4]. Functional methods of MSC preparation for clinical application imply the assessment of morphology, cell survival, differentiation potential, and biological safety [5, 6]. At the same time, additional investigations are needed to study the effect of functional state of the donor’s organism, cultivation conditions, and duration on the quality of cellular products. When classifying the markers of aging, various functional manifestations such as metabolic and immune disorders, genome instability, and epigenome alterations are distinguished [7]. One of the hallmarks of aging is the accumulation of the senescent cells unable to proliferate, resistant to apoptosis, and possessing the characteristic morphological and metabolic phenotype [3, 7].

Inducers of cellular senescence may include the exhaustion of proliferative potential, accompanied by critical telomere shortening; exposure to toxic, genotoxic, and oxidative stress; induction of oncogenes; inflammation; mitochondrial dysfunction; disruption of epigenetic regulatory mechanisms, and other factors [3, 8–11]. It is important to note that the listed factors may also be the secondary effects of cellular aging, and its phenotypic manifestations at the cellular level can vary in a wide range.

Due to its heterogeneity both in vitro and in vivo, there is no specific universal marker of cell senescence [12, 13]. Therefore, the investigation of cellular aging dynamics in general relies on the analysis of several markers, whose combination is inherent to this process [3, 12]. These markers include the induction of senescence-associated β-galactosidase, activation of the cell cycle inhibitors p16^INK4A, p21^CIP1; reduced expression of LMNB1 and HMGB2 proteins, which shape the structure and architecture of the cell nucleus [14–16]. In addition to the analysis of the mentioned markers, functional tests can be performed to assess proliferative potential, presence of the DNA damage or apoptosis markers [12, 17].

The characteristic feature of the senescent cells is secretion of proinflammatory cytokines, chemokines, growth factors, and proteases, which compose specific senescence-associated secretory phenotype (SASP) [18]. The detection of SASP factors serve as an indicator of cell aging, however, their abundance varies significantly and depends, in particular, on functional cause of the cell senescence [11]. Thus, a classic approach to the exploration of cell senescence is based on the analysis of sufficiently wide spectrum of non-exclusive (non-specific) markers and conducting functional tests. At the same time, the perspective integral assessment of cell senescence by predictive models built on the analysis of DNA methylation patterns, transcriptomic data, and cell morphology is presently being developed [19–21]. These predictive models may consider variability of cell aging phenotype and usually depend, to the lesser extent, on separate markers, which makes them a promising analytical tool.

The aim of the study is to analyze manifestations of selected markers of senescence on the models of replicative senescence, stress-induced senescence, and chronological aging of human mesenchymal stem cells. Among the markers we assessed were the level of expression of individual genes whose expression dynamics are associated with aging, global transcriptome alteration during chronological and in vitro cellular aging; telomere length measurement and changes in cell morphology and manifestations of specific cytological markers of cellular aging. In this study, we evaluated the applicability of these markers for assessment of MSC aging as well as the limitations of methods used that could potentially bias the analysis. In addition, we have analyzed the conceptual feasibility of building predictive models for estimating chronological age and duration of in vitro cultivation based on the transcriptomic data and cell morphology analysis.

Materials and Methods

Cell cultures. Cell samples were obtained from the donors with prior written informed consent. The study was approved by the local ethics committee of the Federal Center of Brain Research and Neurotechnologies of the Federal Medico-Biological Agency of Russia (Protocol No.7-5-22 of September 6, 2022).

In the study, MSCs (n=2) isolated from the Wharton’s jelly of the healthy pregnant woman (38–40 weeks of gestation); bone marrow derived MSCs (BM-MSCs) from healthy donors aged 18–25 years (n=3) and donors older than 65 years (n=3) obtained from the mononuclear cell fraction of bone marrow, which was isolated by gradient centrifugation (20 min, 400 g) in the ficoll solution (PanEco, Russia). The cells were cultivated in the DMEM/F12 medium (Servicebio, China) supplemented with 15% fetal bovine serum (Capricorn, Germany) and the antibiotic cocktail of penicillin (100 units/ml) and streptomycin (100 µg/ml) (Gibco, USA). Subculturing was performed at split ratio of 1:4. Using flow cytometry, the expression of the following MSC markers was analyzed: CD29, CD44, CD73, CD90, CD105, CD34, CD45 (FITC/PE; Miltenyi Biotec, Germany), and HLA-DR. The cells exhibited morphology and immunophenotype characteristic of MSCs: CD29⁺, CD44⁺, CD73⁺, CD90⁺, CD105⁺, CD34^–, CD45^–.

For induction of the stress-induced cellular senescence, MSCs were cultured until they reached 60% of confluence, after which the culture medium was replaced with medium containing 200 µM hydrogen peroxide (Dia-m, Russia). After 4 h of incubation, the medium was removed, and MSCs were washed twice with the phosphate-saline buffer. Then the cells were incubated under the standard cultivation conditions for 3 days, after which they were used for further analysis.

Immunostaining and cytochemical analysis of senescence-associated β-galactosidase. Cells were cultured in 96-well plates for confocal microscopy (SPL Lifesciences, South Korea) or on cover glasses precoated with 0.1% gelatin solution (Sigma-Aldrich, USA). Upon reaching the required confluency, the samples were fixed in 4% formaldehyde solution (Sigma-Aldrich, USA). The activity of senescence-associated β-galactosidase was analyzed using the previously described method [22]. For immunostaining, the cells were incubated in the 0.1% Triton X-100 solution (Amresco, USA) for 30 min, after which they were incubated in 1% BSA solution (Sigma-Aldrich, USA) for 1 h. The following primary and secondary antibodies were used for immunostaining: Ki-67 (Cell Signalling Technology, USA or Milteny Biotec, Germany); H3K9me3 (Active Motif, USA); Donkey Anti-Mouse IgG H&L (Alexa Fluor® 488) (Abcam, Great Britain); Goat Anti-Rabbit IgG H&L (Alexa Fluor® 568) (Abcam, Great Britain); Goat anti-Rabbit IgG (H+L) (PE-Alexa Fluor™ 647) (Invitrogen, USA). Hoechst 33342 (Invitrogen, USA) was used for nuclear staining. The samples were analyzed using the Olympus BX 51 fluorescence microscope (Olympus Corporation, Japan) and Nikon A1 scanning laser confocal microscope (Nikon Corporation, Japan).

Gene expression analysis using real-time PCR. The Rizol reagent (diaGene, Russia) was used for RNA isolation following the manufacturer’s protocol. The complementary DNA was synthesized with reverse transcription reagent kit (Biolabmix, Russia). The real-time PCR was performed using BioMaster UDG HS-qPCR SYBR Blue premixes (Biolabmix, Russia). The ACTB and SDHA were used as reference genes for normalization.

Below are the sequences of oligonucleotides used in our work:

ACTB_F ACAGAGCCTCGCCTTTG, ACTB_RCCTTGCACATGCCGGAG;

SDHA_F TTTGATGCAGTGGTGGTAGG, SDHA_R CAGAGCAGCATTGATTCCTC;

p21_F TGGAGACTC TCAGGGTCGAAA, p21_R GGCGTTTGGAGTG GTAGAAATC;

HMGB2_F CTTGGCACGATATGCAGCAA, HMGB2_R CAGCCAAAGATAAACAACCATATGA;

LMNB1_F ACACTTCTGAACAGGATCAACC, LMNB1_R CTGTGACACCAGCGTTTGC;

p16ink4a_F CCCAACGCACCGAATAGTTA,
p16ink4a_R ACCAGCGTGTCCAGGAAG;

IL6_F GTGGCTGCAGGACATGACAA, IL6_R TGAGGTGCCCATGCTACATTT;

IL8_F AAGAGCCAGGAAGAAACCACC, IL8_R CTGCAGAAATCAGGAAGGCTG;

IL1b_F CTGTCCTGCGTGTTGAAAGA, IL1b_R TTGGGTAATTTTTGGGATCTACA;

PAI1-F CTCATCAGCCACTGGAAAGGCA, PAI1-R GACTCGTGAAGTCAGCCTGAAAC;

MCP1_F CTTCTGTGCCTGCTGCTCATA, MCP1_R CTTTGGGACACTTGCTGCTG;

MMP1-F TGGACGTTCCCAAAATCCTG, MMP1-R AAGGGATTTGTGCGCATGTAG;

MMP3-F CTGCTGTTGAGAAAGCTCTG, MMP3-R AATTGGTCCCTGTTGTATCCT.

Measurement of the telomere length using real-time PCR. The real-time PCR was performed with BioMaster UDG HS-qPCR SYBR Blue premixes (Biolabmix, Russia) using primer pairs Tel-F CGGTTTGTTTGGGTTTGGGTTTGGGTTTGGGTTTGGGTT and Tel-RGGCTTGCCTTACCCTTACCCTTACCCTTACCCTTACCCT specific to telomeric DNA repetitive sequence as well as 36B4u CAGCAAGTGGGAAGGTGTAATCC and 36B4d CCCATTCTATCATCAACGGGTACAA specific to the region of acidic ribosomal phosphoprotein P0 gene on chromosome 12. Genome copy numbers and the total length of the telomeric DNA were determined relative to the DNA strand, represented by the plasmid pAL2-T (Eurogene, Russia) bearing 36B4 genomic region or human telomeric DNA fragment encompassing 14 repetitive units with the total length of 84 bp. The following parameters of amplification were used: 50ºС for 5 min, 95ºС for 10 min, 39 cycles of 95ºС for 15 s and 60ºС for 1 min.

Analysis of the next-generation transcriptomic sequencing data and building predictive models based on the transcriptome data. The datasets deposited in the Gene Expression Omnibus repository (GSE139073, GSE145008) were used in our work [23, 24]. Short-read alignment to the reference GRCh38 genome, preprocessing, and detection of gene expression were performed using STAR, SAMtools, and featureCounts programs [25–27]. To eliminate batch effects in the transcriptome data, the ComBat-seq method implemented in the sva package was applied [28]. The statistical edgeR package was employed to analyze the differential gene expression [29].

The predictive models of chronological age and in vitro cultivation duration were created based on the normalized values of gene expression. Genes, whose expression level correlated significantly with the passage or chronological age of the donors, were selected using Spearman and Pearson coefficients of correlation (coefficient >|0.5|, p-value-adjusted <0.05). The regressive models were built using LASSO regression and random forest regressor (RFR) from the Scikit-learn package [30]. Data were divided into two sets: the training set (80%, 84 sequencing samples) and the test set (20%, 22 sequencing samples). For the LASSO regression-based model, automatic tuning of the hyperparameter was applied with LassoCV on the training set. For the RFR model, the base number of tree parameters was used. The model quality was evaluated on the test set which was not involved in the learning process.

Data preparation and training of the neural network segmentation model. In the first stage, using the segment anything image-recognition model followed by manual validation and correction, nuclear masks were generated for microscopic images of cell preparations of umbilical MSCs at different culture passages (passage range — 3–15, a total of 27,500 cells) [31]. At the next stage, images were scaled to the equal resolution and divided by a sliding window with a 246-pixel pitch into the overlapping fragments 256×256 in size. The window step provided an overlap of neighboring image fragments by 10 pixels on each side, which reduced boundary artifacts during the subsequent assembly of the final segmentation map. To increase the model’s robustness to various exposure and contrast variations of the images, augmentation methods were employed. Among the transformations used were horizontal and vertical flips, random adjustments of brightness and contrast, as well as scaling with small shifts. The final dataset comprised 563 examples and was split in an 80 to 20% ratio for training and testing, providing a sufficient number of samples for proper tuning of the network parameters.

For solving the cell nucleus segmentation task, a convolutional neural network architecture, DeepLabV3+ [32] was used. As the backbone network, EfficientNet-b0 [33] pretrained on the ImageNet dataset was selected, providing the models with an initial representation of low-level image features. Training was performed for 40 epochs, allowing the model to reach stable convergence. During this period, the model was trained on a compute node equipped with an NVIDIA A100 GPU, completing the full training cycle in 3.5 h. To minimize the impact of class imbalance (significant differences in nucleus sizes and thin boundaries) and achieve more accurate segmentation, a combined loss function was used, which integrated two components: the BCE-Dice Loss provided high sensitivity to the imbalance between classes (nucleus/background) and accounted for spatial consistency of predictions; the focal loss improved training by reducing the influence of easily classified examples. During training, a learning-rate scheduler was employed, adjusting the learning rate from an initial value of 1e–3 down to 1e–5 after each iteration, ensuring a gradual reduction of the optimization step and promoting stable convergence of the model.

For calculating quantitative morphometric charac­te­ristics of nuclei, a binary mask obtained from the DeepLabV3+ segmentation results was passed to the analysis function. Before the calculations, pixels marked as “border” were excluded from the overall mask, after which sequential erosion and dilation (by 20 pixels) were performed to remove thin artifacts and merge broken contours. The parameters computed for each nucleus: center coordinates (X, Y), area, roundness, semi-major/semi-minor axes of the ellipse and inclination angle, the Hausdorff distance. The coefficient of belonging to the class was also established: class 1 (passages 3–5), class 2 (passages 7–9), class 3 (passages 11–15). The executable scripts of the model are deposited in the GitHub repository (https://github.com/LabADTCellSeg/cellseg).

Results

Senescence-associated changes in gene expression profile. To assess the senescence-associated gene expression dynamics in MSCs we selected genes that might be considered as principle regulators of the cell cycle and nucleus structure, as well as genes encoding components of the proinflammatory phenotype. We examined the expression levels of the following genes — P16INK4a/СDKN2A, P21CIP1/CDKN1A, LMNB1, HMGB2, IL6, IL8/CXCL8, IL1B, SERPINE1/PAI1, MCP1/CCL2, MMP1, MMP3. Expression analysis was performed on a replicative senescence model: independent umbilical cord MSCs subjected to long-term culture (n=2); on a chronological aging model: BM-MSCs from the donors of different ages: 20–25 years (n=2) and older than 65 years (n=2); on the model of stress-induced cellular senescence: umbilical cord MSCs exposed to hydrogen peroxide to induce cellular senescence (n=2).

The genes encoding cyclin-dependent kinase inhibitors (СDKN2A and CDKN1A) displayed a similar dynamics across all examined samples. Their expression increased during prolonged culture, in response to oxidative stress and with increasing donor age (Figure 1). Notably, the increase of expression was more pronounced for CDKN1A, whereas changes of CDKN2A expression were weaker and, in the context of the replicative aging, insignificant. The expression level of the genes encoding the nuclear architectural proteins LMNB1 and HMGB2 consistently decreased both in replicative and stress-induced senescence models as well as during chronological MSC aging (see Figure 1).

The most dramatic changes were observed during replicative senescence, where the decline in expression over the culture period reached about 90% for both LMNB1 and HMGB2. It is noteworthy that genes encoding the SASP components exhibited different dynamics depending on the cause of cellular senescence (see Figure 1). For example, peroxide-induced senescence led to an increase in the expression levels of all studied genes except MCP1, whereas replicative senescence did not activate the genes encoding matrix metalloproteinases MMP1 and MMP3. It is also worth noting that analysis of the BM-MSCs from the donors of different ages did not reveal age-related changes at the expression levels of the cytokine genes IL-6, IL-1β, CXCL8 was weakly activated in the cells of the elderly donors. The SERPINE1/PAI1 gene, encoding PAI-1 protein, appeared to be most stable in the context of aging-associated dynamics of gene expression.

Various manifestations of senescence-associated proinflammatory phenotype depending on the type of cellular senescence are generally expected. Nevertheless, to independently verify the obtained results, we searched for the publicly available high-throughput transcriptome sequencing datasets of cultured MSCs that include donor age and cell passage information in the Gene Expression Omnibus repository. As a result, a dataset (n=37, age 3–85 years, median value 47 years; Table 1) has been prepared and correlation analysis of expression changes in the tested genes with prolonged cultivation and donor age has been performed. According to the obtained results, SERPINE1/PAI1, CDKN1A, and CDKN2A genes demonstrated a positive correlation with culture duration, whereas LMNB1 and HMGB2 showed a negative correlation (Spearman correlation, p<0.05; Table 2). The correlation with the donor age was detected for CDKN1A and LMNB1 genes when no correction for multiple testing was applied (Table 3). Differential expression analysis in BM-MSC samples from donors aged 20–35 years (n=7) and donors aged 60–85 years (n=13) has identified only 50 differentially expressed genes (|logFC|>2, p<0.05) (Table 4). It is highly probable that donor-dependent variability of gene expression patterns in MSCs can be quite substantial, making it difficult to identify transcriptional markers of chronological aging. Nevertheless, the expression of individual genes might be applicable for assessing cellular senescence in vitro.

To conceptually validate the applicability of transcriptomic data analysis for predicting donors’ chronological age or the duration of cell cultivation, we prepared corresponding predictive regression models based on two approaches: LASSO regression and RFR. Genes were selected as predictors based on Pearson and Spearman correlation coefficients, in order to account for both linear and monotonic relationships between features and the target variables — chronological age and cultivation duration. As a result, the models based on LASSO regression and RFR have demonstrated close performance: R²=0.755; MAE=9.858 years, and R²=0.742, MAE=10.060 years, respectively (Table 5). The LASSO regression-based model has demonstrated the highest accuracy of the cell passage prediction on the test sample: R²=0.583; MAE=0.508 passages (Table 6).

Analysis of telomere length for cell senescence assessment. One of the traditional markers of assessing cellular aging is the analysis of telomere length. For this analysis we used the real-time PCR method. Samples of umbilical cord MSCs of different passages and BM-MSC samples from the donors aged 20–25 years (n=2) and older than 65 years (n=2) were analyzed on the sixth cultural passage. This method allowed to detect dynamics of the telomere shortening during MSCs cultivation. Statistically significant differences were observed after seven passages (Figure 2). At the same time, when samples from the donors of various ages were compared, no reliable differences were found.

Analysis of nuclear morphology as a marker of cell senescence. To evaluate the dynamics of the nuclear morphology and the expression of individual protein markers during the cell aging, we have analyzed BM-MSC samples of three donors from the young (18–25 years) and older age groups (over 65 years); umbilical cord MSCs exposed to continuous cultivation; and MSCs treated with hydrogen peroxide to induce senescence. In the cytochemical study of senescence-associated β-galactosidase activity (Figure 3), an increase in its activity was observed during replicative and stress-induced senescence (Figure 3, (b), (c)). When comparing MSCs obtained from the donors of different ages, the differences were not pronounced (Figure 3 (a)).

Immunostaining analysis of the proliferation marker Ki-67 (Figure 4 (а)–(c)) allowed us to detect a relative decrease of the number of Ki-67-positive cells associated with cultivation duration (Figure 4 (a), (d)).

Moreover, we did not find significant differences comparing cell preparations from the donors of various ages (see Figure 4 (c), (d)). However, at the level of nuclear morphology, the reduction of H3K9me3 signal intensity was noted, which agrees with heterochromatin erosion observed in aging, and enlargement of nuclei was observed at later passages and with increase of the donor age. Similar effects were also noted during continuous cultivation and stress-induced senescence of the umbilical cord MSCs (see Figure 4 (а)–(c)).

It is interesting to note that an increasing amount of evidence is accumulating in favor of using changes in nuclear morphology as an independent marker of cellular senescence. From the technical point of view, this analysis looks robust, since it actually requires microscopic analysis coupled with fluorescent nuclear staining. Moreover, available approaches in the machine learning image recognition accelerate processing and collection of the required statistical data. For this reason, to systematically analyze changes in nuclear morphology during aging, we trained a segmentation neural network model that describes nuclear morphology parameters such as area, roundness, and ellipse parameters (see “Materials and Methods”). The model demonstrated high segmentation quality on the test set (Figure 5 (a)). As the primary metric for evaluating model performance, the Intersection over Union (IoU) measure was used. The IoU obtained for the model, equal to 0.88, indicates effective segmentation of cell nuclei and their boundaries even in the presence of noise and variations in the original images. For each nucleus, a set of parameters was calculated:

1) the center coordinates (X, Y) allow the nucleus to be matched with other cellular structures and used for spatial analysis; they are defined as the centroid derived from contour moments;

2) area — a marker of the overall size of the nucleus, an increase in area can be associated with cells transitioning to later passages;

3) roundness characterizes shape compactness, an increase in the value indicates that the nucleus is acquiring a more circular shape, which may be associated with later stages of cellular senescence;

4) the ellipse semi-axes and orientation angle allow assessment of the degree of elongation and the nucleus orientation;

5) the Hausdorff distance serves as a criterion of approximation quality, low values (<10) indicate that the nucleus shape conforms to an ellipse.

Using the developed model, data on the size and shape of MSCs nuclei (Figure 5 (b)–(d)) from the donors of various ages (n=9308), subjected to varying the duration of in vitro cultivation (n=5157), and before and after the induction of stress-mediated cell senescence (n=564) were collected. Cells from donors in the older age group showed a broader range of nuclear sizes at the early stages of culture, whereas with increasing passage the nuclear sizes of cells from donors of different ages converged (see Figure 5 (c)). During prolonged culture the nuclear size also increased significantly, and this increase occurred gradually (see Figure 5 (b)). Stress-induced senescence was associated with the most dramatic increase in nuclear size (see Figure 5 (d)).

At the next step, we evaluated the possibility of predicting the duration of cell cultivation, expressed as the ordinal number of the cultural passage, from nuclear morphology. The developed model classified the cells into three classes: early cultural passages (passages 3–5), intermediate (passages 7–9), and late (passages 11–15) according to the structure of the training dataset. As a result of the model’s operation, the analyzed cells were assigned a membership coefficient for one of the listed classes. The algorithm for calculating class-membership coefficients comprised several stages. At the first stage, a cell region was generated by merging the boundary mask with the class masks; then, based on the segmentation results, the cell contour was extracted and filled. The next stage involved counting overlaps with the class masks to determine the number of pixels that simultaneously lay within the cell region and the mask of the corresponding class. Subsequently, normalization by cell area was performed, whereby for each class the ratio of overlapping pixels to the total cell area was converted into a proportion ranging from 0 to 100%. These proportions were then adjusted so that their sum equaled 100%, and the corrected proportions were taken as the membership-coefficient values for each class. Membership to the primary classes (1, 2, 3) was defined by the highest coefficient value; a mixed class (1/2, 2/3, 3/1) was assigned when the difference between the coefficients of two classes did not exceed 10%; an undefined class was assigned to a cell when the coefficients for all classes were close to each other (all values below 40%). Evaluation of the algorithm on the test set demonstrated adequate prediction of the actual culture passages of the studied cells (Figure 5 (e)). Starting from passage 9, the representation of cells in different classes increased, which may be related to rising morphological heterogeneity of cells associated with aging. With increasing passage number in the test sample, the proportion of cells classified as late-passage cells steadily grew, most likely reflecting the dynamics of accumulation of senescent cells.

Discussion

The key hallmark of cellular senescence is an irreversible arrest of the cell cycle mediated by the activation of cyclin-dependent kinase inhibitors p16^INK4a and p21^CIP1a [3]. In the studied models of MSC cellular senescence, the genes encoding the proteins p16^INK4a, p21^CIP1a also displayed the expected dynamics. Moreover, substantial changes occurred in cell morphology and nuclear architecture. A significant role in the alteration of nuclear and chromatin structure is played by the senescence-associated decrease in the expression level of LMNB1 gene [15, 34]. The reduction of the LMNB1 gene expression level was observed in all examined models of cellular aging and it was most pronounced in the replicative senescence of MSCs. Specific changes in chromatin structure also include the formation of the so-called senescence-associated heterochromatin domains (SAHF and SAHD) together with decondensation of peri/centromeric heterochromatin regions (SADS) and global erosion of heterochromatin [16, 35–37]. These changes are involved in the regulation of both genes comprising the so-called senescence-associated proinflammatory phenotype [14, 35].

One of the factors determining the chromatin structure is a nuclear architectural protein, HMGB2, whose expression declines during cellular senescence and which, in particular, is considered as an early marker of cellular senescence [14, 38]. A decrease in the expression of the HMGB2 gene was also detected with the increase of chronological age of MSC donors, in stress-induced and replicative senescence of MSCs. In the oncogene-induced senescence model HMGB2 has been shown to prevent propagation of heterochromatin in the genome regions containing genes forming the so-called senescence-associated proinflammatory phenotype, thereby helping to maintain their expression [14]. Although such an effect has not been confirmed in the replicative senescence, however, the development of a more permissive chromatin state during aging — caused by disruptions of the machinery that maintains facultative and constitutive heterochromatin, remodeling of the nuclear spatial topology, and activation of intracellular pro-inflammatory signaling pathways in response to DNA damage (cGAS-STING) — are key factors that determine the formation of the SASP [15, 37, 39].

The concept of senescence-associated secretory phenotype unites the complex of proinflammatory cytokines, growth factors, and metalloproteinases [40]. The main factors that constitute the SASP are TNFα, MCP-1, MCP-2, SERPINE1/PAI-1, GM-CSF, GROα, β, γ, IGFBP-7, interleukins IL-1α, IL-6, IL-7, IL-8, chemokine MIP1α, and matrix metaloproteinases MMP-1, MMP-10, and MMP-3 [41]. However, it is important to note that the composition of SASP changes significantly depending on the cause of cellular senescence and the cell type [18]. Interestingly, when comparing MSCs from donors of different ages, no significant dynamics in the expression of genes encoding individual interleukins was detected. At the same time, during culture and under stress-induced senescence, pro-inflammatory SASP factors such as IL-6, CXCL8, IL-1β became activated. This observation partially contradicts the previous reports of increased activity of these SASP factors in MSCs from older donors [42]. It should be emphasized, however, that in the cited study, a convincing difference in expression was demonstrated only for IL-6. Moreover, our analysis of published transcriptomes of BM-MSC samples (n=37) did not find reliable correlation between the expression of the investigated genes encoding individual SASP components and the age. All this together may reflect the heterogeneity and substantial contribution of donor-specific effects that complicate the analysis of age-related changes of gene expression. This is also supported by the performance of the regression model for predicting chronological age, which exhibited a relatively high value of the mean absolute error (R²=0.755; MAE=9.858 years).

It should be noted that previously described predictive models for estimating age from transcriptomic data have demonstrated comparable effectiveness [20, 43, 44]. More accurate similar algorithms generally achieve maximum performance on the narrow age cohorts. Moreover, in the process of model development, the authors excluded multiple available samples from the analysis, since their inclusion significantly worsened the model quality [45]. Thus, the evaluation of gene expression dynamics may be used to the greater extent to analyze cellular senescence in vitro. In this case, typically only limited number of cell lines are investigated under relatively standard cultivation conditions, which is likely to reduce the variability of gene expression profiles inherent to the primary donor cell cultures and samples. Similarly, according to the data obtained by us, the assessment of MSC telomere length, at least on the small sample sets, is also rather applicable for the evaluation of replicative senescence in vitro.

Cytological analysis techniques are widely used to study cellular senescence. Cellular senescence is accompanied by characteristic morphological changes such as flattening, enlargement of the cell and nuclear size, as well as the appearance of specific protein markers, like activation of senescence-associated β-galactosidase [3, 46]. In the present work, the activity of senescence-associated β-galactosidase demonstrated its applicability for the qualitative assessment of both replicative and stress-induced senescence. However, the use of this marker for the evaluation of the functional state of the cells requires standardization of several conditions. First, it is necessary to control the efficiency of the reagent lot used in the work, since the pH shift of these reagents may essentially distort the results. Taking into account the necessity to analyze the freshly prepared cell slides, it is not always feasible in serial experiments conducted in research laboratories. Besides, cell preparations must demonstrate similar cell density, since the elevated confluence can lead to distorted results [47]. Interpretation of the obtained results at the early stages of the cellular senescence may be difficult due to the absence of a fixed threshold value of the β-galactosidase activity, making it difficult to classify cells as positive or negative for this marker. Altogether, this limits the application of this marker for studying cell preparations obtained from donors of various ages.

As an alternative, approaches to assessing cellular aging based on the analysis of several protein markers associated with proliferation, apoptosis, and DNA damage may be considered [17]. At the same time, the applicability of these approaches is again limited by the selection of optimal markers. For example, expression of the protein Ki-67 widely used as a proliferation marker depends on the stage of the cell cycle, while the variant of γH2Ax histone, the marker of the DNA damage, is detected at the late stages of cellular senescence [17, 48].

With the development of the machine learning methods, there was proposed a concept, according to which the processive analysis of cell morphology may serve as an integral metric of cell aging [21, 49, 50]. To estimate the applicability to this approach, we have developed a segmentation neural network model for the automated assessment of the nuclear morphology parameters. The analysis of BM-MSCs from the donors of various ages at the early passages has shown that the range of nucleus sizes was wider in the cell sample of the donors from the older age group. The nucleus size gradually increased during MSC cultivation and the size of the BM-MSC nuclei from the donors of various ages did not demonstrate significant differences between the age groups. The most prominent change of the nuclear morphology was observed in stress-induced cell senescence. The developed model also allows for effective classification of the individual cell passage as a surrogate metric of the cell aging stage for umbilical cord MSC samples. In this connection, it may be supposed to be also employed for the estimation of the functional state of the donor BM-MSC samples if there is a sufficient amount of datasets for the training sample. Besides, the application of similar models to assess the effects directed to the reduction of cell senescence manifestations, such as rejuvenation by partial reprogramming is of great interest [51, 52].

Conclusion

In the presented work, some aspects of phenotypic manifestations of various types of MSC senescence have been studied. At the level of individual gene expression, it has been shown that the change in the expression levels of CDKN1A, LMNB1, HMGB2, and SERPINE1/PAI1 is observed in all investigated models of cellular senescence. At the same time, the analysis of transcriptomic data has demonstrated significant donor-dependent heterogeneity of gene expression profiles, which hampers creation of effective predictive models for the evaluation of chronological age and the duration of the in vitro cultivation. At the same time, an alternative predictive metric of cellular aging — at least in the case of replicative aging — can be changes in nuclear morphology, whose dynamic analysis using neural-network models allows us to estimate the duration of in vitro cultivation. Combining such approaches with other promising metrics, such as epigenetic clock algorithms, gives hope for developing functional algorithms to evaluate the phenomenon of cellular senescence.

Research funding. The work was financially supported by the Russian Science Foundation, Research Project No.22-74-10123.

Conflicts of interest. The authors declare no conflicts of interest.