Ten patients aged 40–60 years and age-matched healthy controls, were recruited for this study. Stool and blood samples were collected on-site at Weifang People’s Hospital. The stool samples were initially stored in ice boxes before being transferred into centrifuge tubes. The collection date was documented, and the samples were and stored at −80 °C until further processing. Peripheral blood mononuclear cells (PBMCs) were isolated from the blood samples for RNA extraction and subsequently stored at −80 °C following reverse transcription. The study was approved by the Ethics Committee of Weifang People’s Hospital (approval number: KYLL20241213-6), and informed written consent was obtained from all participants.

The human umbilical cord was rinsed with phosphate-buffered saline (PBS) and subsequently sectioned into pieces measuring 1 mm³. These tissue blocks were then incubated in DMEM/F12 supplemented with 10% fetal bovine serum (FBS), penicillin (100 U/mL), and streptomycin (100 mg/mL). After approximately 7 days, MSCs emerged from the tissue blocks and proliferated in the surrounding area. Following an additional 3–5 days, the adherent hUC-MSCs were detached using 0.25% trypsin–EDTA. The hUC-MSCs were cultured in 3D FloTrix® mesenchymal stem cell medium (CytoNiche Biotech, Beijing, China) containing penicillin (100 U/mL) and streptomycin (100 mg/mL). The cells were expanded over three to five passages for subsequent use. Flow cytometry analysis was conducted to identify the presence of positive markers CD44 (Biolegend, 338,805), CD73 (Biolegend, 344,003), CD90 (Biolegend, 328,109), and CD105 (Biolegend, 800,505), as well as the absence of negative markers CD34 (Biolegend, 343,603) and CD45 (Biolegend, 304,011) on the hUC-MSCs. The hUC-MSCs were then induced to differentiate into adipocytes (Procell, PD-019), chondrocytes (Procell, PD-018), and osteocytes (Procell, PD-017) using specific stem cell differentiation media (Procell, Wuhan, China), respectively.

The animal experiments were approved by the Animal Ethics Committee of Shandong Second Medical University (approval number: 2022SDL339). Female C57BL/6 mice, aged 6 weeks, were purchased from the Animal Experimental Center of Shandong Second Medical University, and all experiments involving mice were conducted under specific pathogen-free (SPF) conditions. Following a one-week acclimatization period, mice in the model group received subcutaneous injections of BLM (100 μL, 0.5 mg/mL) on the back daily for 4 weeks. The healthy control group (CON) was administered an equivalent volume of PBS. Subsequently, the model mouse was randomly allocated into two groups: the hUC-MSC treatment group, receiving 5 × 10⁵ hUC-MSCs per mouse (MSC), and the disease control group, receiving 100 μL PBS per mouse (SSC). Caudal intravenous therapy was administered during the fifth and seventh weeks. Upon completion of the treatment period (the eighth week), all mice were sacrificed and lung tissues and stool samples were collected for further analysis.

A portion of the lung tissues from the mice was fixed in 4% paraformaldehyde, while the remaining tissues were preserved at −80 °C. Subsequent to dehydration and paraffin embedding, the lung tissues were sectioned using a rotary microtome to a thickness of 5 μm (He et al., 2024). The sections were stained with hematoxylin and eosin (H&E) as well as Masson’s trichrome stain. The expression of α-SMA was assessed by a primary antibody against α-SMA (dilution 1:500, CST, 19245 T) incubated overnight at 4 °C, followed by a horseradish peroxidase (HRP)-labeled goat anti-rabbit secondary antibody (dilution 1:200, Servicebio, GB23303) for 50 min at 25 °C. Similarly, COL1A1 expression was assessed using a primary antibody COL1A1 (dilution 1:100, Affinity, AF7001) overnight at 4 °C and HRP-labeled goat anti-rabbit secondary antibody (dilution 1:200, Servicebio, GB23303) for 50 min at 25 °C. The samples were digitally scanned using the application system EasyScan 6 system (Motic, Xiamen, China) to obtain high-resolution images.

Human blood PBMCs and mice lung tissues were homogenized using TriQuick Reagent (Solarbio, Beijing, China) to facilitate the extraction of total RNA, employing the FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, Nanjing, China). Complementary DNA (cDNA) was synthesized from the total RNA samples by reverse transcription via HiScript II Q-RT SuperMix for quantitative PCR (qPCR) (Vazyme, Nanjing, China). The qPCR assays were performed using Taq Pro Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) on the Applied Biosystems QuantStudio 5 real-time PCR system (Thermofisher, USA). The sequences of the primers used in the assays are detailed in Supplementary Table S1.

The total genomic DNA of stool samples was extracted using the Magnetic Soil and Stool DNA Kit (TianGen, Beijing, China) following the manufacturer’s instructions. The quality and quantity of DNA were verified with NanoDrop and agarose gel. The genomic DNA was randomly sheared into short fragments, and sequencing libraries were generated. The obtained fragments were end-repaired, A-tailed, and further ligated using an Illumina adapter. The fragments with adapters were amplified by PCR, and the size was selected and purified. The library was checked with Qubit, qPCR was used for quantification, and a bioanalyzer was used for size distribution detection. Quantified libraries were pooled and sequenced on the Illumina platform according to the required effective library concentration and data amount. Readfq¹ was used to preprocess raw data from the Illumina PE150 sequencing platform (Novogene, Beijing, China) to obtain clean data for subsequent analysis. MEGAHIT software was used for assembly analysis of clean data, as previously described (Karlsson et al., 2013; Nielsen et al., 2014). DIAMOND software² (Buchfink et al., 2015) was used for the alignment of Unigene sequences with those of bacteria, fungi, archaea, and viruses extracted from NCBI’s NR database.³ LCA algorithm (applied to systematic taxonomy of MEGAN software⁴ was adopted to determine the species annotation information of the sequence; Huson et al., 2011). Based on the results of LCA annotation, the abundance of each sample at each taxonomy (kingdom, phylum, class, order, family, genus, or species) was acquired. Alpha diversity indices included several observed species, Shannon, Chao1, and ACE. Beta diversity analysis was used to evaluate differences among different groups. Group differences were statistically tested using Analysis of Similarities (ANOSIM) via the R vegan package. Linear discriminant analysis effect size (LEfSe) was performed to investigate potential indicator species using the Kruskal–Wallis (KW) sum-rank test within gut microbes enriched explicitly in a specific group (p < 0.05, LDA ≥ 4) (Xiang et al., 2020; Zhao et al., 2019; Segata et al., 2011). Raw sequences data of microbial metagenomics are available in the NCBI Sequence Read Archive (SRA) database under accession number PRJN1124507.⁵

Stool samples from humans (n = 10 per group) and mice (n = 6 per group) were subjected to non-targeted metabolomics analysis. Each sample (100 mg) was individually using liquid nitrogen, and the homogenate was resuspended with prechilled 80% methanol and vortexed well. The samples were then incubated on ice for 5 min and centrifuged at 15,000 × g at 4 °C for 20 min. A portion of the supernatant was diluted to achieve a final concentration containing of 53% methanol using Liquid Chromatograph Mass Spectrometer (LC–MS) grade water, and then centrifuged at 15000 × g at 4 °C for 20 min. Finally, the supernatant was injected into a liquid chromatography-tandem mass spectrometer (LC–MS/MS) system for analysis (Want et al., 2013). Ultra High Performance Liquid Chromatography–Tandem Mass Spectrometry (UHPLC–MS/MS) analyses were performed using a Vanquish UHPLC system (ThermoFisher, Germany) coupled with an Orbitrap Q ExactiveTM HF-X mass spectrometer (Thermo Fisher, Germany) at Novogene Co., Ltd. (Beijing, China). Principal components analysis (PCA) and partial least squares discriminant analysis (PLS–DA) were performed at metaX (Wen et al., 2017). The metabolites with VIP > 1, p < 0.05, and fold change (FC) ≥ 1.5 or ≤ 0.667 were considered differential metabolites. The functions of these metabolites and metabolic pathways were studied using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

The lung tissues from six mice per group were processed for transcriptomic analysis. The integrity of the RNA was evaluated utilizing the RNA Nano 6,000 Assay Kit of the Bioanalyzer 2,100 system (Agilent Technologies, CA, USA). The library preparation workflow included mRNA enrichment, cDNA generation, and end repair to generate blunt ends, adaptor ligation, and PCR amplification. Subsequently, PCR products were purified, and the library quality was assessed using the Agilent Bioanalyzer 2,100 system. Sequencing was performed on an Illumina Novaseq platform (Novogene, Beijing, China), and 150-bp paired-end reads were generated. Clean data (clean reads) were obtained by removing reads containing the adapter, reads containing ploy-N, and low-quality reads from raw data (raw reads). All the downstream analyses were based on these clean, high-quality data. Differentially expressed genes (DEGs) were identified based on a fold change (FC) ≥ 1.5 and p < 0.05. Gene Ontology (GO) and KEGG pathway enrichment analysis of DEGs was implemented using the clusterProfiler R package, in which gene length bias was corrected. GO terms with corrected p < 0.05 were considered significantly enriched by DEGs.

Gut microbiota, metabolomics, and RNA-seq analysis were analyzed and visualized using R software (version 4.3.3, http://www.r-project.org). Principal coordinate analysis (PCoA; packages: vegan, ggplot2) based on the Bray–Curtis distance matrix was performed to compare the differences in microbial community structure. Permutational Multivariate Analysis of Variance (PERMANOVA) (package: vegan, adonis function) based on the Bray–Curtis distance matrix was employed to assess the significant differences in the gut microbiota (Ni et al., 2022; Shuai et al., 2022). Mantel tests (package: corrplot, ggplot2) and Spearman correlation analysis (package: vegan, ggcor, ggplot2) were performed to observe the correlation between the extracellular matrix, collagen-related genes, differential taxa, and metabolites. All statistical differences in this study were considered significant at p < 0.05.

The diversity and structural composition of microbial communities were analyzed based on the results of microbial metagenome sequencing. Compared with healthy controls (HC), patients with SSc (SSC) exhibited a reduction in the α-diversity indices (observed species, Shannon, Chao1, ACE) of the gut microbial community (Figure 1A). Principal coordinate analysis (PCoA) at the species level demonstrated a significant separation between gut microbial communities of the HC and SSC groups (PERMANOVA: p = 0.0134, Figure 1B). Consistent findings were observed through ANOSIM and cluster analysis (p = 0.046, Supplementary Figures S1A,B), suggesting potential dysbiosis in the gut microbial community of SSc patients. Bacillota (65.01%), Bacteroidota (13.69%), and Actinomycetota (6.53%) were the main microbial phylum. Compared with the HC group, the relative abundance of Bacteroidota, Phocaeicola, Eubacterium, Phocaeicola coprocola, Faecalibacterium prausnitzii, and Parabacteroides in the SSC group were significantly decreased (Figures 1C–E; Supplementary Figures S1C,D). By contrast, the relative abundances of Sphingomonadales, Myxococcales, Sphingomonadaceae, and Klebsiella were significantly increased (Supplementary Figure S1D). Linear discriminant analysis effect size (LEfSe) results showed that Parabacteroides, Bacteroidales, Bacteroidaceae, and Phocaeicola plebeius were significantly enriched in the HC group, while Enterobacteriaceae was significantly enriched in the SSC group (Figure 1F). ANOSIM based on functional abundance at the KEGG KO level revealed a significant differentiation in gut microbial function between the HC and SSC groups (R = 0.092, p = 0.046, Figure 1G). Compared to the HC group, carbohydrate metabolism, cellular community, amino acid metabolism, nucleotide metabolism, lipid metabolism, metabolism of cofactors and vitamins were predicted to be associated with alterations of gut microbiota in SSc patients (Supplementary Figure S2). Collectively, these findings suggest that gut microbiota dysbiosis is present in SSc, potentially influencing the progression of the disease through the aforementioned pathways.

qPCR analysis revealed a significant upregulation in the expression levels of TNF-α and IL-6 in the blood of patients with SSc (p < 0.05, Supplementary Figure S3A), suggesting an enhancement of systemic inflammation in these patients. Correlation analysis further demonstrated that the expression levels of inflammatory cytokines TNF-α and IL-6 were significantly associated with the relative abundance of microbial taxa, including Sphingomonadales, Myxococcales, Sphingomonadaceae, and Klebsiella (p < 0.05, Supplementary Figure S3B). Based on these findings, we propose the hypothesis that dysbiosis of the gut microbiota may play a role in the immune and inflammatory responses observed in SSc and could be implicated in its pathogenesis.

To understand the differences in gut microbial metabolic function, we performed an untargeted metabolomics analysis of the collected stool samples. Principal Component Analysis (PCA) and Partial Least Squares Discriminant Analysis (PLS-DA) revealed a significant separation in the gut microbiota metabolite composition between the HC and SSC groups (Figures 2A–E,G), suggesting altered patterns in the composition of gut microbial metabolites in SSc patients. Volcano plot analysis was performed to identify differentially expressed gut microbiota metabolites (p < 0.05, FC = 1.5). In total, 1,318 and 835 metabolites were detected in positive and negative ion modes, respectively. Within the positive ion mode, 219 metabolites were significantly down-regulated, while 109 were significantly up-regulated in the stool samples of SSc patients (Figure 2F). In the negative ion mode, 189 metabolites were significantly down-regulated, and 53 metabolites were significantly up-regulated in the feces of SSc patients (Figure 2H). A heatmap analysis illustrated distinct patterns of differentially expressed metabolites (DEMs) accumulation within each group, with a clear distinction between the HC and SSC groups (Figures 2I,J). Subsequently, we conducted a KEGG enrichment analysis of the differential metabolites. The differentially expressed metabolites were mainly involved in the protein digestion and absorption, aminoacyl-tRNA biosynthesis, and biosynthesis of unsaturated fatty acids (Figures 2K,L). These data reveal that gut microbial metabolism in SSc patients were dysregulated.

As previously outlined, hUC-MSCs were further characterized in accordance with the criteria established by the International Society of Cellular Therapy (Dominici et al., 2006). After approximately 2 days in culture, the MSCs adhered to the culture flask, exhibiting a spindle-shaped morphology. The surface markers of the hUC-MSCs were detected by flow cytometry, revealing high expression levels of CD44 (100.00%), CD73 (100.00%), CD90 (99.99%), and CD105 (99.64%) were highly expressed in hUC-MSCs, while CD45 (0.33%) and CD34 (0.13%) were minimally expressed (Figure 3A). Following adipogenic, chondrogenic, and osteogenic inductions, the hUC-MSCs showed strong capacities of differentiation into adipocytes, chondrocytes, and osteocytes, respectively, as evidenced by Oil Red O staining (Figure 3B), Alcian Blue staining (Figure 3C), and Alizarin Red staining (Figure 3D).

Following a one-week acclimatization period, the C57BL/6 mice were randomly divided into three groups (n = 6 per group): healthy control (CON), bleomycin treatment (SSC), and MSCs therapy (MSC) groups (Figure 4A). Lung tissues from the mice were subsequently collected for pathological staining analysis. Compared to the CON group, H&E and Masson staining revealed a significant increase in inflammatory cell infiltration and collagen deposition in the lung tissues of the SSC group, whereas these symptoms were markedly reduced following MSCs treatment (Figures 4B,C). Immunohistochemistry (IHC) results indicated an upregulation of fibrosis-related proteins α-SMA and COL1A1 in the lung tissues of the SSC group, which was subsequently diminished in the MSC group (Figures 4D,E). These results suggest that MSCs therapy may mitigate pulmonary inflammatory injury and fibrosis in SSc mice.

To further investigate the molecular mechanisms underlying MSCs therapy in SSc-induced pulmonary fibrosis, transcriptome sequencing was performed in CON, SSc, and MSC -treated mice. PCA and PCoA demonstrated a significant divergence in gene expression patterns between the CON and SSC groups (Figures 5A,B), whereas the MSC-treated group exhibited gene expression profiles more closely aligned with the CON group. This observation suggests that MSC treatment may modulate gene expression in SSc mice. Volcano plot analysis was used to find differentially expressed genes (DEGs) (p < 0.05, FC = 2). Compared with the CON mice, 222 genes were significantly down-regulated and 1,244 genes were significantly up-regulated in the SSc mice (Figure 5C). Compared with the SSC mice, 146 genes were significantly down-regulated and 164 genes were significantly up-regulated after MSCs treatment (Figure 5D). Further analysis revealed that MSC treatment reversed the expression of specific genes in the SSC group, with 105 genes initially up-regulated becoming down-regulated, and 17 genes initially down-regulated becoming up-regulated (Figures 5E–G).

GO and KEGG analyses were performed to delineate the significant functional of DEGs induced by SSC and MSCs treatment. GO enrichment analysis showed that genes were significantly up-regulated in the SSC group were primarily associated with the fatty acid metabolic process, signal release, hormone transport, and cellular components, including the extracellular matrix, proteinaceous extracellular matrix, and collagen trimer (Figure 5H). Conversely, genes were significantly down-regulated following MSCs treatment were mainly enriched in the fatty acid metabolic process, regulation of lipid metabolic process, fatty acid oxidation, lipid modification, and regulation of fatty acid oxidation (Figure 5J). The KEGG pathway enrichment analysis indicated that pathways such as the AMPK signaling pathway, PPAR signaling pathway, fatty acid biosynthesis, carbon metabolism, and cytokine−cytokine receptor interaction were significantly up-regulated in the SSC group, but were significantly down-regulated following MSC treatment (Figures 5I,K). These findings suggest that MSCs may alleviate pulmonary fibrosis in SSC mice by regulating these pathways.

Compared with the CON mice, the α-diversity indices (observed species, Shannon, Chao1, and ACE) of the gut microbial community were significantly higher in SSc mice (SSC), while the administrations of MSCs could reduce these changes (Figure 6A; Supplementary Figure S4A). PCoA at both the phylum and species levels demonstrated a significant separation of the gut microbial communities in the SSC group. The microbial communities of the CON and SSC groups were significantly separated, whereas the microbial communities in the MSC group exhibited a divergence from the SSc group, showing a greater similarity to the CON group (PERMANOVA: p = 0.0001, Figure 6B; Supplementary Figure S4B). Consistent findings were observed in the group-paired PCoA analysis (p < 0.05, Supplementary Figures S4C–E). Clustering analysis further revealed that the microbial communities of the CON and MSC groups clustered into a single large branch, whereas the SSC group formed a separate branch (Figure 6C; Supplementary Figure S4F).

Similar to the human gut microbiota, Bacteroidota (54.35%) and Bacillota (13.54%) were the main microbial phylum in mice. At the phylum level, compared with the CON group, the relative abundances of Bacteroidota in the SSC group were significantly decreased, while that of Verrucomicrobiota, Actinomycetota, and Candidatus Saccharibacteria were significantly increased, and then increased or decreased after MSCs treatment (Figures 6D,F). At the genus level, compared with CON group, the relative abundance of Bacteroidesin and Duncaniella in the SSC group was significantly decreased, while that of Akkermansia, Bifidobacterium, Faecalibaculum, and Lactobacillus in SSC group was significantly increased, and then increased or decreased after MSCs treatment (Figures 6E,G). At the species, compared with CON mice, the relative abundance of Akkermansia muciniphila, Bifidobacterium pseudolongum, Erysipelotrichaceae bacterium, Lactobacillus johnsonii, Muribaculaceae bacterium Isolate-110 (HZI), and Faecalibaculum rodentium increased in SSc mice, then decreased with MSCs treated (Supplementary Figures S4G,H). Linear discriminant analysis effect size (LEfSe) results showed that Bacteroidales, Bacteroides, Helicobacter ganmani, Duncaniella freteri, Muribaculaceae bacterium, bacterium_1xD42 87, Bacteroides_sp_CAG_927, and Bacteroides acidifaciens were significantly enriched in the CON group. Erysipelotrichia, Coriobacteriales, Burkholderiales, Akkermansiaceae, Akkermansia, Faecalibaculum, Lactobacillus, Lactobacillus_johnsonii, Faecalibaculum rodentium, Ligilactobacillus murinus, Bifidobacterium pseudolongum, Erysipelotrichaceae bacterium, Limosilactobacillus reuteri, and Akkermansia muciniphila were significantly enriched in the SSC group. Chlamydiia, Chlamydiales, Chlamydiaceae, Chlamydia, and Chlamydia abortus were significantly enriched in the MSC group (Figure 6H). Overall, these results had demonstrated that the structure of gut microbiota in SSc mice was significantly changed compared with that in the control mice, meaning that the gut microbiota of SSc mice is imbalanced. These results are consistent with those of patients with SSc. MSCs therapy could effectively affect the diversity and composition of gut microbiota in SSc mice. As a result, the administrations of MSCs help to alleviate gut microbiota dysbiosis in SSc mice.

To investigate metabolite changes induced by gut microbiota after MSCs treatment, we detected non-targeted metabolites by UHPLC–MS/MS. Similar to the results in human samples, PCA and PLS–DA analysis showed that the gut microbiota metabolite composition between the CON and SSC groups was significantly separated (Figures 7A–D), indicates that the gut microbial metabolism dysregulated in SSc mice. The gut microbiota metabolite composition of MSC group were more similar with CON group, this implies that MSCs treatment may help to altered the dysregulated of gut microbial metabolites in SSc mice. 1,095 and 573 metabolites were detected in the positive and the negative ion modes, respectively. Compared to the CON group, BLM treatment resulted in significant alterations in metabolite levels, with 137 (positive: 71, negative:66) metabolites up-regulated and 275 (positive: 222, negative: 53) metabolites down-regulated (Figures 7E,G). MSCs treatment increased the levels of 128 (positive: 78, negative: 50) metabolites and decreased the levels of 94 (positive: 63, negative:31) metabolites (Figures 7F,H). To refine the list of DEMs, a Venn analysis was conducted across the three groups. MSCs treatment significantly altered the expression trends of 58 metabolites, which were significantly up-regulated or down-regulated in the SSC group (Figures 7I,J). The heatmap analysis demonstrated that the microbial metabolites from the CON and MSC groups clustered into a single large branch, whereas the SSC group formed a distinct branch (Figure 7K). KEGG enrichment analysis results showed that the DEMs were mainly involved in nicotinate and nicotinamide metabolism, pantothenate and CoA biosynthesis, tyrosine metabolism, neuroactive ligand−receptor interaction, inflammatory mediator regulation of TRP channels (Figures 7L,M). These data indicated that MSCs treatment could alleviate metabolic disorders by altering the gut microbial metabolic process in SSc.

We performed Spearman correlation analysis to explore the potential relationships among gut microbiota, metabolites, and gene expression in lung tissue. As shown in the heatmap, we screened 20 genes involved in fatty acid metabolic process, AMPK signaling pathway, PPAR signaling pathway, fatty acid biosynthesis, carbon metabolism, and cytokine−cytokine receptor interaction carbon metabolism from DEGs (Figure 8A, left). Spearman correlation analysis showed that the expression levels of DEGs (including Lep, Nodal, Il17a, Eno1b, Tkt, Eno1, Acacb, and Gys2) were significantly correlated with the relative abundance of differential microbial taxa (including Candidatus Saccharibacteria, Actinomycetota, Bifidobacterium, Faecalibaculum, Lactobacillus, Bacteroides, Bifidobacterium pseudolongum, Faecalibaculum rodentium, and Erysipelotrichaceae bacterium) (Figure 8A left & middle), these results suggest that the gut microbiota may influence fibrosis by affecting the expression of genes in lung tissue. Furthermore, the expression levels of the 20 genes and microbial metabolites (2-Thio-acetyl MAGE, NAOrn 15:0/17:0, AMK, (±)11(12)-DiHET and 16(R)-HETE) also showed a significant correlation (Figure 8A left & right). These results indicate that gut microbiota metabolites influence tissue fibrosis in mice. To further explore the potential relationships among the multi-omics data, we conducted a Mantel test, and the results showedthat, with the exception of Muribaculaceae bacterium Isolate-110 (HZI), there was a significant correlation between differential microbiota and microbial metabolites (Figure 8B). Consequently, the mediating effect of gut microbiota on fibrosis-related pathways in SSc mice may be attributed to its potential regulatory influence on microbial metabolites. In summary, MSCs may alleviated pulmonary fibrosis in SSc mice by regulating gut microbiota and microbial metabolites.