A comprehensive analysis of age-related metabolomics and transcriptomics reveals metabolic alterations in rat bone marrow mesenchymal stem cells

Multivariate statistical analysis

Based on the data of non-targeted metabolomics, we established a multivariate statistical analysis model to reveal the changes of metabolites in old BMSCs and young BMSCs. Intensities were corrected for signal drift and batch effect by fitting a locally quadratic (loess) regression model to the median intensity of pooled QC samples. After correction, the median area of all pooled QC samples was now the same. Metabolites with a coefficient of variation (CV) in QC samples >25% were then filtered out, due to their unstable quantifiability. In the end, a total of 1886 annotated metabolites were obtained, including positive ion modes and negative ion modes, as well as some details information such as m/z, retention time and peak intensity.

Unsupervised principal component analysis (PCA) of the metabolic data was used to observe grouping trends, prominent outliers, and clustering between display groups to provide an initial assessment of metabolic disturbances. PCA score plots showed the distribution between the old BMSCs and young BMSCs as well as QC samples in two dimensions (Figure 1A). According to Supplementary Figure 1, the QC samples were tightly clustered within 95% confidence interval, which proved that the method was steady and the instrument had good repeatability.

Figure 1. Results of multivariate statistical analysis of metabolomics. (A) Multivariate statistical analysis of metabolites-PCA. PC1 and PC2 represent the degree of interpretation of the model with the first and second principal component ranking in the principal component analysis. (B) Multivariate statistical analysis of metabolites-OPLS-DA. The separation of the two classes of observations occurs in the horizontal (t1) direction; The vertical (t1o) direction indicates intra-class variability. (C) Permutations Plot for OPLS-DA Model. For the selected Y variable, the figure shows the R2 and Q2 values of the original model (far right) on the vertical axis; The horizontal axis shows the correlation between the substituted Y vector and the original Y vector of the selected Y. (D) Volcano plots of metabolites. The red dots represent up-regulated metabolites and the blue dots represent down-regulated metabolites, and the grey dots represent that the metabolites were not significantly changed.

We further developed a supervised OPLS-DA methodological model (R2X = 0.946, R2Y = 1, Q2 = 0.991) to better account for changes in metabolites between the old and young cells (Figure 1B). A permutation test of 200 random numbers was used to verify the OPLS-DA model and the result showed that all blue Q2-values to the left were lower than the original points to the right (Figure 1C), which proved the OPLS-DA model was effective and stable. Judging by the OPLS-DA model, there was a distinct difference between the two groups, indicating an apparent disturbance in the senescent BMSCs.

Volcano plots were constructed to determine the difference in metabolites between the old and young cell groups (Figure 1D). In the diagram, each dot represents a metabolite. As shown in the Figure 1D, a total of 130 metabolites showed visible alterations, 92 metabolites were up-regulated (red dots), 38 metabolites were down-regulated (blue dots), and most of the other metabolites were not significantly changed (grey dots). In the volcanic plots, the further away a point was from the origin, the more it contributes to the distinction between the two groups. Thus, the distant points from the origin were candidate biomarkers for senescent BMSCs.

Identification of candidate biomarkers

On the basis of the accuracy of MS data, the exact mass of each feature was submitted to ChemSpider with 4 databases selected (BioCyc; Human Metabolome Database; KEGG; LipidMAPS). According to the predicted value of VIP in the OPLS-DA model and the return p-value of the t-test, 23 metabolites with significant changes were screened as candidate biomarkers (VIP > 1 and p-value <0.05), which were displayed using a heatmap (Figure 2). Details of 23 differential metabolites are shown in Table 1. As shown in the heatmap, the differential biomarkers can distinguish the old BMSCs group from the young group, indicating that the biomarkers we have obtained were reliable.

Figure 2. Heatmap of metabolites with significant differences. Each row represents a sample and each column represents a metabolite. Red is the high expression level, blue is the low expression level.

Metabolic pathway analysis

After identifying the candidate biomarkers, we enriched 15 disordered metabolic pathways by using MetaboAnalyst 4.0: glycerophospholipid metabolism, linoleic acid metabolism, biosynthesis of unsaturated fatty acids, arachidonic acid metabolism, caffeine metabolism, taurine and hypotaurine metabolism, phenylalanine metabolism, alpha-linolenic acid metabolism, glycosylphosphatidylinositol (GPI)-anchor biosynthesis, fructose and mannose metabolism, sphingolipid metabolism, lysine degradation, cysteine and methionine metabolism, glycine, serine and threonine metabolism and pyrimidine metabolism (Figure 3). Details of these disordered metabolic pathways are shown in Table 2.

Figure 3. The enriched pathway analysis of metabolites. There are five metabolic pathways with significant changes: (A) Glycerophospholipid metabolism; (B) Linoleic acid metabolism; (C) Biosynthesis of unsaturated fatty acids; (D) Arachidonic acid metabolism; (E) Caffeine metabolism.

Transcriptomics profiling of the BMSCs

The results of the transcriptional group showed that there was a remarkable difference in gene transcription levels between the old BMSC group and the young BMSC group. A total of 16205 RNAs were detected, and 14574 meaningful RNAs (FPKM > 0.5) were screened. We further confirmed the screening conditions of |log2(foldchange)|> 1.5 and p-value < 0.05. Based on this condition, 590 RNAs with significant differences were considered to be DEGs. Compared with the young group, the expression of 311 RNAs were down-regulated and 279 RNAs were up-regulated in the old BMSCs group. The volcanic plots visually showed the variation trend of these RNAs (Figure 4A). We proved further that the 590 DEGs identified in the form of a heatmap (Figure 4B), and we found that DEGs could separate the old BMSC group from the young group, indicating that the DEGs we obtained were reliable.

Figure 4. Results of statistical analysis of transcriptome differences. (A) Volcano plots of mRNAs. The volcanic map is drawn with -log10 (p-value) as the vertical axis and log2 (Fold Change) as the horizontal axis, which shows the changes of mRNA expression in the old BMSC group and the young BMSC group. The red dots represent up-regulated mRNAs and the blue dots represent down-regulated mRNAs, and the grey dots represent that the mRNAs were not significantly changed. (B) Heatmap of DEGs. Each row represents a sample and each column represents a mRNA. Red is the high expression level, blue is the low expression level.

GO and KEGG enrichment analyses of DEGs

GO and KEGG enrichment analysis are based on different angles to better explain and understand the biological significance of DEGs. GO enrichment analysis focuses on the functional enrichment of DEGs. From the results of GO enrichment analysis, we found that up-regulated DEGs were enriched in 270 significant items (Supplementary Table 1), and down-regulated DEGs were enriched in 73 significant entries (Supplementary Table 2). Among the GO items into which up-regulated DEGs were enriched, we selected the top 10 entries of BP (biological process), CC (cellular component) and MF (Molecular Function) to draw bar charts (Figure 5A). Among the GO entries enriched by down-regulated DEGs, top 15 of BP, all CC and all MF entries were drawn as bar charts (Figure 5B).

Figure 5. Gene ontology and KEGG pathway enrichment results. (A) GO enrichment analysis of up-regulated DEGs. (B) GO enrichment analysis of down-regulated DEGs. (C) KEGG enrichment analysis of up-regulated DEGs. (D) KEGG enrichment analysis of down-regulated DEGs. Abbreviations: BP: Biological Process; CC: Cellular Component; MF: Molecular Function.

KEGG enrichment analysis focuses on gene pathway analysis. From the results of pathway enrichment, we found that 279 up-regulated genes enriched in 170 KEGG pathways (Supplementary Table 3) and 311 down-regulated genes enriched in 240 KEGG pathways (Supplementary Table 4). It was enriched in 38 significantly up-regulated pathways (Figure 5C) and 38 down-regulated pathways (Figure 5D).

Protein-protein interaction network

KEGG pathways related to lipid metabolism were selected, and protein-protein interaction (PPI) network was further constructed between the des enriched in these pathways (Figure 6). Scd, Mboat2, Agpat5, Dgat2, Fads2, Lpin1, Acaa2, Lpcat3, Pcyt2 and Pla2g4a were strongly correlated with other proteins (connected with >3 proteins).

Figure 6. PPI network construction. Each node represents one differentially expressed protein. Each edge represents regulation. Red, up-regulated expressed protein; blue, down-regulated expressed protein.

Integrative transcriptomic and metabolomic molecular profiling analysis

Furthermore, we analyze the changes of metabolic pathways in the process of cellular senescence by integrating metabolomic and transcriptome data, so as to better understand the mechanisms of BMSCs senescence from a metabolic point of view. 52 KEGG metabolic pathways (Figure 7A) were enriched by integrating data, of which 5 metabolic pathways (p-value < 0.05 and Impact > 0.1) had visible changes in the process of cellular senescence, and 29 metabolic pathways (p-value > 0.05 but Impact > 0.1) were potentially associated with cellular senescence. In order to understand the metabolic response of senescent BMSCs at the transcriptome level, correlation analysis heat map of the connection between differential metabolites and DEG was established (Figure 7B). The DEGs involved in lipid metabolism exhibited the most dramatic changes in BMSC senescence. The differential metabolites in lipid metabolism were extremely correlated and statistically significant with DEGS.

Figure 7. Multiple omics combined analysis results. (A) The heatmap of differential metabolites and genes. (B) The enriched pathway analysis of integrative data. There were five metabolic pathways with significant changes: (a) Glycerophospholipid metabolism; (b) Taurine and hypotaurine metabolism; (c) Glycerolipid metabolism; (d) Drug metabolism - cytochrome P450; (e) Drug metabolism - other enzymes.

Age-associated changes in BMSCs and validation of mRNA expression

In order to verify the results of data analysis, we first obtained young and old BMSCs. With the aging of BMSCs, the cells gradually changed from long spindle to flat, lost vitality, and there were clear particles in the cytoplasm (Figure 8A). The analysis results indicated that compared with young BMSCs, the cell aspect ratios in old BMSCs decreased (Figure 8B) and the cell areas increased (Figure 8C). And the number and ratio of SA-β-gal positive cells in old BMSCs were obviously elevated by SA-β-gal staining (Figure 8D, 8E). Then the classical aging evaluation marker pl6^INK4A mRNA was up-regulated in old BMSCs (Figure 8F). These data above displayed that old BMSCs presented senescent alterations. Further, RT-qPCR results demonstrated that consistent with the trend of high-throughput sequencing analysis, the expressions of Scd, Scd2, Dgat2, Fads2, Lpin1, Gpat3, Acaa2 and Pla2g4a in old BMSCs were down-regulated compared with the young group, and Scd2 decreased most significantly, while the expressions of Lpcat3 and Pcyt2 were up-regulated, which was opposite to the results of high-throughput sequencing (Figure 8G).

Figure 8. Detection of age-related changes in BMSCs and verification of gene expression. (A) Morphological alterations were observed under a phase-contrast microscope (scale bar = 100 μm). (B) The cell areas were obviously increased in the old BMSCs, and (C) the cell aspect ratios were clearly decreased. (D–E) SA-β-gal staining (scale bar = 100 μm). Senescent cells are stained blue. Compared with young cells, the proportion of positive cells in old BMSCs was significantly elevated. (F) RT-qPCR analyses of mRNA expression of the age-related factor pl6^INK4A. BMSCs obtained from aged rats expressed elevated levels of pl6^INK4A. (G)RT-qPCR detected the expression of the DEGs. Data indicate the mean ± SD, n = 3. ^*P < 0.05, ^**P < 0.01, ^***P < 0.001 vs. Young (Y).

Scd2 over-expression ameliorates the senescence of BMSCs

Since Scd2 gene was identified to be most significantly reduced in old BMSCs, further study was investigated whether BMSC senescence could be influenced by enforcing Scd2 expression. For this purpose, senescent BMSCs were transduced with lentivirus expressing Scd2 (LV-Scd2) and the lentiviral vector (LV-Vector). Fluorescence microscopy showed that Scd2 was successfully over-expressed in old BMSCs (Figure 9A). Scd2 expression at protein level was obviously up-regulated in the LV-Scd2 group compared with that in the LV-Vector group, which was confirmed by Western blot (Figure 9B). Furthermore, SA-β-gal staining and quantitative analysis demonstrated that SA-β-gal activity in Scd2-suffficient BMSCs was dramatically diminished in the LV-Scd2 group (Figure 9C). And P16^INK4a mRNA expression was markedly down-regulated after Scd2 repletion (Figure 9D). Accordingly, the above data manifested that Scd2 replenishment can attenuate BMSC senescence.

Figure 9. Scd2 over-expression alleviates senescence-associated variations in BMSCs. (A) EGFP expression under the fluorescence microscope (scale bar = 200 μm). (B) Determination of Scd2 protein expression levels by Western blot to demonstrate the transduction efficiency. (C) SA-β-gal staining (scale bar = 100 μm) and quantification of β-gal positive cells. (D) Gene expression of the senescence-related factor p16^INK4a. Data indicate the mean ± SD, n = 3. ^*P < 0.05, ^**P < 0.01, vs. LV-Vector.

Materials and reagents

Healthy male Wistar rats, SPF grade, purchased from Animal Experimental Center, Basic Medical College of Jilin University, P.R. China. Liquid chromatography-mass spectrometry (LC-MS) grade ammonium acetate, ammonium hydroxide and MeOH were purchased from Sigma-Aldrich (USA). Acetonitrile was purchased from J.T.Baker (USA). KAPA Stranded RNA-Seq Library Prep Kit and TruSeq SR Cluster Kit v3-cBot-HS were purchased from Illumina (USA). NEBNext^® Poly(A) mRNA Magnetic Isolation Module was purchased from NEB (USA). RiboZero Magnetic Gold Kit (Human/ Mouse/Rat) was purchased from Epicentre (USA).

BMSC isolation and culture

Primary BMSCs, which were isolated from the bone marrow of young (1–2 months, n = 3) and old (15–18 months, n = 3) male Wistar rats (Permit Number: SYXK 2018-0001) by whole bone marrow attachment method under aseptic conditions, were plated into 10-cm culture dishes in complete medium containing 89% Dulbecco's Modified Eagle Medium with nutrient mixture F-12(DMEM-F12) (HyClone, USA) supplemented with 10% FBS (Gibco, USA) and 1% penicillin streptomycin (HyClone, USA) [43]. The young and old BMSCs were incubated at 37°C with 5% CO2 for 24 h in a humidified atmosphere, the medium was replaced by half amount. After that, the liquid was changed every 2–3 days. When reaching 80% confluence, 0.25% trypsin-EDTA (Gibco, USA) was added to digest the cells and the subculture was continued at a ratio of 1:3 BMSCs at passage 3, young BMSCs were used as the control group and old BMSCs as the experimental group for subsequent experiments.

Sample preparation and intracellular metabolites extraction

The culture media of young and old cells were removed, and the cells were washed with pre-cooled PBS. Then a small amount of PBS was added to the petri dish, and the cells were carefully scraped and collected into the centrifuge tube. After that, centrifugation was carried out at room temperature, 1000 rpm, for 5 min. The upper liquid was removed as far as possible without residue and the lower cell was precipitated. H2O was added into the cells at 4°C and ultrasound was performed for 10 min. MeOH and CAN were added in a 1:1 ratio to extract the spare sample 200 mL, then eddy the sample for 30s, ultrasound for 10 min. To remove the protein, the sample was incubated at −20°C for 1 h, and then centrifuged at 4°C at 20,000 rpm for 15 min. The supernatant produced by centrifugation was discarded in a vacuum concentrator and evaporated to dry. ACN and H2O were added to the precipitation at the ratio of 1:1, and then the mixture was mixed, eddy for 30s, and ultrasonic for 10 min. The insoluble fragments were removed by centrifugation at 20000 rpm at 4°C for 15 min. The supernatants were transferred to HPLC vials and stored at −80°C prior to LC/MS analysis [44]. QC samples were prepared by pooling 10 μL from each sample. The extraction of QC samples was the same as sample preparation.

LC-MS measurements

Samples were separated on an amide column, using mobile phase A consisting of water mixed with 25 mM ammonium acetate and 25 mM ammonium hydroxide and mobile phase B ACN. The injection volume was 4 μL and flow rate was 0.4 ml/min. The generic HPLC gradient was listed in Supplementary Table 5. Then, MS analysis was carried out on the Q-Exactive MS/MS in both positive and negative ion modes.

The relevant tuning parameters for the probe were set as listed: aux gas heater temperature, 400°C; sheath gas, 40; auxiliary gas, 13; spray voltage, 3.5 kV for positive mode and negative mode. The capillary temperature was at 350°C, and S-lens at 55. A DDA method was built as follows: Full scan range: 60 to 900 (m/z); resolution for MS1 and ddMS2: 70,000 and 17,500 respectively; maximum injection time for MS1 and ddMS2: 100 ms and 45 ms; automatic gain control (AGC) for MS1 and ddMS2: 3e6 and 2e5; isolation window: 1.6 m/z; normalized collision energies (NCE): 10, 17, 25 or 30, 40, 50. Build a full scan method as follows: Full scan range: 60 to 900 (m/z); resolution: 140,000; maximum injection time: 100 ms; automatic gain control (AGC): 3e6 ions. TIC plots for sequentially selected QC samples were displayed as bellow (Supplementary Figure 1).

Metabolomics analysis

Raw files were submitted to Thermo Compound Discover 2.1, (CD), and processed with Untargeted Metabolomics workflow with minor modification to find and identify the differences between samples. The quantification and annotation results from CD were analyzed with R script. In the analysis of multivariate data, the normalized spectral data were imported into SIMCA-P software (Version 12.0, Umetrics AB, Umea, Sweden), and the importance of low-level metabolites was increased by Pareto scaling, while the noise was not significantly amplified. Unsupervised principal component analysis (PCA) was used to observe grouping trends, highlight outliers, and show clustering between groups [45].

On this basis, the relative integral of significant metabolites determined in the OPLS-DA model was used to quantitatively analyze the metabolite levels. Variations were statistically calculated by utilizing one-way analysis of variance (ANOVA) followed by Tukey’s multiple comparison test with the Bonferroni correction using the R software (Version x64 4.0.2) and Rstudio (Version x64 3.6.2). Differential metabolites were identified with statistical significances of ^*P < 0.05.

RNA isolation and sequencing

When the growth of young and old cells converged and fused to 80%–90%, the culture medium was abandoned and cells were washed with PBS once or twice; 1–2 ml Trizol was added to the adherent area of 10–15 cm² cells. After repeated batting for several times, the visible cell layer was completely dissolved and transferred to the RNA-free EP tube. Agarose gel electrophoresis was used to detect the integrity of the total RNA in the sample, and quantitative analysis was carried out using Nano DropNd-1000. After that, RNA was inspected. Fragmented mRNA enrichment or rRNA removal; The first strand of cDNA was generated by reverse transcription, and the double strand of cDNA was synthesized by adding dUTP. Various polymerases were added to repair the ends of the double-stranded cDNA and A was added, followed by ligase Illumina specific connector; The final RNA sequencing library was obtained by PCR amplification and purification. Agilent 2100 BioAnalyzer was used for quality control, and the final quantification of the library was carried out using qPCR method.

The sequencing library of the mixed samples was deformed by Na OH to generate single stranded DNA, which was diluted into 8p M concentration and amplified in situ on Tru Seq SR Cluster Kit V3-C BOT-HS. The end of the generated fragment was sequenced for 150 cycles.

Transcriptomics analysis

The software StringTie was used to compare the results to the known transcriptome, calculate the transcriptional abundance. Differentially expressed mRNAs were detected by using the Ballgown [46]. The unit of expression amount was expressed by FPKM, and the threshold value of gene or transcript expression amount was FPKM greater than or equal to 0.5. Genes or transcripts whose mean value of FPKM exceeded 0.5 in each group are considered to be expressed in the group for statistical analysis. In addition, R packages were clusterProfiler (3.16.1) and ggplot2(3.3.2) to perform Gene ontology (GO) and KEGG pathway enrichment to describe the properties of genes and gene products of various organisms, and to determine the different biological pathways involved in mRNAs [47]. The DEGs were analyzed using STRING (http://string-db.org) database and a protein-protein interaction (PPI) network was visualized using Cytoscape (version 3.6.2).

Integration of metabolomic and transcriptomic analysis

Metabolome and transcriptome data were visualized and clarified by using the network analysis module of online analysis software Metaboanalyst 4.0 (https://www.metaboanalyst.ca/). In addition, the metabolite differences and abundant pathways from transcriptome data were able to be clearly observed and identified.

Analysis of the expression levels of mRNAs using RT-qPCR

The reliability of high-throughput RNA sequencing data was further confirmed. The expression levels of the top eight up-regulated mRNAs were measured using real-time quantitative polymerase chain reaction (RT-qPCR) with TransStart Top Green qPCR SuperMix (TRANS, China) in a 7300 Real-Time PCR System (ABI, Vernon, CA, United States). The primer sequences used are shown in Supplementary Table 6. SYBR Green Supermix was used for RT-qPCR. The expression was determined using the threshold crossing point (Ct) as calculated by ΔΔCt.

Detection of senescence-associated β-galactosidase (SA-β-gal) activity

The young and old BMSCs in good growth condition were detected with the Senescence β-Galactosidase Staining Kit (Beyotime Biotechnology, China). After washing with phosphate buffered saline (PBS) and fixing, the cells were stained overnight at 37°C in the staining solution containing X-Gal. The staining was observed under inverted fluorescence microscope, and the positive cell rate (number of blue stained cells/total cells × 100%) in different areas of the culture dish was calculated.

Lentiviral transduction of BMSCs

The optimal infection conditions of senescent BMSCs were determined according to the cell density and virus titer. The cells were inoculated at a density of 1.5×10⁵ cells/well in a 6-well plate. When about 50% confluence reached, cell complete culture medium, Scd2 over-expression lentivirus and virus infection enhancement solution were added in a certain ratio. After 12 h incubation, the original medium was discarded and replaced with new cell complete culture medium. Then EGFP expression was monitored under a fluorescence microscope and transduction efficiency was determined by western blot after subsequent 48–72 h.

Western blot analysis

After protein extraction using RIPA lysis buffer, the total protein content was determined using the BCA protein assay kit (Beyotime). Then, 30 μg of protein extracted from each sample was resolved by 10% SDS-PAGE gels and transferred by electroblotting onto PVDF membranes (Millipore, Billerica, CA, USA). The blotted membranes were incubated with 5% skim milk for 1–2 h at RT and then detected overnight at 4°C with anti-Scd2 (1:500 dilution, Santa Cruz, USA) and anti-β-actin (1:2000 dilution, Proteintech, USA) diluted in TBST. Incubation with horseradish peroxidase-conjugated anti-rabbit IgG secondary antibody (1:2000 dilution, Beyotime Biotechnology, China) was performed for 1–2 h followed by washing with TBST. Finally, protein blots were visualized using an enhanced Electro-Chemi-Luminescence detection system (Amersham Biosciences, Inc., Piscataway, NJ, USA). β-actin was used as an internal standard.