A sample of LG-ESS was collected immediately after surgery, some of which were fixed with 10% formalin fixation, and some were snap frozen. The formalin fixed sample was further embedded in paraffin for immunohistochemical staining. The snap frozen sample was digested and dispersed to form a single-cell suspension. Single cell capture, reverse transcription, amplification, and library construction were performed using the 10X Genomics platform according to the manufacturer’s instructions. Libraries were then sequenced, and data was processed using CellRanger software to generate a single cell transcriptome data matrix. The study was approved by the Ethics Committee of Zunyi Medical University (No. 2020-1-013), with informed consent obtained from all participants.

Single-cell transcriptome sequencing data for five healthy control (HC) endometrial samples (raw fastq files) were downloaded from the NCBI SRA database, accession number SRP349751 (23). Bulk transcriptomic datasets used for validation were sourced from the GEO database (GSE128630, GSE119041, and GSE85383). The GSE128630 dataset encompasses 75 ESS cases (24). GSE119041 includes 13 leiomyosarcoma (LMS), 16 ESS, 26 UUS, and three YWHAE-FAM22 endometrial stromal sarcomas (YFAM) samples, along with 14 benign leiomyoma controls (LM) (25). GSE85383 contains 9 LG-ESS, 4 HG-ESS, 8 UUS, and 4 uterine LMS (26). The clinical characteristics of these samples were extracted from the corresponding publications.

DoubletFinder (version 2.0.3) (27) was used for each sample to identify and remove doublets, with the doublet rate set to 0.03 for 10x genomic sequencing.

Data integration of Seurat objects after doublet removal was performed using the Seurat (R package,version 4.4.0) (28), filtering out abnormal cells (mito_percent < 15% & nFeature > 200 & nFeature < 6000) and logarithmizing the UMIcount matrix. In order to eliminate experimental batch effects, the harmony package (version 1.0.3) was used for the integration of the data (29).

The top 2000 highly variable genes were selected for principal component analysis (PCA), with principal components determined through the JackStraw procedure. After determining the selected number of principal components (dims = 1:30) during the integration process, we clustered the cells with a resolution of 1 and performed dimensionality reduction with t-SNE for visualization. Cell subpopulations were named based on SingleR (version 2.2.0) annotation and classic marker expression (23, 30).

We performed RNA velocity analysis based on bam files generated from CellRanger Count (version 7.2.0) analysis using the scVelo software (version 0.2.5) (31). scVelo infers temporal dynamic changes in gene expression by estimating splicing dynamics and unspliced/spliced ratios of individual cells. Finally, based on the RNA velocity results, we can reveal the dynamic changes in cell status and their potential differentiation pathways in the LG-ESS tumor sample.

For trajectory analysis, Monocle 2 (version 2.18.0) (32) was used to investigate dynamic changes in cell states based on the expression of variable genes in the target population (cells expressing ≥10 genes with an average expression >0.5). Genes with significant temporal variation were identified using the differentialGeneTest function. Dimensionality reduction was performed with the DDRTree algorithm, and a minimum spanning tree was constructed using the plot_cell_trajectory function. To visualize gene expression changes along pseudotime, the top 100 dynamically expressed genes were clustered using the plot_pseudotime_heatmap function. Functional enrichment analysis for each gene cluster was conducted using KOBAS 3.0.

In order to investigate the interactions between different cell types in the tumor microenvironment of LG-ESS, we performed intercellular communication analysis using CellPhoneDB2 (database version v2.0.0) (33), a Python-based computational analysis tool, to identify interaction networks.

To quantify the cell type composition within our bulk transcriptome samples, we utilized the BisqueRNA (version 1.0.5) algorithm (34). This tool allowed us to leverage the gene expression profiles derived from identified subpopulations in single-cell RNA sequencing data. Through this deconvolution approach, we were able to infer the relative abundances of different cell types present in the mixed tissue samples.

We performed gene set enrichment analysis (GSEA) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on the differential gene expression results using the clusterProfiler (version 4.8.3) (35) and kobas (version 3.0) software, respectively.

Survival analysis: We used the cutoff values (median expression or best cutoff value) of the expression of specific genes in the GSE128630 cohort and the GSE119041 cohort to divide the cohort into two groups and used the survival package to compare the differences in survival rates between the groups and to plot Kaplan-Meier (KM) curves.

Univariate Cox regression analysis: Cox regression modeling was performed using the survival package of the R software, and one-way Cox regression was performed with the expression of the gene of interest as a covariate.

Multivariate Cox regression analysis: Cox regression modeling was performed using the survival package in R software, incorporating relevant clinical factors such as age at diagnosis, tumor grade, and BMI, alongside the expression levels of FGF12 and KLHL29 as covariates. This approach allowed us to assess the independent prognostic significance of FGF12 and KLHL29 while adjusting for the effects of these clinical variables.

The paraffin-embedded tissue samples (FFPE) are cut into 3-5 micron sections. For H&E staining, the sections are stained in an acid hematoxylin solution for about 8 minutes and then in an eosin solution for about 2.5 minutes using a slide stainer.

A 7-color multiplex immunofluorescence staining was performed using the OPAL™ multiplexing method. The staining protocol for FFPE tissue sections was optimized for the simultaneous detection of 6 antibodies and DAPI as a nuclear stain. The fixed tissues were cut into 4µm sections by Rotary Microtome (Leica RM2255, Germany). The sections were deparaffinized, rehydrated, subjected to heat-induced epitope retrieval and incubated with primary and secondary antibodies. The antibodies were visualized using the fluorescent tyramide from the Opal 6-Plex Manual Detection Kit (NEL861001KT, Akoya Biosciences). The process of epitope retrieval and staining was repeated sequentially for different primary antibodies and fluorescent tyramide combinations. The following primary antibodies with different dilutions were used: CD10 (#110M-16, Cell Marque) with 1:40 dilution, Desmin (#243M-14, Cell Marque) with 1:50 dilution, CD4 (#104R-14, Cell Marque) with 1:100 dilution, CD8 (#108R-14, Cell Marque) with 1:75 dilution, CD68 (#168M-94, Cell Marque) with 1:75 dilution, FOX-P3 (#14-4777-82, eBioscience). Antibodies were visualized with the following tyramide dyes from the Opal Detection kit (NEL861001KT, Akoya Biosciences): Opal Polaris 480, Opal 520, Opal 570, Opal 620, Opal 690, and DIG-Opal 780. Sections were mounted with ProLong^® Diamond Antifade Mountant (P36961, Thermo Fisher Scientific). Multiplex-stained slides were imaged using a PhenoImager Fusion system (Akoya Biosciences).

The log-rank test was used for survival analysis. For other intergroup comparisons, the Wilcoxon rank-sum test was used. P<0.05 is considered to be statistically significant.

To investigate the cellular composition of LG-ESS, we collected a fresh tissue sample from an LG-ESS patient during surgery and performed single-cell transcriptomic sequencing. Concurrently, we rigorously curated and integrated this dataset with single-cell transcriptomic data from five healthy controls (HCs) sourced from a public database, ensuring stringent quality control throughout the analysis ( Figure 1A ).

Following batch correction ( Supplementary Figure S1A ), we dimensionally reduced and clustered 38,288 cells, identifying 12 primary cell subpopulations ( Figure 1B ): non-ciliated epithelial cells (Epi, EPCAM⁺ SNTN^-), ciliated epithelial cells (Cili_Epi, EPCAM⁺ SNTN⁺), Str (stromal fibroblasts, DCN⁺ MME⁺ LUM⁺⁾, smooth muscle cells (Smooth_muscle_cell, ACTA2⁺⁾, proliferative stromal cells (pStr, MKI67⁺⁾, pericytes (Peri, NOTCH3⁺ RGS5⁺⁾, endothelial cells (Endo, PECAM1⁺⁾, lymphatic endothelial cells (Lymph, CCL21⁺⁾, macrophages (Mac, CD163⁺⁾, T cells (Tcell, CD3D⁺⁾, natural killer cells (NK, GNLY⁺⁾, and mast cells (Mast_cell, TPSB2⁺⁾ ( Figure 1C ). Supplementary Figure S1B displays the top five highly expressed genes in each cell subpopulation, and Supplementary Figure S1C shows the quantification of all cell subpopulations from each sample. Furthermore, differential gene analysis of all subpopulations between LG-ESS and HC revealed that multiple subpopulations exhibited significant differences in gene expression, particularly within the Str subpopulation ( Figure 1D ). We speculate that changes in the gene expression profile of this subpopulation may be intricately linked to the pathogenesis of LG-ESS.

Previous studies have indicated that most LG-ESS exhibit positive expression of CD10 (protein of MME gene) and WT1, with elevated levels of estrogen receptor expression (20, 36, 37). Through H&E staining ( Figure 2A ) and multiplex immunofluorescence (mIF) staining ( Figures 2B, C ; Supplementary Figure S2A ) of adjacent tissue sections, we found that the tumor cell-rich regions of the LG-ESS sample displayed positive CD10 expression ( Figure 2B , right middle) and were located adjacent to smooth muscle cell populations (α-SMA⁺ DES⁺) ( Figure 2B ). Additionally, the tumor sample contained various other cell types, including epithelial cells (EPCAM⁺), pStr (MKI67⁺) CD4⁺T cells (CD4⁺), CD8⁺ T cells (CD8A⁺), macrophages (CD68⁺), and regulatory T cells (FOXP3⁺), which were further validated by our single-cell data ( Figure 2D ). Unlike most epithelial-derived solid tumors (including endometrial carcinoma (22)), LG-ESS originates from mesenchymal fibroblasts rather than endometrial epithelial cells (4). Consistent with this notion, we found that MME⁺(CD10) was predominantly expressed in the Str subpopulation, with almost no expression in epithelial cells.

To explore the functional roles of tumor cell subpopulations in LG-ESS, we re-clustered the Str subpopulation ( Figure 2E ), resulting in eight distinct clusters. We observed that Cluster2 and Cluster4 were primarily composed of LG-ESS, and notably, LG-ESS tumor cells (MME⁺) were mainly located in Cluster4 ( Figure 2F ). Consequently, we independently clustered the LG-ESS samples and further differentiated the Str subpopulation into two subgroups: Str1, characterized by MME⁺ tumor cells (predominantly in Cluster4), and Str2, defined by MME- cells (predominantly in Cluster3) ( Figures 2G, H ). KEGG enrichment analysis indicated that the highly expressed genes in the Str1 tumor cell subgroup were significantly enriched in several classic cancer-related pathways, including focal adhesion, choline metabolism in cancer, Wnt signaling pathway, and PI3K-Akt signaling pathway ( Figure 2I ). This aligns with previous findings demonstrating significant activation of the Wnt signaling pathway in LG-ESS (38).

Subsequently, we utilized BisqueRNA to estimate the proportions of cells in a large number of transcriptomic samples based on the top 50 most highly expressed genes across all LG-ESS cell subpopulations, and we explored the relationship between subpopulation proportions and the prognosis of LG-ESS. We found that an increased proportion of Str1 in LG-ESS was significantly associated with shorter overall survival and progression-free survival of patients ( Supplementary Figure S2B ). At the same time, we observed no significant differences in the proportion of Str1 across different sarcoma types ( Supplementary Figure S2C ).

To further elucidate the origins of the Str1 tumor cell subpopulation in LG-ESS, we conducted RNA velocity analysis on both Str1 and Str2 subpopulations to determine their lineage relationships. Surprisingly, our findings indicated that Str1 and Str2 may share a common ancestral origin, contradicting our initial hypothesis that Str1 arose from Str2 ( Figure 3A ).

To validate this finding, we utilized the Str subpopulation from HC as a normal reference and performed RNA velocity analysis on all Str cells from both HC and LG-ESS samples. We discovered that Str1 (predominantly corresponding to Cluster4) and Str2 (primarily corresponding to Cluster2) indeed originated from the same ancestral cluster, Cluster3. Moreover, Cluster4 originating from Cluster3 through the intermediate Cluster6 ( Figure 3B ). Simultaneously, we conducted pseudo-time analysis on Clusters3, Clusters4, and Clusters6 using Monocle2. Consistent with the RNA velocity analysis, our results confirmed that Cluster4, which primarily comprises the Str1 tumor cell subpopulation in LG-ESS, predominantly develops from Cluster3 and Cluster6 ( Figure 3C ).

To explore the characteristic changes occurring during the differentiation of tumor cells, we clustered the genes that exhibited significant changes throughout the differentiation process ( Figure 3D ). KEGG enrichment analysis revealed that genes that were significantly upregulated during differentiation were primarily associated with common cancer-related pathways, such as focal adhesion, proteoglycans in cancer, Wnt signaling pathway, and Jak-STAT signaling pathway. Conversely, genes that were significantly downregulated during differentiation were predominantly linked to metabolic processes and oxidative phosphorylation pathways ( Figure 3E ).

Upon investigating the biological alterations within tumor cell subpopulations of ESS as compared to HC, our enrichment analyses of differentially expressed genes (LG-ESS vs. HC) within the Str-Cluster4 subpopulation ( Figure 2E ) indicated significant gene expression modifications. These modifications predominantly clustered around pathways integral to oncogenesis and tumor progression, including growth hormone synthesis, secretion, and action, along with cell cycle regulation, focal adhesion, Wnt, and MAPK signaling pathways ( Figure 4A ). To investigate the underlying drivers of these pathways, we performed Cox regression and survival analysis on these significantly differentially expressed genes based on the gene expression data of ESS patients from the GEO database, focusing on genes uniquely expressed in Str1 ( Figure 4B ). Notably, high expression of FGF12 and KLHL29 in ESS was significantly associated with shorter overall and progression-free survival ( Figures 4C-E ). Through multivariate cox regression analysis, we found that FGF12 and KLHL29 remain important independent prognostic factors affecting overall survival ( Supplementary Figure S3A ) and progression-free survival ( Supplementary Figure S3B ). In addition, an independent cohort of patients further validated the association between high FGF12 gene expression and shortened overall survival (OS) in uterine sarcoma patients ( Supplementary Figure S3C ), suggesting that FGF12 could potentially be a prognostic biomarker for ESS.

Moreover, a comparative analysis of FGF12 and KLHL29 expression across different sarcoma subtypes and leiomyosarcomas (LMS) revealed a significant difference in KLHL29 expression between HG-ESS and LG-ESS and LMS. In contrast, FGF12 expression showed no significant difference. This indicates that KLHL29 may serve as a potential biomarker for distinguishing HG-ESS from LG-ESS ( Supplementary Figure S3D ).

To further elucidate the biological functions of FGF12 and KLHL29 in ESS, we categorized ESS samples into high- and low-expression groups based on their gene expression levels and conducted GSEA enrichment analysis. We observed that the MYOGENESIS pathway was significantly activated in both the FGF12 ^high and KLHL29 ^high groups. Conversely, immune-related pathways, such as the INFLAMMATORY RESPONSE and INTERFERON GAMMA RESPONSE, were notably suppressed in the FGF12 ^high group, while the ESTROGEN RESPONSE EARLY pathway was significantly suppressed in the KLHL29 ^high group ( Figures 4F, G ).

Previous studies have shown that in epithelial-derived tumors, fibroblasts within the tumor microenvironment may be induced into cancer-associated fibroblasts (CAFs) through interactions with tumor cells (39, 40), and macrophages may undergo M2 polarization under the influence of tumor-derived factors (41). These phenotypic and functional changes in stromal and immune cells within the tumor microenvironment, in turn, contribute to tumor progression.

Through intercellular communication analysis, we identified a strong interaction between Str1 and Str2 as well as smooth muscle cell subpopulations in LG-ESS, particularly mediated by secreted growth factors ( Figures 5A, B ). We hypothesize that these interactions may lead to functional alterations in the subpopulations, thereby promoting tumor progression. To further explore these potential functional changes, we performed KEGG enrichment analysis on the highly expressed genes in the Str2 and smooth muscle cell subpopulations in LG-ESS. Interestingly, the genes predominantly clustered in pathways associated with tumorigenesis and metastasis, such as focal adhesion, proteoglycans in cancer, PI3K-Akt signaling, and ECM-receptor interactions ( Figures 5C, D ). Differential gene expression analysis between LG-ESS and healthy controls further supported these findings, indicating that the functions of Str2 and smooth muscle cells in LG-ESS are altered under the influence of tumor cells ( Figures 5E, F ).