Spatial and single-cell transcriptomics landscape of adenomyosis

The funders played no role in the design, experiment conduction, data analysis or manuscript preparation.

A total of 28 participants were enrolled in this study, including 22 patients diagnosed with adenomyosis based on medical imaging and histopathological examination (Sup. Table S1) . Among the adenomyosis cases, eight patients had not received hormonal treatment within the three months prior to surgery (Non-GnRHa-treatment Adenomyosis group, abbreviated as NGAM), whereas 14 patients underwent pre-surgical GnRHa treatment to manage symptoms and reduce uterine size, thereby facilitating surgical intervention (GnRHa-treatment Adenomyosis group, abbreviated as GTAM). In addition, six control participants who underwent hysterectomy for non-adenomyosis indications, and had not received hormonal treatments for at least three months prior to surgery (CTRL group). Consequently, tissue samples from patients in the NGAM group represented the natural proliferative phase of the menstrual cycle, while those from the GTAM group reflected a pharmacologically modified proliferative phase. All control samples were also obtained during the natural proliferative phase. Notably, high-throughput sequencing was performed on samples from 15 of the 28 participants, wheras the remaining samples were used for validation analyses. Furthermore, we integrated previously published scRNA-seq datasets[ 14 , 15 , 16 , 17 , 18 ] from samples in the proliferative phase of the menstrual cycle, including data from 13 healthy uterine samples and three patients with adenomyosis without prior medication, to validate our findings (Sup. Table S3) . To elucidate the cellular and molecular mechanism underlying adenomyosis, we performed scRNA-seq on the entire uterine tissue samples, spaning from endometrium to perimetrium. Samples were collected from 11 individuals, including three CTRL patients, three NGAM adenomyosis patients and five GTAM adenomyosis patients ( Fig. 1 , Fig. 1 and Sup. Table S1) . After quality control, a total of 100,404 high-quality single cells were obtained, with a median 8,381 unique transcripts and 2,841 detected genes per cell (Supplementary Fig. S2a) . Canonical correlation analysis (CCA) was used to integrate individual samples and identify common sources of variation. For data normalization, we applied a global-scaling method to normalize feature expression levles for each cell based on the total expression, followed by log-transformation. Based on the expression of known markers and the cell-type prediction with published references[ 15 , 28 ], we identified six major cell clusters: epithelial cells ( EPCAM and CDH1 ), endothelial cells ( VWF and PECAM1 ), smooth muscle cells (SMCs; ACTA2 and MYH11 ), fibroblasts ( PDGFRA and APOD ), myeloid cells ( PTPRC, CD68 and CD14 ), and lymphocytes ( PTPRC, CD3D, CD3E and CD2 ; Fig. 1 , Fig. 1 and Supplementary Fig. S2b ). The cell-type composition was examined at the patient level (Supplementary Fig. S2c and 2d), and clustering performance was evaluated using the adjusted rand index (ARI) and normalized mutual information (NMI)[ 11 ]. Using the pre-integration cell-type composition as the ground truth, CCA integration achieved excellent clustering performance (Supplementary Fig. S2e) . Among these identified cell types, fibroblasts were the most abundant, accounting for approximately 30 % of all cell (Supplementary Fig. S2f) . Endothelial cells and SMCs comprised 28 % and 27 %, respectively. Immune cells, including myeloid (3 %) and lymphoid (6 %) populations, represented a smaller fraction, collectively constituting about 9 % of total cell population (Supplementary Fig. S2f) . Together, these analyses established a comprehensive single-cell atlas of the human uterus across control and adenomyosis conditions, providing a robust foundation for dissecting cell-type-specific alterations and the molecular events driving disease progression. Fig. 1 Single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (Geo-seq) analysis of control and adenomyosis samples. (a) Schematic depicting samples acquisition, processing and bioinformatics analysis. (b) Block chart summarizing sampling details for each patient. (c) Uniform manifold approximation and projection (UMAP) plot showing 100,404 single cells from control (CTRL, n  = 3) and adenomyosis individuals (NGAM, n  = 3; GTAM, n  = 5), clustered into six main groups. (d) Each subpopulation was identified with the indicated markers. (e) UMAP plot showing 300 spatially resolved tissue specimens collected from uterine tissues of 12 individuals (CTRL, n  = 4; NGAM, n  = 2; GTAM, n  = 6). (f) Schematic digram depicting the relative position of 300 tissue specimens. Pie charts indicate the proportion of each cell type within regions, deconvoluted using CIBERSORT. (g) Bar chart showing the distribution of cell-type proportion in scRNA-seq dataset across different uterine regions. Con Gland, control endometrial gland; Con Mesen, control endometrial mesenchyma; Con Myo, control myometrium; Myomesen, control myometrial mesenchyma; Eut Gland, eutopic endometrial gland; Eut Mesen, eutopic endometrial mesenchyma; JZ, adenomyosis junctional zone; Ect Gland, ectopic endometrial gland; Ect Mesen, ectopic endometrial mesenchyma; Outer Myo, adenomyosis outer myometrium; Outer Mesen, adenomyosis myometrial mesenchyma; CTRL, Control group; NGAM, Non-GnRHa treatment adenomyosis group; GTAM, GnRHa treatment adenomyosis group. Single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (Geo-seq) analysis of control and adenomyosis samples. (a) Schematic depicting samples acquisition, processing and bioinformatics analysis. (b) Block chart summarizing sampling details for each patient. (c) Uniform manifold approximation and projection (UMAP) plot showing 100,404 single cells from control (CTRL, n  = 3) and adenomyosis individuals (NGAM, n  = 3; GTAM, n  = 5), clustered into six main groups. (d) Each subpopulation was identified with the indicated markers. (e) UMAP plot showing 300 spatially resolved tissue specimens collected from uterine tissues of 12 individuals (CTRL, n  = 4; NGAM, n  = 2; GTAM, n  = 6). (f) Schematic digram depicting the relative position of 300 tissue specimens. Pie charts indicate the proportion of each cell type within regions, deconvoluted using CIBERSORT. (g) Bar chart showing the distribution of cell-type proportion in scRNA-seq dataset across different uterine regions. Con Gland, control endometrial gland; Con Mesen, control endometrial mesenchyma; Con Myo, control myometrium; Myomesen, control myometrial mesenchyma; Eut Gland, eutopic endometrial gland; Eut Mesen, eutopic endometrial mesenchyma; JZ, adenomyosis junctional zone; Ect Gland, ectopic endometrial gland; Ect Mesen, ectopic endometrial mesenchyma; Outer Myo, adenomyosis outer myometrium; Outer Mesen, adenomyosis myometrial mesenchyma; CTRL, Control group; NGAM, Non-GnRHa treatment adenomyosis group; GTAM, GnRHa treatment adenomyosis group. We next compared the proportions of various cell types among CTRL, NGAM and GTAM groups.The endothelium showed an increasing trend in the NGAM group, suggesting enhanced angiogenic activity within adenomyotic lesions (Supplementary Fig. S2 g) . Differentially expressed gene (DEG) analysis between CTRL and NGAM groups revealed increased expression of genes, such as ICAM1 and CCL2 in adenomyosis, consistent with previous reports[ 29 , 30 ] (Supplementary Fig. S2h and Sup. Table S4) . The upregulated genes were enriched in pathways related to tissue development and regulation of cell migration (Supplementary Fig. S2i) . Notably, genes including CCL2 , SOD2 , YBX3 and LMNA were upregulated in the NGAM group (Supplementary Fig. S2h). This pattern suggests that adenomyosis is characterized by the activation of inflammatory processes and epithelial-mesenchymal transition (EMT)[ 31 ]. Consistent with previous studies[ 7 , 8 ], most patients included in the scRNA-seq analysis also suffered from comorbidities, such as endometriosis or leiomyomas. To clarify the potential impact of these complications on our conclusions and enhance the stringency of our study, we integrated previously published scRNA-seq datasets from healthy uterine tissues[ 14 , 15 , 16 , 17 , 18 ], yielding a combined dataset of 204,627 cells. These cells were divided into three groups: the control group, the adenomyosis group (without complications), and the complications group (adenomyosis patients with additional complications). Cells were categorized into six subclusters based on defined marker genesr (Supplementary Fig. S3a-c) . No significant differences in overall cell-type composition were observed between the adenomyosis and complications groups (Supplementary Fig. S3d) . We then calculated Kolmogorov-Smirnov distances among the three groups for each cluster and found that the distance between the adenomyosis and complication groups was significantly smaller than that between the control group and either the adenomyosis or complications groups (Supplementary Fig. S3e and 3f) . Using the CINEMA-OT tool[ 19 ], we further assessed the influence of the adenomyosis and complications on specific cell types. The analysis revealed that fibroblast were the most affected clusters (Supplementary Fig. S3g-i) . Consequently, our downstream analysis primarily focuses on epithelial and immune cells. Because the ectopic endometrium in adenomyosis differs spatially from the eutopic endometrium, anatomical localization is essential for characterizing the pathophysiology of the disease. To systematically dissect the spatial organization of the endometrium, we performed spatial transcriptomic analysis using laser microdissection-based Geo-seq on 300 specimens obtained from uterine tissues of 12 individuals, including four controls, two adenomyosis patients without GnRHa treatment, and six adenomyosis patients treated with GnRHa ( Fig. 1b and Sup. Table S1) . Multiple spatial regions were collected, including the endometrial gland, endometrial mesenchyma, myometrium, outer myometrial mesenchyma in the CTRL group; and the eutopic endometrial gland, eutopic endometrial mesenchyma, junctional zone, ectopic endometrial gland, ectopic endometrial mesenchyma, outer myometrium, outer myometrial mesenchyma in the NGAM and GTAM groups ( Fig. 1b , Supplementary Fig. S4a and 4b) . Sampling was guided by precise porphological information to ensure accure anatomical localization. It is important to note that these spatial sampling procedures were applied exclusively to the spatial transcriptomic dataset. Alignment of the raw reads to the human reference genome revealed an average detection of approximately 15, 000 unique genes across CTRL, NGAM and GTAM groups (Supplementary Fig. S4c) . Geo-seq data were normalized using transcripts per million (TPM) to account for sequencing depth and transcript length, followed by log transformation. Uniform manifold approximation and projection (UMAP) based on Geo-seq profiles clearly delineated the gland, mesenchyma and myometrium regions, indicating that the spatial transcriptomic data faithfully captured the anatomical architecture of the uterine wall ( Fig. 1e ) . These results demonstrate that the Geo-seq–based spatial transcriptomic framework provides a high-resolution molecular map of the uterine tissues, enabling precise characterization of region-specific transcriptional programs. To maximize the strengths of both techniques and identify adenomyosis-associated cellular states that might not be detected by each method alone, we integrated the scRNA-seq and Geo-seq datasets through deconvolution analysis with CIBERSORT and spatialDWLS algorithms[ 21 , 22 ] ( Fig. 1f and Supplementary Fig. S4d) . We compared the proportions of major cell type across different anatomical regions and observed significant differences in cellular composition. As expected, gland regions were enriched in epithelial cells, myometrium regions in SMCs, and mesenchyma regions in fibroblasts, endothelial cells, lymphocytes and myeloid cells ( Fig. 1g ) . Consistent with the known effect of low oestrogen level following GnRHa treatment, a significant decrease in epithelial cells was observed in the endometrial glands of GTAM group (Supplementary Fig. S4e) . Notably, both myeloid and lymphoid populations exhibited marked alterations in the endometrial glands, mesenchyma and myometrium of adenomyosis tissues (Supplementary Fig. S4e) . Because scRNA-seq alone cannotdistinguish eutopic from ectopic endometrial glands, owing to the absence of major differences in overall cell-type distribution among the CTRL, NGAM and GTAM groups (Supplementary Fig. S2g) , we conducted spatially resolved Geo-seq analysis across distinct uterine regions in adenomyosis. Collectively, these results demonstrated that integrating scRNA-seq and Geo-seq enables the identification of tissue structure-assoicate cell types and provids a powerful framework for dissecting the molecular architecture and microenviromental niches underlying adenomyosis. Adenomyosis is characterized by the presence of endometrial glands within the myometrium, in which epithelial cells represent the predominant cell type ( Fig. 1g ) . To gain insight into the molecular alterations associated with adenomyosis, we focused our analysis on the epithelial components. Based on the scRNA-seq data, we identified three epithelial subtypes by assigning their characteristic marker genes ( Fig. 2 , Fig. 2 ) . These epithelial clusters corresponded to those previously described in the endometrium[ 16 , 32 ], including early proliferative epithelial cells (EarlyPro, PGR , ESR1 and IHH ), late proliferative epithelial cells (LatePro, GNG11 and CREB3L1 ), and ciliated epithelial cells ( DNAH9 , TPPP3 , PIFO and FOXJ1 ) ( Fig. 2b ) . Furthermore, nearly all unciliated epithelium could be classified as glandular cells ( FOXA2 and SMAD9 ), rather than luminal cells ( PTGS1, LGR5 and WNT7A ) based on established markers [ 33 ] ( Fig. 2b ) . All clusters exhibited similar proportions of cells from the CTRL, NGAM and GTAM groups, indicating an unbiased representation of epethelial subpopulations (Supplementary Fig. S5a) . Fig. 2 Epithelial cell heterogeneity in control and adenomyosis. (a) UMAP visualization of all epithelium subpopulations. (b) Bubble chart showing marker genes across epithelium subpopulations. (c) GSEA of epithelial cells based on scRNA-seq data comparing CTRL and NGAM groups. (d) Violin plots showing the area under the recovery curve (AUC) scores of marker genes associated with epithelial and mesenchymal signatures, as well as Wnt-catenin pathway, in scRNA-seq data comaring CTRL and NGAM groups. A two-tailed Mann–Whitney U test was conducted. (e) GO enrichment analysis comparing epithelial cells between eutopic and ectopic endometrial gland based on Geo-seq data. (f) Box plots showing the proportion of ciliated cells estimated by CIBERSORT deconvolution among control (n = 17), eutopic (NGAM n = 12; GTAM n = 28) and ectopic (NGAM n = 12; GTAM n = 20) endometrial glands (Kruskal-Wallis test with post-hoc analysis). (g) Representative immunofluorescence images showing co-localization of endometrial epithelial cells (Pan-cytokeratin, green) and ciliated cells (DNAH9, red) in ectopic endometrium (n = 8), but not in control (n = 5) or eutopic (n = 6) endometrium. Nuclei were counterstained with DAPI (cyan). Scale bar, 50  µm. (h) Quantification of multi-immunofluorescence (g) using HALO software (Kruskal-Wallis test with post-hoc analysis). (i) Ciliated cells identified in both our scRNA-seq dataset and that of Fonseca et al.[ 25 ] using established ciliated cell markers. (j) GSEA analysis of ciliated cells in adenomyosis showing higher enrichment of gene programs associated with the EMT pathway. P values were adjusted using the Benjamini–Hochberg method. Epithelial cell heterogeneity in control and adenomyosis. (a) UMAP visualization of all epithelium subpopulations. (b) Bubble chart showing marker genes across epithelium subpopulations. (c) GSEA of epithelial cells based on scRNA-seq data comparing CTRL and NGAM groups. (d) Violin plots showing the area under the recovery curve (AUC) scores of marker genes associated with epithelial and mesenchymal signatures, as well as Wnt-catenin pathway, in scRNA-seq data comaring CTRL and NGAM groups. A two-tailed Mann–Whitney U test was conducted. (e) GO enrichment analysis comparing epithelial cells between eutopic and ectopic endometrial gland based on Geo-seq data. (f) Box plots showing the proportion of ciliated cells estimated by CIBERSORT deconvolution among control (n = 17), eutopic (NGAM n = 12; GTAM n = 28) and ectopic (NGAM n = 12; GTAM n = 20) endometrial glands (Kruskal-Wallis test with post-hoc analysis). (g) Representative immunofluorescence images showing co-localization of endometrial epithelial cells (Pan-cytokeratin, green) and ciliated cells (DNAH9, red) in ectopic endometrium (n = 8), but not in control (n = 5) or eutopic (n = 6) endometrium. Nuclei were counterstained with DAPI (cyan). Scale bar, 50  µm. (h) Quantification of multi-immunofluorescence (g) using HALO software (Kruskal-Wallis test with post-hoc analysis). (i) Ciliated cells identified in both our scRNA-seq dataset and that of Fonseca et al.[ 25 ] using established ciliated cell markers. (j) GSEA analysis of ciliated cells in adenomyosis showing higher enrichment of gene programs associated with the EMT pathway. P values were adjusted using the Benjamini–Hochberg method. To investigate the biological characteristics of endometrial glands from patients with or without adenomyosis, we performed DEG and GSEA based on scRNA-seq data from the CTRL and NGAM groups. Notably, samples from adenomyosis patients were significantly enriched in gene signatures assocaited with EMT and the Wnt/β-catenin signaling pathway, further supporting the invagination theory of adenomyosis pathogenesis ( Fig. 2c ) . EMT plays a critical role in the development of adenomyosis, and abnormal activation of the Wnt/β-catenin pathway can induce EMT[ 34 ]. To validate these findings, we quantified epithelial and mesenchymal specificity scores in the scRNA-seq dataset ( Fig. 2d ) . Epithelial cells from adenomyosis samples showed a lower epithelial score and a higher mesenchymal score. Consistent with these results, the Wnt/β-catenin pathway was highly activated in adenomyosis ( Fig. 2d ) . We further validated these results by immunohistochemistry. Compared to the CTRL group, glandular cells in the NGAM group displayed elevated expression of the mesenchymal markervimentin and the key Wnt/β-catenin pathway component β-catenin, accompanied by reduced expression of the epithelial marker E-cadherin (Supplementary Fig. S5b and 5c) . Following GnRHa treatment, GSEA revealed no significant changes in the EMT and Wnt/β-catenin pathways (Supplementary Fig. S5d) . Gene signature specificity scores indicated that the GTAM group maintained a lower epithelial score and a higher mesenchymal score, whereas the Wnt/β-catenin pathway score remained unchanged (Supplementary Fig. S5e) . Thus, the mechanism by which GnRHa influences EMT remains unclear and warrants further investigation. Next, we focused on the transcriptome alterations between eutopic and ectopic endometrial glands based on Geo-seq data. We found that ectopic endometrial gland samples were significantly enriched for biological processes such as axoneme assembly, cilium movement, and microtubule-based movement ( Fig. 2e , Supplementary Fig. S5f and Sup. Table S5). By integrating the scRNA-seq and Geo-seq datasets, we observed that ciliated cells were prodeminantly localized in the ectopic endometrial glands of adenomyosis, irrespective of whether GnRHa treatment was administered ( Fig. 2f ) . Ciliated cells were previously described and verified to exist in the epithelium of the normal endometrium[ 14 ]. It has been reported that the differentiation of ciliated cells is induced by oestrogen and that the Wnt/β-catenin pathway is closely associated with ciliated cell differentiation[ 14 ]. We further confirmed the enrichment of ciliated cells in ectopic endometrium through multiplex immunofluorescence staining ( Fig. 2 , Fig. 2 , Supplementary Fig. S6 and Sup. Table S6) . Furthermore, we performed a comprehensive comparison between our adenomyosis dataset and previously published scRNA-seq dataset on endometriosis[ 25 ]. This published dataset also demonstrated a remarkable enrichment of signatures associated with ciliated epithelial cells[ 25 ]. By integrating both datasets corresponding to the proliferative phase of the menstrual cycle, we were able to distinguish ciliated cells from other epithelial cells in both adenomyosis and endometriosis samples ( Fig. 2i and Supplementary Fig. S5g) . Within these ciliated cells, we identified the specific expression of canonical marker genes FOXJ1 , PIFO and TPPP3 , further validating their cellular identity (Supplementary Fig. S5h) . A closer examination of gene expression profiles revealed distinct transcriptomic differences between ciliated cells in adenomyosis and those in endometriosis. In adenomyosis, ciliated cells exhibited higher enrichment of gene programs associated with the EMT pathway, suggesting a potential role of EMT-related mechanisms in the invagination process of the endometrium ( Fig. 2j ) . These findings highlight the unique molecular characteristics and potential functional roles of ciliated cells in adenomyosis compared with endometriosis. Taken together, by applying a deconvolution approach to integrate scRNA-seq with spatial Geo-seq data and immunofluorescence staining, we identified a ciliated cell population specifically associated with ectopic endometrial glands in adenomyosis. Futhermore, GnRHa administration did not reduce the abundance of ciliated cells or suppress the Wnt/β-catenin signaling pathway, which may contribute to disease recurrence after the discontinution of GnRHa treatment. Therefore, targeting the differentiation or formation of ciliated cells to prevent the establishment of ectopic lesions may represent a promising therapeutic strategy, potentially providing an alternative to conventional GnRHa therapy. Numerous studies have reported an aberrant immune microenvironment in adenomyosis[ 35 ]. To further investigate immune infiltration, we analyzed immune cell populations using scRNA-seq. This analysis revealed nine lymphocyte subpopulations and six myeloid subpopulations across all samples, including T cells, B cells, natural killer (NK) cells, macrophages, monocytes, dendritic cells (DC) and mast cells. All of these immune cell types were present in patients from CTRL, NGAM and GTAM groups ( Fig. 3 a-c) . CD4 + T cells were characterized by high expression of CD3D , CD3E , and CD4 , whereas CD8 + T cells expressed CD8A , CD8B and GZMK at high levels. The remaining T cells were categorized as other T-cell subsets. Among the NK clusters, some corresponded to NCAM1 + FCGR3A − KLRD1 + decidual NK subsets, specifically dNK1 and dNK2[ 36 ]. In addition, a proliferative NK subset with high MKI67 expression, designated as dNKp, was also identified. A small proportion of peripheral blood cells, such as peripheral blood NK cells (pbNK) marked by FCGR3A [ 36 ], were likely introduced during decidual tissue preparation. DCs were defined by high expression of CLEC9A , XCR1, CLEC9A, BATF3 and CADM1 . Two macrophage populations were identified: one with elevated IL1B expression (designated Macro-IL1B), and the other with high LYVE1 expression (designated Macro-LYVE1). Monocytes showed specific expression of S100A8 and S100A9 . Other clusters included neutrophil cells ( FCGR3B + CD79A + ), B cells ( MS4A1 + CD79A + ), mast cells ( CPA3 + KIT + ), and mucosal-associated invariant T (MAIT) cells enriched for SLC4A10 expression. To characterize the dynamic immune cell landscape in adenomyosis, we evaluated the distribution of lymphoid and myeloid cells across all samples. Our scRNA-seq data revealed significant differences in immune cell composition among the CTRL, NGAM and GTAM groups ( Fig. 3d and Supplementary Fig. S7a) . Compared with the CTRL group, a marked alteration was observed in CD4 + T cells within the NGAM group ( Fig. 3d ) . Using deconvolution analysis with CIBERSORT[ 21 ] , we further identified a pronounced enrichment of CD4 + T cells in the NGAM group ( Fig. 3e ) . Notably, elevated CD4 + T cell levels in the eutopic endometrium of patients with adenomyosis have been associated with an increased risk of infertility[ 37 ]. Following GnRHa treatment, the abundance of CD4 + T cells returned to levels comparable to those in the CTRL group, suggesting that a reduction in CD4 + T cells within the eutopic endometrium may improve clinical pregnancy outcomes. Fig. 3 The complex immune microenvironment in adenomyosis. (a) UMAP visualization of all lymphocyte and myeloid cell subpopulations. (b) Expression profiles of the marker genes of all lymphocyte and myeloid subpopulations. (c) UMAP visualization of CTRL, NGAM and GTAM groups. (d) Box plot showing the proportion of CD4 + T cells and Macro-LYVE1 cells among the CTRL ( n  = 3), NGAM ( n  = 3) and GTAM ( n  = 5) groups in scRNA-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (e) Deconvolution analysis using CIBERSORT showing enrichment of CD4 + T cells and Macro-LYVE1 cells in the NGAM group. Sample sizes were as follows: CTRL ( n  = 18), NGAM_eut ( n  = 12), NGAM_ect ( n  = 12), GTAM_eut ( n  = 20), and GTAM_ect ( n  = 20). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (f) Volcano plots showing DEGs in CD4 + T cells between CTRL and NGAM (left), and between GTAM and NGAM samples (right) (FDR  0.5). (g) Venn plot showing that most genes upregulated in the NGAM group were significantly downregulated in the GTAM group (Fisher’s exact test). (h) Box plot showing expression levels of CD44 and TGFB1 among the CTRL ( n  = 18), NGAM ( n  = 24) and GTAM ( n  = 40) groups in Geo-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (i) Reactome pathway analysis based on DEGs in (f, left) indicating T-cell activation in the NGAM group. P values were calculated using the hypergeometric distribution and adjusted using the false discovery rate method. (j) Heatmap showing AUC scores of marker genes related to M1 phenotype, proliferation, angiogenesis and phagocytosis of all the macrophages in scRNA-seq data among the CTRL, NGAM and GTAM groups. (k) GSEA analysis showing that IL-2/STAT5 signaling, TGF-β signaling, inflammatory response, and TNFα signaling via NF-κB pathways were upregulated in the NGAM group and downregulated following GnRHa treatment. P values were adjusted using the Benjamini–Hochberg method. NGAM_eut, eutopic endometrial mesenchyma of NGAM group; NGAM_ect, ectopic endometrial mesenchyma of NGAM group; GTAM_eut, eutopic endometrial mesenchyma of GTAM group; GTAM_ect, ectopic endometrial mesenchyma of GTAM group. The complex immune microenvironment in adenomyosis. (a) UMAP visualization of all lymphocyte and myeloid cell subpopulations. (b) Expression profiles of the marker genes of all lymphocyte and myeloid subpopulations. (c) UMAP visualization of CTRL, NGAM and GTAM groups. (d) Box plot showing the proportion of CD4 + T cells and Macro-LYVE1 cells among the CTRL ( n  = 3), NGAM ( n  = 3) and GTAM ( n  = 5) groups in scRNA-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (e) Deconvolution analysis using CIBERSORT showing enrichment of CD4 + T cells and Macro-LYVE1 cells in the NGAM group. Sample sizes were as follows: CTRL ( n  = 18), NGAM_eut ( n  = 12), NGAM_ect ( n  = 12), GTAM_eut ( n  = 20), and GTAM_ect ( n  = 20). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (f) Volcano plots showing DEGs in CD4 + T cells between CTRL and NGAM (left), and between GTAM and NGAM samples (right) (FDR  0.5). (g) Venn plot showing that most genes upregulated in the NGAM group were significantly downregulated in the GTAM group (Fisher’s exact test). (h) Box plot showing expression levels of CD44 and TGFB1 among the CTRL ( n  = 18), NGAM ( n  = 24) and GTAM ( n  = 40) groups in Geo-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (i) Reactome pathway analysis based on DEGs in (f, left) indicating T-cell activation in the NGAM group. P values were calculated using the hypergeometric distribution and adjusted using the false discovery rate method. (j) Heatmap showing AUC scores of marker genes related to M1 phenotype, proliferation, angiogenesis and phagocytosis of all the macrophages in scRNA-seq data among the CTRL, NGAM and GTAM groups. (k) GSEA analysis showing that IL-2/STAT5 signaling, TGF-β signaling, inflammatory response, and TNFα signaling via NF-κB pathways were upregulated in the NGAM group and downregulated following GnRHa treatment. P values were adjusted using the Benjamini–Hochberg method. NGAM_eut, eutopic endometrial mesenchyma of NGAM group; NGAM_ect, ectopic endometrial mesenchyma of NGAM group; GTAM_eut, eutopic endometrial mesenchyma of GTAM group; GTAM_ect, ectopic endometrial mesenchyma of GTAM group. To minimize potential sampling biases and further validate our finding, we integrated immune cell data from several published scRNA-seq studies of the healthy uterus during the proliferative phase of the menstrual cycle[ 14 , 15 , 16 , 17 , 18 ] (Supplementary Fig. S7b-d and Sup. Table S3) . Correlation analysis revealed that the immune microenvironment in adenomyosis (NGAM group) differed significantly from that of the control group (Supplementary Fig. S7e) . Minor variation among different control datasets may be attributed to potential batch effect. After GnRHa treatment, the immune microenvironment in patients with adenomyosis closely resembled that of the control group (Supplementary Fig. S7e) . Furthermore, based on the integrated data, we compared CD4 + T cell transcriptomic profiles among the CTRL, NGAM and GTAM groups. DEGs that were upregulated in the NGAM compared to the CTRL group showed a significant decrease following GnRHa treatment (Fisher’s exact test, P < 2.2e -16 ). This finding indicates that most of these genes became downregulated in the GTAM group when compared to the NGAM group, and vice versa ( Fig. 3 , Fig. 3 , Supplementary Fig. S7f and Sup. Table S7) . These results suggest a restoration of the CD4 + T cells immune gene program in adenomyosis following GnRHa treatment. Interpretatively, both the scRNA-seq and Geo-seq data revealed high expression of CD44 , TGFB1 , CXCR4 and IL7R , which have been previously reported to be involved in the activation of CD4 + T cells[ 38 ]. These genes were markedly downregulated following GnRHa treatment ( Fig. 3 , Fig. 3 and Supplementary Fig. S7g) . Enrichment analysis based on upregulated genes indicated significant enrichment in CD28 stimulation and IL4/IL13 signalings pathways, whereas the analysis of downregulated genes showed enrichment in the innate immune system and PD1 immunosuppression pathways, suggesting enhanced CD4 + T cell activation in patients with adenomyosis ( Fig. 3i ) . Macrophages have been recognized as key players in the pathogenesis of endometriosis[ 39 ]; therefore, we sought to investigate their role in adenomyosis. We observed an enrichment of Macro-LYVE1 cells in adenomyosis tissue without GnRHa treatment ( Fig. 3 , Fig. 3 ) . These macrophages that highly express LYVE1 are known to faciliate tissue repair and regeneration by clearing cellular debris, supporting extracellular matrix remodeling, and promoting angiogenesis[ 40 ]. Following GnRHa treatment, the proportion of Macro-LYVE1 cells was significantly reduced ( Fig. 3 , Fig. 3 ) . Consisitent with this finding, analysis of integrated published macrophage datasets revealed that the GTAM group exhibited a markedly lower angiogenesis signature, even lower than that of the CTRL group ( Fig. 3j ) . Moreover, the NGAM group displayed a stronger M1 macrophage signature, indicative of a pro-inflammatory state. DEG and GSEA analysis further revealed enrichment of inflammation-related pathways, such as IL/STAT5, TNFβ, inflammatory and TNFα signaling pathway in adenomyosis tissues ( Fig. 3k ) . Intriguingly, these inflammation-related pathways were significantly suppressed following GnRHa treatment ( Fig. 3k ) . In summary, CD4 + T cells in adenomyosis exhibit high expression of activation-associated genes, where are downregulated after GnRHa treatment, thereby reducing CD4 + T cell activation and potentially improving clinical pregnancy outcomes. Concurrently, GnRHa treatment decreases the proportion of angiogenesis-associated Macro-LYVE1 macrophages, suppressing these processes and potentially reducing disease activity. Collectively, these findings suggest that GnRHa treatment may effectively manage adenomyosis by modulating both CD4 + T cell activation and macrophage-mediated immune responses. Adenomyosis is characterized by the excessive proliferation of fibrous tissue[ 4 ], the presence of endometrial islands embedded within hypertrophic smooth muscle cells, and an invasive process of the endometrial basalis that resembles malignant tumor metastasis. To elucidate the diversity of mesenchymal cells in adenomyosis, we examined their heterogeneity and identified 12 distinct cellular populations. These included post-capillary venous endothelial cells(PCV; enriched for RGCC and SGK1 ), capillary endothelial cells (CapiEC; enriched for ACKR1 , SELP and CCL14 ), artery endothelial cells (ArteryEC; enriched for BMX , SEMA3G , NEBL and IGFBP3 ), lymphatic endothelial cells (LEC; enriched for FABP4 , CCL21 and FABP5 ), fibroblast (enriched for TMEM196 and PLA2G7 ), C7 + fibroblasts (Fibro-C7; enriched for C7 )[ 14 ], MKI67 + fibroblast (Fibro-MKI67; enriched for MMP11 , MKI67 , TOP2A and BIRC5 ), pericytes (enriched for RGS5 ), RGS6 + vascular smooth muscle cells (vSMC; enriched for RGS6 ), MCAM + vSMC (enriched for RERG , MUSTN1 and MCAM ), CLSTN2 + vSMC (enriched for ECRG4 , INPP4B and CLSTN2 ), and DES + SMC (enriched for CLDN1 , GADD45G and DES ) ( Fig. 4 , Fig. 4 ) . No significant batch effect was observed among the groups ( Fig. 4c ) . Fig. 4 Enrichment of PCV and pericytes in adenomyosis. (a) UMAP visualization of all mesenchymal cell subpopulations. (b) Bubble chart displaying marker genes across fibroblasts, SMCs and endothelium cell subpopulations. (c) UMAP visualization of all cells among the CTRL, NGAM and GTAM groups. (d) Box plot showing the proportion of PCV and pericytes among the CTRL ( n  = 3), NGAM ( n  = 3) and GTAM ( n  = 5) groups in scRNA-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (e) Results of e CIBERSORT analysis showing the proportions of PCV and pericytes among the CTRL, NGAM and GTAM groups. Sample sizes for left panel were as follows: CTRL ( n  = 18), NGAM_eut ( n  = 12), NGAM_ect ( n  = 12), GTAM_eut ( n  = 20), and GTAM_ect ( n  = 20). Sample sizes for right panel were as follows: CTRL ( n  = 29), NGAM_JZ ( n  = 12), NGAM_outer ( n  = 12), GTAM_JZ ( n  = 27), and GTAM_outer ( n  = 35). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (f) Volcano plots showing DEGs in pericyte between CTRL and NGAM groups (FDR  0.5). (g) Box plot showing the expression levels of the ADAMTS4 and ADAMTS9 among the CTRL ( n  = 69), NGAM ( n  = 72) and GTAM ( n  = 159) groups in Geo-seq data (Kruskal-Wallis test with post-hoc analysis). (h) GSEA of pericytes in adenomyosis showing higher enrichment of gene programs associated with inflammatory pathways. P values were adjusted using the Benjamini–Hochberg method. (i) Cell-cell communication network between pericytes and other endothelial cells. (j) Regulon analysis using the SCENIC tool showing activation of SOX7 and CREB3L2 regulons and inactivation of CHURC1, MAZ and SMARCC2 regulons in the NGAM group. (k) Feature plots showing increased SOX7 and CREB3L2 regulon activity among the CTRL, NGAM and GTAM groups. (l) Box plot showing the expression levels of the SOX7 and CREB3L2 among the CTRL ( n  = 69), NGAM ( n  = 72) and GTAM ( n  = 159) groups in Geo-seq data (Kruskal-Wallis test with post-hoc analysis). NGAM_eut, eutopic endometrial mesenchyma of NGAM group; NGAM_ect, ectopic endometrial mesenchyma of NGAM group; GTAM_eut, eutopic endometrial mesenchyma of GTAM group; GTAM_ect, ectopic endometrial mesenchyma of GTAM group; NGAM_JZ, junctional zone of NGAM group; NGAM_outer, outer myometrium of NGAM group; GTAM_JZ, junctional zone of GTAM group; GTAM_outer, outer myometrium of GTAM group. Enrichment of PCV and pericytes in adenomyosis. (a) UMAP visualization of all mesenchymal cell subpopulations. (b) Bubble chart displaying marker genes across fibroblasts, SMCs and endothelium cell subpopulations. (c) UMAP visualization of all cells among the CTRL, NGAM and GTAM groups. (d) Box plot showing the proportion of PCV and pericytes among the CTRL ( n  = 3), NGAM ( n  = 3) and GTAM ( n  = 5) groups in scRNA-seq data. P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (e) Results of e CIBERSORT analysis showing the proportions of PCV and pericytes among the CTRL, NGAM and GTAM groups. Sample sizes for left panel were as follows: CTRL ( n  = 18), NGAM_eut ( n  = 12), NGAM_ect ( n  = 12), GTAM_eut ( n  = 20), and GTAM_ect ( n  = 20). Sample sizes for right panel were as follows: CTRL ( n  = 29), NGAM_JZ ( n  = 12), NGAM_outer ( n  = 12), GTAM_JZ ( n  = 27), and GTAM_outer ( n  = 35). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (f) Volcano plots showing DEGs in pericyte between CTRL and NGAM groups (FDR  0.5). (g) Box plot showing the expression levels of the ADAMTS4 and ADAMTS9 among the CTRL ( n  = 69), NGAM ( n  = 72) and GTAM ( n  = 159) groups in Geo-seq data (Kruskal-Wallis test with post-hoc analysis). (h) GSEA of pericytes in adenomyosis showing higher enrichment of gene programs associated with inflammatory pathways. P values were adjusted using the Benjamini–Hochberg method. (i) Cell-cell communication network between pericytes and other endothelial cells. (j) Regulon analysis using the SCENIC tool showing activation of SOX7 and CREB3L2 regulons and inactivation of CHURC1, MAZ and SMARCC2 regulons in the NGAM group. (k) Feature plots showing increased SOX7 and CREB3L2 regulon activity among the CTRL, NGAM and GTAM groups. (l) Box plot showing the expression levels of the SOX7 and CREB3L2 among the CTRL ( n  = 69), NGAM ( n  = 72) and GTAM ( n  = 159) groups in Geo-seq data (Kruskal-Wallis test with post-hoc analysis). NGAM_eut, eutopic endometrial mesenchyma of NGAM group; NGAM_ect, ectopic endometrial mesenchyma of NGAM group; GTAM_eut, eutopic endometrial mesenchyma of GTAM group; GTAM_ect, ectopic endometrial mesenchyma of GTAM group; NGAM_JZ, junctional zone of NGAM group; NGAM_outer, outer myometrium of NGAM group; GTAM_JZ, junctional zone of GTAM group; GTAM_outer, outer myometrium of GTAM group. Next, we analyzed the cell composition of each group and observed a significant enrichment of PCV and pericytes populations in adenomyosis based on scRNA-seq analysis ( Fig. 4d and Supplementary Fig. S8a). This finding suggests the presence of a distinct microenvironment characterized by abundant vascularization in adenomyosis. By integrating the scRNA-seq and Geo-seq datasets, we found that PCV cells predominantly increased in ectopic lesions, whereas pericytes were elevated in both the junctional zone and outer myometrium. Notably, following GnRHa treatment, the proportion of PCV decreased, while the abundance of pericytes remained largely unchanged ( Fig. 4e ) . Pericytes are an essential components of the vascular wall and play critical roles in angiogenesis, vascular stabilization, regulation of permeability, and immune modulation within the microenvironment[ 41 ]. To gain deeper insights into the functional changes of pericytes under different conditions, we performed DEG analysis focused on this cell population ( Fig. 4f and Sup. Table S8) . The results revealed significant upregulation of genes associated with pericyte migration and vascular stabilization, such as ADAMTS4 [ 42 ] and ADAMTS9 [ 43 ], in the NGAM group ( Fig. 4 , Fig. 4 ) , suggesting that pericytes in this group may exhibit enhanced capabilities for extracellular matrix remodeling and vascular support. Furthermore, GSEA indicated that pericytes from the NGAM group showed activation of inflammation-related signaling pathways, including TNFα signaling and inflammatory response, as well as enrichment in proliferation-associated pathways such as G2M checkpoint and MYC targets ( Fig. 4h ) . These findings suggest that pericytes in the NGAM group exist in a more active state, charactericed by heightened inflammatory and proliferative activity, highlighting the potential influcence of microenvironmental alterations on pericyte function. In addition, we observed increased ligand-receptor-mediated cell–cell communication signals in pericyte from the NGAM group compared with those from the CTRL and GTAM groups. Notably, several key ligand–receptor pairs were significantly enriched, including LAMC1 - (ITGA1/2/5/6/9 + ITGB1) , NAMPT-INSR , CCL2-ACKR1 , VEGFA-VEGFR1/2 , and TGFB1-(TGFBR1 + TGFBR2) ( Fig. 4i ) . Consistently, pericytes in the NGAM group exhibited high expression of ligands such as LAMC1 , NAMPT , CCL2 , VEGFA and TGFB1 (Supplementary Fig. S8b) . Previous studies have reported that these ligand-receptor interations, including LAMC1-(ITGA1/2/5/6/9 + ITGB1) [ 44 ], NAMPT-INSR [ 45 ], CCL2-ACKR1 [ 46 ], VEGFA-VEGFR1/2 [ 47 ], and TGFB1-(TGFBR1 + TGFBR2) [ 48 ], are closely associated with tumor invasion and metastasis. Their upregulation in adenomyosis suggests a potential involvement in the disease pathogenesis, possibly through mechanisms related to vascular remodeling, inflammation, and tissue invasion. Finally, we performed transcriptional regulatory networkanalysis using the SCENIC tool[ 20 ]. The results revealed activation of SOX7 and CREB3L2 regulons, along with inactivation of the CHURC1 , MAZ and SMARCC2 regulons in the NGAM group compared with the CTRL and GTAM groups ( Fig. 4j ) . SOX7 , a transcription factor associated with cell proliferation[ 49 ], was primarily activated in the PCV cluster, while CREB3L2 , a transcription factor involved in angiogenesis[ 50 ], was mainly activated in pericytes from patients in the NGAM patients ( Fig. 4k ) , consistent with the Geo-seq data ( Fig. 4l and Supplementary Fig. S8c) . These findings underscore the pivotal roles of pericytes and PCV cells in the pathogenesis of adenomyosis, suggesting that GnRHa treatment may exert its therapeutic effects primarily by inhibiting abnormal angiogenesis. The diagnostic criteria for adenomyosis typically include a thickened junctional zone[ 51 , 52 ]. While existing studies have mainly focused on the morphological alterations of this region, there remains a lack of in-depth molecular characterization. Leveraging the advanced Geo-seq technology, we identified the transcriptional differences between the junctional zone and outer myometrium in both control subjects and patients with adenomyosis. Intriguingly, the outer myometrium was enriched for angiogenesis-related pathway, whereas the junctional zone showed enrichment in pathways related to cell adhesion and response to chemical stimulus ( Fig. 5 , Fig. 5 and Sup. Table S9) . Previous studies describing the pathological changes of the junctional zone have largely remained at the morphological level, providing limited insight into the the underlying molecular mechanisms. To address this gap, we performed differential enrichment analysis comparing the junctional zone between the CTRL and NGAM groups. Notably, the junctional zone of the NGAM group exhibited significant enrichment of EMT, angiogenesis and TNFα signaling pathways (Supplementary Fig. S9a), suggesting their potential involvement in the pathogenesis of adenomyosis. By integrating scRNA-seq and Geo-seq data, we further observed a marked enrichment of mast cells in the NGAM group, particularly within the junctional zone ( Fig. 5c ) . This observation was validated through multiplex immunofluorescence staining, which confirmed higher mast cell density in both the junctional zone and outer myometrium of the NGAM group compared with controls (Supplementary Fig. S9b and 9c) . Mast cells, identified by the expression of CPA3 , were significantly enriched in these regions in NGAM samples compared with both the CTRL and GTAM groups. As key pro-inflammatory cells, mast cells release various inflammatory mediators that enhance the activity of matrix metalloproteinases and fibroblast growth factors[ 53 ], thereby promoting myometrial fibrosis in adenomyosis. Following GnRHa treatment, the proportion of mast cells was markedly reduced, further suggesting that GnRHa may mitigate disease activity by suppressing immune responses in adenomyosis ( Fig. 5c ) . Fig. 5 Changes in the junctional zone of adenomyosis. (a) DEGs analysis between the outer myometrium and the junctional zone based on Geo-seq data(FDR  0.5). (b) GSEA based on DEGs identified in (a). (c) Boxplot showing the proportions of mast cells calculated usingCIBERSORT among the CTRL, NGAM and GTAM groups. Sample sizes for right panel were as follows: CTRL ( n  = 29), NGAM_JZ ( n  = 12), NGAM_outer ( n  = 12), GTAM_JZ ( n  = 27), and GTAM_outer ( n  = 35). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (d) Integration our dataset with other published datasets. (e) UMAP visualization of immune cells acorss different uterine region. (f) Boxplot showing the proportions of mast cells in control myometrium, junctional zone and outer myometrium based on the integrated scRNA-seq datasets (Kruskal-Wallis test). (g) Cell-cell communication strength between mast cells and other subsets. Numbers indicate the interaction strength between each pair of cell types in each group, calculated using the CellChat algorithm. The top three interaction strengths in each group are marked in red, with corresponding clusters labeled. Each cluster corresponds to those shown in (e). (h) Heatmap showing the ligand-receptor interaction strength between mast cells and Marco-IL1B cells. (i) Violin plot showing the expression levels of TGFB1 and ADAM17 in mast cells, and HMOX1 , MT1E and MT2A in Macro-IL1B cells among the control myometrium, junctional zone and outer myometrium, based on Geo-seq data. P values were calculated using one-way ANOVA, followed by Bonferroni correction. CT, myometrium from control patients; JZ, junctional zone of adenomyosis patients. Changes in the junctional zone of adenomyosis. (a) DEGs analysis between the outer myometrium and the junctional zone based on Geo-seq data(FDR  0.5). (b) GSEA based on DEGs identified in (a). (c) Boxplot showing the proportions of mast cells calculated usingCIBERSORT among the CTRL, NGAM and GTAM groups. Sample sizes for right panel were as follows: CTRL ( n  = 29), NGAM_JZ ( n  = 12), NGAM_outer ( n  = 12), GTAM_JZ ( n  = 27), and GTAM_outer ( n  = 35). P values were calculated using the Kruskal–Wallis test, followed by Dunn’s post hoc pairwise comparisons with Bonferroni adjustment. (d) Integration our dataset with other published datasets. (e) UMAP visualization of immune cells acorss different uterine region. (f) Boxplot showing the proportions of mast cells in control myometrium, junctional zone and outer myometrium based on the integrated scRNA-seq datasets (Kruskal-Wallis test). (g) Cell-cell communication strength between mast cells and other subsets. Numbers indicate the interaction strength between each pair of cell types in each group, calculated using the CellChat algorithm. The top three interaction strengths in each group are marked in red, with corresponding clusters labeled. Each cluster corresponds to those shown in (e). (h) Heatmap showing the ligand-receptor interaction strength between mast cells and Marco-IL1B cells. (i) Violin plot showing the expression levels of TGFB1 and ADAM17 in mast cells, and HMOX1 , MT1E and MT2A in Macro-IL1B cells among the control myometrium, junctional zone and outer myometrium, based on Geo-seq data. P values were calculated using one-way ANOVA, followed by Bonferroni correction. CT, myometrium from control patients; JZ, junctional zone of adenomyosis patients. To validate the enrichment of mast cells in the myometrium of adenomyosis, we integrated our scRNA-seq dataset with other published datasets[ 8 ] collected from various uterine regions, including the junctional zone and outer myometrium ( Fig. 5 , Fig. 5 , and Supplementary Fig. S10a) . Consistent with our Geo-seq data ( Fig. 5c ) , we observed a higher tendency of abundance of mast cells in the junctional zone compared with the control myometrium ( Fig. 5f ) . Next, we analyzed ligand–receptor–mediated cell–cell communication using scRNA-seq data to assess the interaction strength between mast cells and other cell clusters across the three groups ( Fig. 5g ) . In the control and outer myometrium, mast cells primarily interacted with CD8 + T cells (cluster C15 and C16) and NK cells (cluster C18 and C19) ( Fig. 5g ) . However, in the junctional zone, mast cells exhibited stronger interactions not only with lymphoid cells but also with myeloid cells. The top three interacting clusters in the junctional zone were CD8-IL7R (cluster C16), Macro-IL1B (cluster C25) and monocytes (cluster C28). We specifically focused on macrophages expressing the pro-inflammatory gene IL1B and used the NicheNet tool[ 27 ] to predict ligand-receptor interactions that could drive gene expression changes in both mast cells and Macro-IL1B cells. Through this analysis, we identified potential ligands and active target genes ( Fig. 5h ) , and further examined the expression of these ligands in mast cells and the target genes in Macro-IL1B cells. Our findings revealed that TGFB1 and ADAM17 were highly expressed in mast cells, whereas HMOX1 , MT1E and MT2A were predominantly expressed in Macro-IL1B cells within the junctional zone ( Fig. 5i ). Notabley, a strong regulatory potential was observed between TGFB1 and HMOX1/MT2A ( Fig. 5h ) . To validate these findings, we performed multiplex immunofluorescence staining to assess the co-expression of HMOX1 and IL-1β-positive macrophages in the junctional zone of the NGAM group (Supplementary Fig. S10b and 10c) . Strikingly, cells co-expressing CD68 and IL-1β, representing the Macro-IL1B population, exhibited high levels of HMOX1 in the junctional zone of NGAM samples. Together, these findings suggest that the molecular mechanism underlying the thickening of the junctional zone may involve an inflammatory response, primarily characterized by an increase in mast cells and interactions with Macro-IL1B cells through TGFB1 - HMOX1/MT2A interactions. GnRHa is widely applied in the treatment of adenomyosis[ 54 ]. However, studies investigating the molecular mechanisms underlying the effects of GnRHa in adenomyosis remain limited[ 55 ]. Our above results demonstrate the effects of GnRHa treatment on immune inflammation and angiogenesis across the entire uterus layer ( Fig. 3 , Fig. 4 , Fig. 5 ) . Here, we compared cellular interactions and transcriptional activity between GTAM and NGAM groups to provide further molecular insights supportting the clinical application of GnRHa treatment. Through DEG analysis based across all cells, we confirmed severaldifferentially upregulated genes, such as WFDC2 , PTGDS and SCGB1D2 in GTAM group, which are known to be influenced by hormones[ 56 ] ( Fig. 6a and Sup. Table S10) . The upregulated genes in the GTAM group were primarily enriched in pathways related to oxidative phosphorylation and female pregnancy, consistent with the increased clinical pregnancy rates observed in adenomyosis patients following GnRHa treament[ 57 ] ( Fig. 6b ) . Specifically, the expression of CRIP1 , GNAS , and NDUFA13 , which were highly expressed in the GTAM group according to scRNA-seq, was further confirmed by the Geo-seq data ( Fig. 6 , Fig. 6 ) . CRIP1 primarily functions in oestrogen-sensitive tissues[ 58 ], and GNAS plays a crucial role in regulating PI3K/AKT signaling[ 59 ]. These results suggest that GnRHa treatment may modulate oestrogen sensitivity, thereby promoting embryo implantation. NDUFA13 , a mitochondrial gene involved in cell metabolism[ 60 ], is critical for both embryo implantation and the maintenance of early pregnancy. Fig. 6 The effect of GnRHa in adenomyosis. (a) DEGs analysis based on all cells between NGAM and GTAM groups in scRNA-seq data (FDR  0.5). (b) GO enrichment analysis based on DEGs identified in (a). (c) Spatial expression patterns of CRIP1 , GNAS and NDUFA13 in NGAM and GTAM groups. (d) Boxplot showing expression levels of CRIP1 , GNAS and NDUFA13 in NGAM ( n  = 72) and GTAM ( n  = 159) groups. P values were calculated by a two-tailed Mann–Whitney U test. (e) Number of inferred interactions between NGAM and GTAM groups determined by CellChat analysis. (f) Signaling pathways showing significant differences between NGAM and GTAM groups. (g) Intercellular communication networks of the GRN signaling pathway between NGAM and GTAM groups. (h) Violin plot showing expression levels of GRN between NGAM and GTAM groups based on scRNA-seq data. P value was calculated by a two-tailed Mann–Whitney U test. (i) Heatmap displaying transcription factors differentiall expressed between NGAM and GTAM groups. FOSL1 , FOSL2 , and RELB were significantly downregulated in the GTAM group. (j) Boxplot showing expression levels of FOSL1 , FOSL2 , and RELB in NGAM ( n  = 72) and GTAM ( n  = 159) groups based on Geo-seq data. P value was calculated by a two-tailed Mann–Whitney U test. The effect of GnRHa in adenomyosis. (a) DEGs analysis based on all cells between NGAM and GTAM groups in scRNA-seq data (FDR  0.5). (b) GO enrichment analysis based on DEGs identified in (a). (c) Spatial expression patterns of CRIP1 , GNAS and NDUFA13 in NGAM and GTAM groups. (d) Boxplot showing expression levels of CRIP1 , GNAS and NDUFA13 in NGAM ( n  = 72) and GTAM ( n  = 159) groups. P values were calculated by a two-tailed Mann–Whitney U test. (e) Number of inferred interactions between NGAM and GTAM groups determined by CellChat analysis. (f) Signaling pathways showing significant differences between NGAM and GTAM groups. (g) Intercellular communication networks of the GRN signaling pathway between NGAM and GTAM groups. (h) Violin plot showing expression levels of GRN between NGAM and GTAM groups based on scRNA-seq data. P value was calculated by a two-tailed Mann–Whitney U test. (i) Heatmap displaying transcription factors differentiall expressed between NGAM and GTAM groups. FOSL1 , FOSL2 , and RELB were significantly downregulated in the GTAM group. (j) Boxplot showing expression levels of FOSL1 , FOSL2 , and RELB in NGAM ( n  = 72) and GTAM ( n  = 159) groups based on Geo-seq data. P value was calculated by a two-tailed Mann–Whitney U test. In the cellular interaction analysis, we observed that cell–cell communications became stronger after GnRHa treatment, both in terms of the number and strength of interactions ( Fig. 6e ) . Interaction pathways such as SPP1, SELE and ANGPT were upregulated in the NGAM group, whereas SEMA3, CD22, THY1, GALECTIN and GRN were enriched in the GTAM group ( Fig. 6f ) . GRN , which is involved in cell proliferation and tissue repair[ 61 ], showed extensive intercellular communication through its signaling pathway among multiple cell types in the GTAM group ( Fig. 6g ) . This was further supported by the high expression of GRN observed in the GTAM group through scRNA-seq data ( Fig. 6h ), which may be associated with enhanced clinical pregnancy rates following GnRHa treatment. Finally, we analyzed the impact of GnRHa treatment on transcriptional activity. By comparing the GTAM and NGAM groups (comprising six patients, three from our study and three from Chen et al. [ 8 ]), we identified several significantly different transcriptional profiles. Notably, the regulon activity of FOSL1 , FOSL2 and RELB was significantly reduced after GnRHa treatment, consistent with their lower expression levels in the GTAM group as observed in Geo-seq ( Fig. 6 , Fig. 6 ) . FOSL1 and FOSL2 are members of the activator protein-1 (AP-1) family, while RELB belongs to the nuclear factor-kappa B (NF-κB) family. Both AP-1 and NF-κB are crucial in initiating inflammatory responses by inducing the expression and secretion of chemokines and cytokines that attract and activate immune cells[ 62 ]. Altogether, thses findings suggest that the mechanism by which GnRHa treatment affects adenomyosis involves the suppression of immune responses and activation of granulin pathway. Future validation experiments will be necessary to confirm these findings and further elucidate the underlying mechanisms.

This study was approved by the Ethics Committee of the Institutional Review Board at the Sixth Affiliated Hospital of Sun Yat-sun University and was conducted in accordance with all relevant ethical regulations concerning research involving human participants. Written informed consents were obtained from all participants. Our study enrolled a total of 28 participants, comprising six individuals serving as controls and 22 patients with adenomyosis recruited from Sixth Affiliated Hospital of Sun Yat-sun University in China. The pathogenesis of adenomyosis, characterized by abnormal proliferation, migration, and invasiveness of endometrial cells, was investigated during the proliferative phase of the menstrual cycle, when active cell growth coincides with dysregulated disease-related biological processes. This design aimed to capture relevant transcriptional changes and analyze aberrant cellular behaviors underlying disease progression. Therefore, we included all the participants with proliferative phase. The control group consisted of six subjects in the proliferative phase who underwent hysterectomy for benign conditions unrelated to adenomyosis, such as symptomatic uterine leiomyoma or high-grade squamous intraepithelial lesions. Adenomyosis was diagnosed based on both imaging examinations and histopathological evaluation of hematoxylin and eosin (H&E)-stained tissue sections. Among the 22 adenomyosis patients, eight did not receive any hormonal treatment prior to surgery, whereas the remaining patients were treated with a pre-surgical therapeutic regimen of GnRHa without oestrogen add-back therapy. This pharmacological approach aimed to suppress ovulation and induce a hypoestrogenic state to alleviate symptoms such as excessive menstrual bleeding and dysmenorrhoea. Although this intervention was designed to improve clinical outcomes, the responses varied among patients. A subset of the cohort exhibited suboptimal efficacyor experienced adverse effects associated with reduced oestrogen levels. When conservative medical management failed to achieve satisfactory results, hysterectomy was performed. None of the patients included in the study received any other hormonal, antibiotic, or anti-inflammatory medications within three months prior to surgery. During tissue sampling, the entire uterine layer, from the endometrium to the perimetrium was collected from patients who underwent total hysterectomy. Any macroscopic fibroids or endometrial polyps were carefully avoid, maintaining a minimum distance of one centimetre from these regions. The uterine tissue was then evenly divided into three portions: one for scRNA-seq, one for Geo-seq, and one for paraffin embedding. Lastly, we excluded patients whose H&E-stained slides revealed any microscopic fibroids or endometrial polyp lesions to eliminate potential confounding factors and to ensure the reliability of comparative analyses. scRNA-seq was performed on surgical biopsy specimens from 11 individuals, including three without adenomyosis (CTRL group) and eight with adenomyosis (no GnRHa administration, NGAM group, n  = 3; GnRHa administration, GTAM group, n  = 5). Geo-seq was conducted on tissue samples from 12 individuals, comprising four without adenomyosis (CTRL group) and eight with adenomyosis (NGAM group, n  = 2; GTAM group, n  = 6). Additionally, multiplexed immunofluorescence and immunohistochemistry analyse was performed on tissues from six control individuals and 16 with adenomyosis for further validation. Sup. Table S1 summarizes the treatment parameters, including GnRHa administration (Leuprorelin, initiated between days 1 and 5 of the menstrual period), dosage (3.75 mg), number of administrations (ranging from one to three), the interval between surgery and the last menstrual period (Days since LMP), as well as demographic data for all participants, providing a comprehensive overview of the clinical context. For molecular and immunofluorescence analysis, tissue histology was first confirmed using formalin-fixed, paraffin-embedded (FFPE) H&E-stained slides. Fresh biopsy specimens were preserved in MACS Tissue Storage Solution (#130–100-008, Miltenyi Biotec) and kept on ice until further processing. Uterine tissues from patients with adenomyosis and control subjects were embedded in optimum cutting temperature (OCT) compound, snap-frozen in liquid nitrogen, and labeled for orientation. After the OCT blocks were fully frozen, the tissues were cryosectioned at a thickness of 20 μm. The Geo-seq procedure was performed as described previously[ 9 ]. Briefly, the cells of interest are harvested and lysed for first-strand cDNA synthesis, and a modified Smart2-seq protocol[ 10 ] is used to amplify the cDNA. cDNA samples that pass quality control are then used for library construction and next-generation sequencing on the BGISEQ-2000 (Beijing Genomics Institute) platform using a 150-bp paired-end mode. From the control samples, endometrial gland, endometrial mesenchyma, myometrium, and myometrial mesenchyma were collected. From the adenomyosis samples, eutopic endometrial gland, eutopic endometrial mesenchyma, adenomyosis junctional zone, ectopic endometrial gland, ectopic endometrial mesenchyma, outer myometrium, and outer myometrial mesenchyma were captured by laser microdissection. Multiplexed immunofluorescence was performed on 4 μm-thick FFPE whole-tissue sections using standard primary antibodies sequentially paired with the TSA 2-color kit (abs50012-20 T, Absinbio, China), followed by nuclear counterstaining with DAPI. For instance, deparaffinized slides were first incubated with anti-DNAH9 (#11210R, Yaji Biotechnology) overnight at 4 °C, and subsequently treated with horseradish peroxidase-conjugated (HRP)-conjugated anti-rabbit or anti-mouse secondary antibody (#A0011-60, Absinbio) for 10 min. Signal development was achieved using TSA fluorophore 570 (#323272, Bio-Techne, China) with a strictly observed 10 min. After washing in TBST buffer, slides were transferred to preheated citrate solution (90 °C) and subjected to heat treatment in a microwave at 20 % of the maximum power for 15 min. Slides were then cooled in the same solution to room temperature. Between each step,the slides were thoroughly washed with Tris buffer. Detailed information about all antibodies used is provided in Sup. Table S2 . Following staining, slides were incubated with DAPI (D1306, Thermo Fisher Scientific), rinsed with distilled water, and manually cover-slipped. After air-drying, fluorescence images were captured using an Aperio Versa 8 tissue imaging system (Leica) and analyzed with Halo software (Indica labs, Version 3.3). Immunohistochemistry was performed on 4 μm-thick FFPE whole-tissue sections. After deparaffinization, antigen retrieval was performed in a microwave oven using EDTA buffer (pH 9.0). Endogenous peroxidase activity was blocked with 3 % H 2 O 2 for 20 min. The sections were then incubated with the following primary antibodies: E-cadherin (1:500, GB12083, Servicebio), Vimentin (1:500, 60330–1-lg, Proteintech), β-Catenin (1:500, GB15016, Servicebio). After incubation with a HRP–conjugated secondary antibody (#G1301, Servicebio), the signals were visualized using diaminobenzene (#G1212, Servicebio) and counterstained with haematoxylin. After air-drying, images were captured using a Pannoramic MIDI III Digital Scanner (3DHISTECH) and analysed with ImageJ software (Version 1.54i). Fresh uterine tissues were dissociated into single cells using a dissociation solution (0.35 % collagenase IV, 5.2 mg/mL papain, 120 U/mL DNase I). The digestion was performed in a 37℃ water bath with shaking at 100 rpm for 20 min, and was terminated by adding phosphate-buffered saline containing 10 % fetal bovine serum. The resulting cell suspension was passed through 70 μm nylon cell strainer and centrifuged at 300g for 5 min at 4℃. Red blood cells were removed using a Red Blood Cell Lysis Solution (#130–094-183, Miltenyi Biotec), and dead cells were removed with the Miltenyi ® Dead Cell Removal Kit (#130–090-101, Miltenyi Biotec). Viable cells were counted using Countess II Automated Cell Counters (Thermo Fisher Scientific), and the cell concentration was adjusted to 700–1200 cells/μL. Single-cell suspensions were loaded onto the 10 × Genomics Chromium platform using the Single-Cell 3′ kit (V3, 10 × Genomics) following the manufacturer’s instructions. Library construction was performed according to the standard protocol, and sequencing was carried out on an Illumia NovaSeq 6000 sequencer (Illumina) at a depth of 20,000 reads per cell. The data processing workflow is shown in Supplementary Fig. S1a . FASTQ files were processed using the Cell Ranger count pipeline (v.7.0.0, 10 × Genomics) with default parameters to demultiplex reads, align them to the human reference genome (GRCh38), and perform filtering and quantification, thereby generating feature-barcode matrices. The resulting count matrices were further processed using the Seurat package (v.4.3.0) to exclude low-quality cells through the following steps. First, We loaded the count matrices into Seurat objects using the “Read10X” function and included features detected in at least three cells (min.cells = 3) for each sample. Second, to remove cells with extensive mitochondrial contamination, mitochondrial gene expression was fitted to a normal distribution, and a false discovery rate (FDR) threshold of 10.9 %) were subsequently excluded ( Supplementary Fig. S1b ). Third, to remove low-quality cells or empty droplets, we visualized the number of features (nFeature_RNA) and the number of counts (nCount_RNA) using a scatter plot to determine appropriate thresholds for quality control. Cells with fewer than 360 detected features were excluded, while no threshold was applied to nCount_RNA, as no cells exhibited abnormally low total counts ( Supplementary Fig. S1c ). Finally, to remove cell doublets or multiplets, which may exhibit an aberrantly high gene count, we employed the DoubleFinder package (v.2.0.3) and applied a cutoff at the 92.5th percentile of the doublet score. Cells exceeding this threshold were identified as potential doublets and subsequently excluded from downstream analysis. Following these quality control steps, we obtained a final dataset comprising 100,404 single cells for subsequent analyses. To integrate individual samples and identify common sources of variation, we performed canonical correlation analysis (CCA) and evaluated the results using the adjusted rand index and normalized mutual information[ 11 ], two commonly used metrics for assesing single-cell clustering performance. The data were then normalized using a global-scaling method with a scaling factor of 10,000 to partially mitigate technical effects, implemented via the“NormalizeData” function. After normalization, a pseudo-count of one was added to each value followed by log-transformation. Highly variable genes (top 2,000) were identified using the “FindVariableFeatures” function. Subsequently, the gene expression matrix was scaled and centered using the “ScaleData” function. For dimensionality reduction and visualization, we performed principal component analysis (PCA) using the “RunPCA” function based on the highly variable 2,000 genes. The top 30 principal components were then used for uniform manifold approximation and projection (UMAP) analysis to visualize the clusters in two-dimensional space. Finally, clustering was performed using the Leiden community detection algorithm to identify distinct cell populations based on their transcriptional profiles. We conducted a series of optimization trials for each major cell type, testing resolutions between 0.4 and 1.4 and employing correlation analysis to assess cluster separation. The final resolutions used for epithelial cells, immune cells, fibroblasts, and smooth muscle cells (SMC) were 0.6, 0.8, 0.6, and 0.8, respectively. The assignment of major cell types to each identified cluster was performed through a rigorous two-step approach. Initially, marker genes for each cluster were identified using Wilcoxon rank-sum test in a one-versus-all manner with the “FindAllMarkers” function with a minimum fold-change of 0.5–2 between groups. By comparing each individual cluster against all others, we identified genes exhibiting significant differential expression. Subsequently, cell types were assigned to the respective clusters by examing the top-ranked differentially expressed genes for the presence of well-characterized, cell-type-specific markers. Referring to pulished literature and the established marker gene databases CellMarker 2.0[ 12 ], we inferred the most probable cell type for each cluster based on these differentially expressed markers. For the annotation of immune–cell subcluster, we employed an unbiased cell type recognition method, SingleR (v.2.0.0)[ 13 ], which leverages reference transcriptomic datasets of known cell types for annotation. SingleR was first applied to verify the consistency of the predicted annotations for each query cluster across multiple reference transcriptomic datasets integrated within the SingleR framework. Upon confirming this consistency, the corresponding cell type annotations were confidently assigned to each cluster. We collected published datasets[ 14 , 15 , 16 , 17 , 18 ] and integrated them with our scRNA-seq dataset for comparative analysis. Analysis, including preprocessing, integration, normalization, clustering, and visualization was performed using methods consistent with those described above. Notably, for datasets derived from healthy uterine tissue, only immune cells were extracted and integrated with our dataset. The correlation between samples was assessed using Spearman correlation analysis. The Kolmogorov-Smirnov statistic, computed using the “ks.test” function, was used to quantify differences among the control group, the adenomyosis group (without complications), and the complications groups (adenomyosis patients with other complications). The primary differences among the three groups were observed between the adenomyosis and control groups. To further investigate differences between the adenomyosis and complications groups, we employed CINEMA-OT[ 19 ] to evaluate the impact of complications on adenomyosis. By calculating the Silhouette score (using the “silhouette_samples” function in sklearn.metrics) for clusters corrected by CINEMA-OT, we found that the fibroblasts exhibited the greatest divergence between the two groups. Differential expression analyse comparing NGAM vs CTRL and GTAM vs NGAM was performed on scRNA-seq data using the presto (v.1.0.0) package with default parameters. P values were computed with the Mann–Whitney U test and adjusted using the Benjamini–Hochberg method. Genes with an adjusted P  0.5 were considered significantly differentially expressed. The raw reads from Geo-seq underwent a series of pre-processing steps. First, quality control was performed using FastQC (v.0.11.9) to assess read quality. Adapter sequences were trimmed sequentially using Trim Galore (v.0.6.7). The resulting clean reads were aligned to the human reference genome (GRCh38) using HISAT2 (v.2.2.1). Samtools (v.1.17) was used to index and merge BAM files. Transcript expression levels were estimated with StringTie (v.2.1.7), which utilizes GTF annotation and BAM alignment to quantifytranscript abundance. The prepde.py3 script from StringTie was used to extract count dataand calculate transcripts per million (TPM) values for each gene. TPM normalization accounts for both gene length and read count, providing a standardized measure of gene expression. Differentially genes expression (DGE) between patients with and without adenomyosis was performed on the Geo-seq data using the edgeR package (v.3.40.2) with default parameters. Differentially expressed genes were defined as those with an adjusted P -value  0.5. To evaluate the adenomyosis-associated gene programs, we performed gene set enrichment analysis. The “enrichGO” function from the clusterProfiler package (v.4.6.2) was used to assess the fold enrichment score and P value of differentially expressed genes within various Gene Ontology (GO) categories, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. P values were calculated using the hypergeometric test and adjusted with the Benjamini–Hochberg method. An adjusted P < 0.05 was considered statistically significant. This analysis allowed us to determine whether specific gene sets or pathways were significantly enriched among the differentially expressed genes. In addition, we performed fast gene set enrichment analysis using fGSEA package (v.1.11.1), focusing on identifying common pathways using the hallmark gene sets in the Molecular Signatures Database (MsigDB). P values for the GSEA analysis were adjusted using Benjamini–Hochberg method and adjusted P < 0.05 was regarded as statistically significance, indicating the enrichment of specific gene sets or pathways associated with adenomyosis. To evaluate the enrichment of specific gene sets or other relevant gene programs among the highly expressed genes in each cell, we employed the AUCell package (v.1.30.1)[ 20 ], which is specifically designed for scRNA-seq data analysis. AUCell identifies cells with active gene regulatory networks by quantifying the enrichment of a gene set within individual cells. It calculates a gene-set activity score by measuring the area under the recovery curve (AUC) across the ranked expression values of top expressed genes (5 % in default) within each cell. This score reflects both the proportion of expressed genes in the query gene set and their relative expression levels compared to other genes in the same cell. Notably, this approach is independent of gene-expression units and normalization methods, providing a robust and unbiased measure of gene-set activity. We focused on assessing the activity scores of gene sets associated with epithelial-mesenchymal transition, endothelial-mesenchymal transition, and the WNT signaling pathway in each cell. To compare the cellular profile of patients with adenomyosis and controls in each Geo-seq sample, we applied the CIBERSORT[ 21 ] (v.1.03) and spatialDWLS[ 22 ] (v.1.0.2) algorithms, which estimate cell-type composition at each tissue location using Geo-seq transcriptomics data and scRNA-seq signature matrix. For the CIBERSORT algorithms, the scRNA-seq signature matrix (denoted as B) was generated using the “FindAllMarker” function in the Seurat package and served as a reference. It comprised K cell types, each characterized by a set of G cell type-informative genes. In the Geo-seq transcriptomic data, we denote the G-by-N gene expression matrix as X, where G represents the same set of informative genes and N represents the number of spatial locations. The N-by-K cell-type composition matrix was denoted as V, wich each row of V representing the proportions of the K cell types within a spatial location. To estimate V, the CIBERSORT algorithm employs a non-negative matrix factorization model linking the three matrices: X = B V T + E where E is a G-by-N residual error matrix. The CIBERSORT algorithm applied aν-support vector regression (ν-SVR)-based machine learning method to estimate matrix V, which has been shown to outperform other approaches[ 21 ]. This analysis provided spatially resolved insights into the cellular heterogeneity and potential differences in cell-type proportions between the two groups, thereby contributing to a better understanding of the cellular dynamics underlying adenomyosis. To validate the deconvolution results obtained by CIBERSORT, we further applied the spatialDWLS algorithm using the default parameters in the official tutorials[ 22 ]. Transcription factor activity was predicted using pySCENIC (v.0.12.1) with default parameters[ 23 ]. In brief, co-expression modules were identified from the normalized RNA-seq expression matrix. The correlations between transcription factors and their target genes were determined using GRNBoost2[ 24 ]. Regulons were subsequently inferred based on publicly available motif-binding databases. scRNA-seq data of endometriosis was downloaded from NCBI Gene Expression Omnibus (GEO) under accession number GSE213216 [ 25 ]. To assign cell types within this dataset, we used the markers EPCAM and CDH1 to identify epithelial cells. CCA was performed to identify shared sources of variation between datasets. To ensure inclusion of only relevant clusters, we excluded any clusters that did not express marker genes associated with canonical keratins. This refinement allowed us to focus on clusters exhibiting epithelial characteristics. Finally, ciliated epithelium cells were identified using well-established markers such as FOXJ1 , PIFO and TPPP3 . Cell-cell interaction analysis was performed using CellChat(v.1.6.1) [ 26 ]. To ensure data integrity, we performed quality inspection and normalization on both the control and adenomyosis datasets. Following the official workflow and default parameter settings, the two datasets were independently loaded into CellChat. Potential ligand-receptor interactions between cell types were identified using the built-in CellChatDB.human database, which contains curated information on known human receptor-ligand interactions. By leveraging this reference, we screened for specific interactions and evaluated their potential biological significance in the context of adenomyosis. To quantify the strength and significance of ligand-receptor interactions, several CellChat functions were employed. The “computeCommunProb” function was used to calculate communication probabilities between cell types, while the “computeCommunProbPathway” function identified potential signaling pathways associated with the interactions. The “aggregateNet” function was then used to integrate and visualize the resulting interaction network. To further explore ligand-receptor interactions that may contribute to transcriptional changes within cell clusters of interest, we applied NicheNet(v.2.2.0) [ 27 ] to the single-cell matrix following the official tutorials. Statistical analysis was performed using R software (v.4.2.2). For comparisons between two groups, a two-tailed Mann–Whitney U test was applied. For comparisons among multiple groups, the Kruskal–Wallis test (the nonparametric version of one-way ANOVA) was used unless othervise specified. When the Kruskal-Wallis test indicated a significant difference between groups ( P < 0.05), post-hoc pairwise Wilcoxon rank-sum tests were conducted, followed by Bonferroni correction. Differentially expressed genes in the scRNA-seq data were identified using the “FindAllMarkers” function in the Seurat packages (v.4.3.0).

All procedures performed in this study were conducted in accordance with the ethical standards of the institutional and national research committees and with the Helsinki Declaration of 1975, as revised in 2008 (5). Informed consent was obtained from all patients included in the study. The study was approved by the Ethics Committee of the Sixth Affiliated Hospital, Sun Yat-sen University. The ethical approval number is 2020ZSLYEC-281.

Studies on the pathogenesis of adenomyosis are constantly evolving. However, due to the challenges in obtaining clinical specimens and the difficulties in in vitro modeling, the exact mechanisms remain elusive. Previous researches on the endometrium of women with adenomyosis have typically focused on the expression of a limited set of genes[ 63 ] or included small patient cohorts[ 7 , 8 ]. To address this knowledge gap, our study employed scRNA-seq to profile the entire uterine layer, providing a more comprehensive transcriptional analysis of adenomyosis. Our dataset includes a relatively larger sample size, and importantly, it also includs individuals who had received GnRHa pre-surgical treatment. Thus, this study offers novel insights into the cellular mechanism underlying adenomyosis and the effect of GnRHa treatment on the condition. Epithelial cells play a crucial role in the function and pathology of the endometrium. scRNA-seq data analysis of adenomyosis conducted by Liu et al. [ 64 ] revealed an inflammation feature in the endometrium of adenomyosis, consistent with our findings. This observation also aligns with the high comorbidity rate of adenomyosis and endometritis, which often leads to recurrent miscarriages. Additionally, our integrated spatial analysis identified a higher proportion of ciliated epithelial cells in ectopic endometrial glands, regardless of GnRHa administration. Previous studies have also identified ciliated cells in endometriosis, showing distinct transcriptome profiles compared to other epithelial cells and exhibiting high activity of WNT pathway targets[ 14 , 16 , 25 ]. Our study confirms the presence of ciliated cells in both in situ and ectopic endometrium, but with significantly higher proportions observed in ectopic endometrial glands, as validated by immunofluorescence. Moreover, DEGs in epithelial cells from adenomyosis are significantly enriched in pathways related to ciliated cells movement, suggesting that ciliated cells may play an important role in the proliferation of ectopic endometrium. Given that adenomyosis is a hormone-dependent disease and oestrogen (E2) signaling has been shown to drive ciliogenesis in endometrial organoid[ 65 ], we speculate that ciliated cells may proliferate in ectopic lesions and contribute to the proliferation of ectopic endometrium. Experimental evidence elucidating the precise mechanisms through which ciliated epithelial cells contribute to the pathogenesis of adenomyosis is currently lacking. As a next step, we plan to incorporate in vivo mouse models and in vitro organoid systems to conduct functional validation studies. Imbalanced inflammation and apoptotic signals have been reported in adenomyosis[ 31 , 66 ]. The ectopic endometrium may activate inflammatory and apoptotic responses, potentially leading to infertility in patients with adenomyosis. In addition to the inflammatory status of the immune microenvironment, GnRHa effectively targets inflammatory cells, such as CD4 + T cells and Macro-LYVE1 cells, thereby alleviating symptoms in patients with adenomyosis. The immune system plays a vital role in maintaining the uterine cavity microenvironment. Macrophages and T cells are modulated by this microenvironment and, in turn, contribute to tissue homeostasis[ 67 ]. Macrophages are key players in both innate and adaptive immune responses, as they are involved in angiogenesis, antigen presentation to T cells, and the secretion of cytokines/chemokines. In the normal endometrium, macrophages represent approximately 10 % of the total endometrial immune cell population[ 68 ]. While a recent study suggested that macrophages present an immunotolerant phenotype in the endometriosis microenvironment[ 32 ], few studies have focused on the macrophage function within the context of the adenomyosis myometrial microenvironment. In the NGAM group, we observed higher M1 phenotype and angiogenesis scores of macrophages compared to the CTRL group. Several lines of evidence suggest that macrophages are not only “trapped” in ectopic endometrial lesions but also locally activated[ 69 ]. Lousse et al. confirmed that the NF-κB transcription factor is activated in macrophages from endometriotic patients, and responsive genes that regulate macrophage functions undergo transactivation, supporting angiogenesis and tissue remodeling[ 70 ]. We also found higher phagocytosis and proliferation scores of macrophages in the GTAM group, suggesting that GnRHa treatment may reduce ectopic lesions by activating macrophage function. Further studies are needed to fully elucidate the role of macrophages in adenomyosis. The integrity of the uterine junctional zone is crucial for preventing the eutopic endometrium from invading the myometrium[ 51 ]. Multiple pregnancies and iatrogenic abortions may damage this junctional zone, predisposing individuals to the development of endometriosis and adenomyosis[ 71 ]. Inflammation at the endo -myometrial interface may further facilitate endometrial invasion into the myometrium, contributing to the formation of adenomyosis lesions[ 35 ]. Taking advantages of Geo-seq technology, we conducted a detailed transcriptional analysis of the junctional zone. Comparative profiling between the junctional zone and the outer mayometrium revealed a marked enrichment of cell adhesion-related pathways in the junctional zone, supporting the hypothesis that endometrial invagination into the myometrium represents a key pathogenic mechanism in adenomyosis. Structural disruption of the junctional zone may therefore predispose to ectopic implantation of eutopic endometrial tissue, driving disease progression. Notably, the junctional zone of the NGAM group exhibited pronounced enrichment of EMT, angiogenesis, and TNFα via NF-κB signaling pathways compared with healthy controls. These findings further reinforce the hypothesis that chronic inflammation at the endo -myometrial interface plays a critical role in the pathogenesis of adenomyosis. Clinically, patients with adenomyosis exhibit a significantly thickened junctional zone compared with healthy controls. Moreover, cycle-dependent contraction waves originating from the junctional zone are closely associated with the pathological symptoms of adenomyosis [ 51 ]. In this study, we sought to elucidate the underlying mechanisms of junctional zone involvement in adenomyosis. We observed a notable tendency for mast cell accumulation within the junctional zone compared with the outer myometrium in the NGAM group. This trend was further supported by the integration of published datasets[ 8 ]. Additionally, our analyses revealed an enrichment of mast cells in the junctional zone and their interactions with the pro-inflammatory Macro-IL1B macrophages. EMT has been associated with myocyte fibrosis[ 72 ], and fibrosis, angiogenesis, and activation of the TNFα pathway may represent key molecular mechanisms underlying junctional zone thickening. Beyond these processes, mast cells appear to play an essential role. They can release histamine, nerve growth factors and prostaglandins, substances known to stimulate nerve fiber proliferation and increase neural sensitivity[ 73 ]. Thus, the accumulation of mast cells in the junctional zone of adenomyosis patients may contribute to abnormal contractions. Interestingly, following GnRHa treatment, the observed reduction in mast cells within the junctional zone may partially explain the improved clinical pregnancy rates among patients with adenomyosis. Together, these findings underscore the pivotal role of mast cells in the pathophysiology of adenomyosis. Although GnRHa is widely used for the treatment of adenomyosis, most studies have focused on its role in suppressing oestrogen and progesterone levels[ 5 ], suggesting that its therapeutic effect primarily arises from hormonal suppression. However, GnRHa has also been shown to improve clinical pregnancy rates, and its underlying molecular mechanism remain unclear. In this study, by integrating scRNA-seq and Geo-seq data, we identified molecular mechanisms through which GnRHa may enhance clinical pregnancy rates, involving diverse cell–cell interactions. CRIP1 , GNAS and NDUFA13 were highly expressed in the GTAM group. CRIP1 primarily functions in oestrogen-sensitive tissues[ 58 ], and GNAS plays a crucial role in regulating PI3K/AKT signaling[ 59 ], suggesting that GnRHa may modulate oestrogen sensitivity and thereby enhance clinical pregnancy outcomes. NDUFA13 , a mitochondrial gene involved in cell metabolism[ 60 ], is critical for both embryo implantation and the maintenance of early pregnancy. Cell-cell communication analysis revealed stronger interactions involving SELE and ANGPT in the NGAM group, whereas SEMA3 , CD22 , THY1 , GALECTIN and GRN were enriched in the GTAM group. Granulin signaling promotes cell proliferation and tissue repair[ 74 ]. We found that GRN was highly expressed in the GTAM group, suggesting that GnRHa may influence clinical pregnancy outcomes through multiple cell–cell interactions, including granulin signaling. However, further functional assays are required to validate these findings and elucidate the precise mechanisms underlying the therapeutic effects of GnRHa in adenomyosis. Adenomyosis is a complex pathology driven by intricate molecular mechanisms. In this study, we presented a comprehensive cellular landscape of adenomyosis by integration Geo-seq and scRNA-seq technologies, providing a valuable framework for dissecting molecular changes in tissue architecture. In summary, adenomyosis may result from multiple pregnancies or iatrogenic factors that damage the junctional zone, leading to the invasion of ciliated cells into the myometrium. The ectopic endometrium can subsequently trigger immune response involving CD4 + T cells and macrophages, while also promoting angiogenesis characterized by the enrichment of PCV and pericytes. Administration of GnRHa significantly alleviates the immune and angiogenic microenvironment by reducing immune cells infiltration, while enhancing clinical pregnancy rates through upregulation of the granulin signaling pathway. Furthermore, our analysis highlights the therapeutic potential targeting ciliated cells in adenomyosis, although in-depth investigation into the regulatory mechanisms using diverse experimental models are still required. One limitation of our study lies in the presence of comorbidities in both the adenomyosis and control groups, which may introduce potential confounding effects. Patients with adenomyosis often present with comorbidities such as endometriosis and leiomyoma, whereas our control group included patients with leiomyoma or cervical epithelial lesions. These comorbidities are of high incidence in the general population [ 75 , 76 ]. Moreover, obtaining tissue samples from patients without any comorbidities remains particularly challenging in real-world clinical scenarios. To address these concerns, we conducted additional analyses by integrating our dataset with published data, encompassing a total of 204,627 cells categorized into control, adenomyosis without complications, and adenomyosis with complications groups. Using advanced analytical methods, including cell composition analysis, expression-based Kolmogorov-Smirnov distance, and CINEMA-OT[ 19 ], an algorithm specifically designed to assess confounding sources in scRNA-seq data, we found that the differences between adenomyosis with and without complications were minimal. This suggests that potential confounding effects from these comorbidities have only a limited impact on our overall findings. It should also be noted that our clinical samples were collected from a single center, and the sample size remains relatively small, which may to some extent constrain the robustness of our findings. Although we integrated external scRNA-seq datasets to enhance the robustness of our conclusions, future multi-center studies with larger cohorts are warranted to further validate and generalize our findings.

Adenomyosis is characterized by the presence of endometrial glands and stroma-like tissue within the myometrium, leading to abnormal uterine bleeding, dysmenorrhoea, chronic pelvic pain and infertility[ 1 ]. Around 35 % of reproductive-age women suffer from adenomyosis worldwide[ 2 ]. Symptoms variability complicates diagnosis and often leads to misdiagnosis. Advances in diagnostic imaging, such as transvaginal ultrasound scan and magnetic resonance imaging, have been accompanied by a rise in reported incidence, particularly among younger women[ 3 ]. Several hypotheses about the pathogenesis of adenomyosis have been proposed, including endometrial invagination into the myometrium, injury at the endometrial-myometrial interface, and stem cell-driven de novo lesions[ 4 ], but none of these theories explains all phenotypes. This gap in pathophysiology hinders early intervention, and highlights the need for further research. Current management of adenomyosis includes both pharmacological and surgical interventions. Most medications focus on symptom relief but do not effectively reduce lesion size. Gonadotropin-releasing hormone agonists (GnRHa) can effectively shrink lesions and are commonly used in patients undergoing assisted reproductive technology (ART). Studies have shown pre-ART GnRHa treatment has been associated with improved clinical pregnancy rates in adenomyosis[ 5 ]; however, it often induces hypoestrogenic symptoms such as hot flashes. These side effects limit treatment adherence and preclude long-term use. Advancements in single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics enable high-resolution interrogation of cellular heterogeneity and tissue architecture[ 6 ]. Recently scRNA-seq studies have identified functional differences across adenomyosis regions, offering valuable insights into cellular heterogeneity. One study suggested endometrial fibroblasts as progenitors of adenomyosis, with WNT signaling implicated in disease biology[ 7 ]. Another study highlighted the interplay between the “invagination” and “metaplasia” hypotheses, emphasizing endothelial heterogeneity and dysregulated angiogenic signaling involving vascular endothelial growth factor and angiopoietin pathways[ 8 ]. These findings underscore the complexity of adenomyosis and the importance of cellular heterogeneity. However, current studies are limited by small cohorts and insufficient spatial resolution. Moreover, it remains unclear whether GnRHa exerts distinct mechanisms in the eutopic versus ectopic endometrium, a gap in knowledge that could be critical for elucidating the molecular pathways underlying the improvement of clinical pregnancy rates with GnRHa treatment. Here, we integrate scRNA-seq and spatial transcriptomics (Geo-seq)[ 9 ] to attain a comprehensive understanding of cell types and niches in the adenomyosis uterus tissue. We also include patients who received pre-surgical GnRHa, the most common preoperative therapy, to assess region-specific treatment responses. By assigning cells to anatomical partitions, inferring ligand-receptor interactions, and localizing sources of niche factors, our data sheds light on invagination and the tissue architecture of adenomyosis.

The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.

X.L., G.L., X.Y. conceived the study. S.S. provided advice on study design. G.P. and G.C. designed the experiments. G.P., G.C. and J.X. performed experiments. C.L. and J.H. analyzed the sequencing data under the guidance from S.S. C.L. performed the data analysis, wrote the results and methods, and generated the majority of the figures. J.H. analyzed the Geo-seq data and wrote the manuscript. B.X. evaluated the impact of complications on adenomyosis. L.L. and J.H. performed immunofluorescence. Y.L. provided technical support of frozen slicing. M.L., Z.Z., J.T., C.Z. and S.C. participated in the interpretation of the results. All authors read and approved the manuscript.

All raw sequencing data are available in the National Omics Data Encyclopedia (NODE) ( https://www.biosino.org/node ) with the accession number (OEP004165) or through the URL: https://www.biosino.org/node/project/detail/OEP004165 .