Validation of the datasets and identification of DEGs

To investigate intra-group data repeatability, PCA was employed to reveal the data distribution in each sample to ensure that the data structure (young and old) was appropriate (accurate and reliable). As shown in Fig. 1A, the evaluation showed that the data from GSE39035 could be unambiguously divided into two groups based on age. Likewise in the GSE97311 dataset (Fig. 1B), the data could be divided into young and old groups according to age.

Fig. 1

Principal component analysis (PCA) and Screening of DEGs. A. PCA of DEGs between samples of MSCs from young (< 65 years of age) and old (> 65) donors in the GSE39035 dataset. B. PCA of DEGs between samples of MSCs from young (< 65 years of age) and old (> 65) donors in the GSE97311 dataset. C. Heatmap of DEGs in samples of MSCs from young (< 65 years of age) and old (> 65) donors in the GSE39035 dataset. D. Heatmap of DEGs in samples of MSCs from young (< 65 years of age) and old (> 65) donors in the GSE97311 dataset. Red rectangles represent a high level of expression, and blue rectangles represent a low level of expression. Group1 from young (< 65 years of age) and Group2 from old (> 65) donors. E. Venn diagram of intersecting co-expressed upregulated genes (DEGs) from the GSE39035 and GSE97311 datasets. logFC > 1, and adj.P < 0.05. F. Venn diagram of intersecting co-expressed downregulated genes (DEGs) from GSE39035 and GSE97311 datasets. logFC <-1, and adj.P < 0.05

There were 132 DEGs identified in the GSE39035 dataset, including 36 upregulated (UP) genes and 96 downregulated (DOWN) genes. Heatmap and volcano plot analyses were used to visualize these DEGs (Fig. 1C and S1A). In the GSE97311 dataset, 205 DEGs were identified, including 126 UP genes and 79 DOWN genes (Fig. 1D and S1B).

Using a Venn Diagram online, 337 intersecting genes of the two datasets were obtained, including two co-expressed UP and eight co-expressed DOWN genes (Fig. 1E and F). The details of these co-expressed genes are also shown in Table 2.

Table 2 Common upregulated and downregulated DEGs in GSE39035 and GSE97311Identification of the tissue/organ‑specific expressed genes

To identify the tissue/organ-specific expression of the DEGs in GSE39035 and GSE97311 datasets, a total of 33 tissue/organ-specific expressed genes were screened by BioGPS according to the criteria described previously (Table 3). We observed that the highest number of genes were expressed in the haematopoetic/ immune system (5/33, 15%). The second group included heart, adipocytes, and colorectal adenocarcinoma, each of which included 4 genes (4/33, 12%). The digestive system and the placenta were each represented by 3 genes (3/33, 9%) while the nervous system and genital system were represented by 2 genes (2/33, 6%). Finally, the endocrine system, respiratory system, urinary system, tongue, skeletal muscle and retina were each represented by one specifically-expressed genes (1/33 each, 3%).

Table 3 Distribution of tissue/organ-specific expressed genes identified by BioGPSEnrichment analysis

To further explore the biological function of the DEGs, the R package was used to perform GO enrichment analyses. In the biological process (BP) group, unregulated DEGs were primarily enriched in categories of extracellular structure organization, extracellular matrix organization and pattern specification process (Fig. 2A; Table 4). In contrast, downregulated DEGs were mostly enriched in embryonic organ development, embryonic organ morphogenesis and urogenital system development (Fig. 2B; Table 5). In the cellular component (CC) group, the collagen-containing extracellular matrix featured enrichment of both upregulated and downregulated genes, indicating that these genes may play different roles, respectively. Whereas the secretory granule lumen and cytoplasmic vesicle lumen showed upregulated genes and the downregulated DEGs were mostly enriched in the integral components of both the presynaptic and the postsynaptic membrane (Fig. 2A and B). In the molecular function (MF) group, upregulated genes were primarily enriched in sulfur compound binding, signaling pattern recognition receptor activity and pattern recognition receptor activity, while downregulated DEGs were mostly enriched in extracellular matrix structural constituent, DNA-binding transcription repressor activity, RNA polymerase II-specific and metallocarboxypeptidase activity (Fig. 2A and B).

Fig. 2

GO analysis results of co-expressed genes. Results of GO enrichment analysis of the DEGs that were upregulated (UP), or downregulated (DOWN). GO, Gene Ontology; MF, Molecular Function; BP, Biological Process; CC, Cellular Components; DEGs, differentially expressed genes

Table 4 Significant enriched GO analysis of UP DEGsTable 5 Significant enriched GO analysis of Down DEGs

In addition to GO enrichment analysis, functional enrichment analysis of DEGs were further performed using Reactome pathway database. The results showed that the upregulated DEGs were significantly enriched in the WNT signaling pathway (Fig. 3). PORCN (Porcupine O-Acyltransferase), one of the upregulated DEGs both in GSE39035 and GSE97311 datasets, is related to the WNT signaling pathway. Gene Ontology (GO) annotations related to PORCN include acyltransferase activity and palmitoleoyltransferase activity. It is well known that WNT pathway is critical for stem cell differentiation and development. Moreover, aberrant WNT signaling is associated with the development of numerous cancers [33]. LGK974, also shown on Fig. 3, is a PORCN-inhibitor that interferes with the WNT/β-catenin pathway ligands, and can be used for the treatment of WNT-dependent cancers [34].

Fig. 3

Reactome pathway enrichment analysis

PPI Network analysis and hub gene identification

To explore the biological characteristics of these DEGs, a PPI network was created using the STRING database. Cytoscape was used to present the PPI network modules. As shown in Fig. 4A and B, TBX15 (T-Box Transcription Factor 15), IGF1 (Insulin Like Growth Factor 1), GATA2 (GATA Binding Protein 2), PITX2 (Paired Like Homeodomain 2), SNAI1 (Snail Family Transcriptional Repressor 1), VCAN (Versican), C3 (Complement C3), COMP (Cartilage Oligomeric Matrix Protein), ADAMTS2 (ADAM Metallopeptidase With Thrombospondin Type 1 Motif 2), CBX2 (Chromobox 2), HAND2 (Heart And Neural Crest Derivatives Expressed 2), IRX3 (Iroquois Homeobox 3), SULF1 (Sulfatase 1) proteins interact with other proteins by > 4, which was the central node of the protein interaction network. These genes are the most important genes in PPI network and may play important roles in the aging of MSCs. Finally, we chose the top six genes for further analysis.

Fig. 4

A. PPI network of DEGs in samples of MSCs from young (< 65 years of age) and old (> 65) donors. PPI, protein- protein interaction; DEGs, differentially expressed genes. B. Top 13 nodes of PPI networks of DEGs in the samples of MSCs from young (< 65 years of age) and old (> 65) donors

Validation of diagnostic value of hub genes

Analysis using area under the curve (AUC), which combines both sensitivity and specificity and so can be used to describe the intrinsic effectiveness of diagnostic tests [35] was applied. To validate the diagnostic value of the top six hub genes from the above analyses, we constructed ROC curves and calculated the corresponding AUC of these gene expression levels in the GSE39035 and GSE97311 datasets. Figure 5A shows the results of GSE39035: the AUC for TBX15, IGF1, GATA2, PITX2, SNAI1, VCAN were 1.000 [95% confidence interval (Cl), 1.000–1.000; P < 0.05], 0.938 (95% Cl, 0.824-1.000; P < 0.05), 0.750 (95% Cl, 0.488-1.000; P < 0.05), 1.000 (95% Cl, 1.000–1.000; P < 0.05), 0.891 (95% Cl, 0.700-1.000; P < 0.05) and 0.984 (95% Cl, 0.941-1.000; P < 0.05). Figure 5B shows the AUC curve in the groups classified as young and old in the GSE97311 datasets. The AUC for six hub genes were 0.500 (95% Cl, 0.000–1.000; P < 0.05), 0.750 (95% Cl, 0.328-1.000; P < 0.05), 1.000 (95% Cl, 1.000–1.000; P < 0.05), 0.917 (95% Cl, 0.686-1.000; P < 0.05), 0.750 (95% Cl, 0.234-1.000; P < 0.05), 0.833 (95% Cl, 0.468-1.000; P < 0.05). These results indicate that TBX15, IGF1, GATA2, PITX2 SNAI1 and VCAN may serve as potential biomarkers for evaluating the age of MSCs donors or the quality of MSCs.

Fig. 5

ROC curve of the hub genes in samples of MSCs from young (< 65 years of age) and old (> 65) donors. Diagnostic value of top 6 hub genes with ROC curves in GSE39035 dataset (A) and in GSE97311 dataset (B). AUC area under the ROC curve

Comparison of hub genes expressions between normal and tumor tissues

To evaluated hub gene (TBX15, IGF1, GATA2, SNAI1 and VCAN) expression at the protein level, the IHC images were download from the HPA database (http://www.proteinatlas.org/). As shown in Fig. 6A-E, normal ovarian, lung, prostate, breast and liver had negative or moderate IHC staining, while tumor tissues had strong staining. Furthermore, we also noted that these genes also showed high expression in other tumor tissues, such as colorectal and pancreatic (not shown). Some hub genes (TBX15 and VCAN) are not well known, but have been reported to be highly expressed in many tumors in recent years, such as gliomas [36], ovarian cancer [37]. The other hub genes, IGF1, GATA2 and SNAI1, have been reported to be widely expressed and involved in the development of many tumors [38,39,40]. Therefore, the results indicate that these hub genes may be involved in the development of several tumors, which is consistent with the results of our previous bioinformatics analysis.

Validation of hub genes expression in senescent MSCs

To further assess whether hub genes can be used as diagnostic markers for senescence in MSCs donors, we cultured human bone marrow-derived MSCs in vitro and induced cellular senescence with different concentrations of H2O2. As shown in Fig. 7A, the viability of MSCs was significantly reduced in a concentration-dependent manner after H2O2 treatment. We chose 100 µM H2O2 to treat MSCs for 2 h and extracted total mRNA to detect changes in the expressions of senescent markers p16 and p21. As shown in Fig. 7B, we observed that expressions of both senescence markers, p16 and p21, was significantly upregulated in MSCs after H2O2 treatment. Moreover, by comparing the relative telomere lengths in MSCs and MSC-H2O2 groups, we further confirmed that the relative telomere length was significantly shortened in MSCs after H2O2 treatment (Fig. 7C). These results indicated that MSCs entered senescence state following induction by H2O2. Subsequently, we used qPCR analysis to detect the expression of the six hub genes in the H2O2-induced MSCs at the transcriptional level (Fig. 7D). Compared with the control MSCs group, the expression of all six hub genes, including TBX15 (p < 0.01), IGF1 (p < 0.01), GATA2 (p < 0.01), PITX2 (p < 0.05), SNAI1 (p < 0.01) and VCAN (p < 0.01) (p < 0.01), was significantly upregulated in senescent MSCs. These results suggest that the expressions of these hub genes are elevated in senescent MSCs, which is consistent with the foregoing results in this paper. It indicates that the expression of some potentially tumor-associated genes may be upregulated in senescent MSCs donors compared to young donors, and thus there may be a higher safety risk of using MSCs from old donors to treat diseases.

Fig. 6

Immunohistochemistry images of hub genes in normal and tumor tissues extracted from the HPA (http://www.proteinatlas.org/). TBX15 protein expression was significantly higher in ovarian cancer tissue than normal tissue. IGF1 protein expression was significantly higher in lung cancer tissue than normal tissue. GATA2 protein expression was significantly higher in prostate cancer tissue than normal tissue. SNAI1 protein expression was significantly higher in breast cancer tissue than normal tissue. VCAN protein expression was significantly higher in liver cancer tissue than normal tissue