Single-Cell Profiling Uncovers the Roles of Endometrial Fibrosis and Microenvironmental Changes in Adenomyosis

Adenomyosis (AM) is defined as the infiltration of ectopic endometrial glands and stroma into the myometrium, causing myometrial inflammation and fibrosis that result in dysmenorrhea, pelvic pain, uterine bleeding and subfertility. 1 Although AM is a well-known hormone-sensitive disease, progestogenic agents are not always effective. 2 AM is a benign disease, but it also has some pathophysiological features similar to malignant tumors including abnormal proliferation, angiogenesis, apoptosis, migration and invasion. 3 However, little is known about the mechanism associated with the pathogenesis of AM, and clinical treatment has always been a challenge. 4 The AM myometrium exhibits hyperperistalsis and dysperistalsis that induce chronic tissue microtrauma at the endometrial-myometrial junctional zone (EMJZ). Repeated tissue injury and repair (TIAR) leads to increased fibrosis. Accumulating evidence suggests that fibrosis plays a role in AM development. 5 Fibrogenesis can be broadly divided into four stages: initiation of tissue injury, inflammation and activation of fibroblasts, extracellular matrix (ECM) synthesis and deposition of ECM. 6 Although progressive accumulation of ECM components is a key feature leading to AM fibrosis, the pathological mechanisms underlying fibrogenesis remain unclear. Transforming growth factor β1 (TGF-β1) is induced upon tissue microtrauma, which in turn regulates the proliferation of fibroblasts at the EMJZ to help tissue repair. 7 Notably, TIAR is characterized by an interplay of macrophages, platelets and their secreted cytokines, which induces chronic inflammation at the EMJZ and promotes endometrial attachment and infiltration. 8 Furthermore, macrophages play a key role in fibrotic development. Classical M1 macrophages with pro-inflammatory function can transition to the M2 phenotype with tissue remodeling that produces ECM products. 9 Activated macrophages produce TGF-β1 and then regulate the TGF-β1/Smad3 signaling pathway in the endometrium, resulting in collagen production. 10 In addition to macrophages, mast cells are critical sentinel cells that can promote pathological conditions by releasing profibrotic mediators including TGF-β1, fibroblast growth factor (FGF). 11 , 12 More importantly, TGF-β1 signaling plays a critical role in the pathogenesis of AM, suggesting the crucial role of fibrosis in AM development. In recent years, single-cell RNA sequencing (scRNA-seq) has become an incredibly powerful tool and has been extensively used to discover novel cell types and provide insight into the molecular mechanisms of disease pathogenesis. 13 , 14 In this report, we applied scRNA-seq to comprehensively compare the transcriptional profiles of nine cell types between eutopic endometrium (AM_EM) and ectopic endometrium (AM_EC) from one AM patient.

After obtaining endometrial samples from total hysterectomies of AM patient used for sequencing, we performed HE staining to examine the histological pathology, and the gland was found to significantly invade the muscle layer in the AM_EC group with unclear structure ( Figure S1 ), indicating the AM characteristics. To explore the gene expression of the uterine endometrium, we performed scRNA-seq using cells isolated from the AM_EM and AM_EC groups ( Figure 1A ). The raw scRNA-seq data were filtered with strict quality control to reduce the bias caused by low-quality sequencing ( Figure S2A and B ). A total of 28,284 single cells were obtained, of which 8331 cells from the AM_EM group and 12,816 cells from the AM_EC group passed stringent quality control filters for further analysis ( Figure S2C , Table 3 ). All data reached a median depth of 26,247–40,066 reads per cell and 1478–1530 genes per cell. After quality control, a graph-based clustering approach and t -SNE with Seurat were used to classify principal cell clusters. A total of 19 cell clusters were detected ( Figure 1B ). Interestingly, cluster 0 was mainly present in AM_EC but not AM_EM ( Figure 1C ). Differential expression analysis was also performed between AM_EC and AM_EM ( Figure 1D and Figure 1 ); cluster-specific marker genes were also generated to determine the cellular identity of each cluster. Table 3 Summary of Single-Cell Sequencing by 10 × Genomics AM_EM AM_EC Sequencing Number of Reads 525,188,993 398,326,112 Valid Barcodes 97.1% 97.2% Reads Mapped Confidently to Transcriptome 62.0% 53.2% Reads Mapped Confidently to Exonic Regions 66.2% 58.1% Reads Mapped Confidently to Intronic Regions 19.0% 26.2% Sequencing Saturation 52.6% 40.2% Q30 Bases in Barcodes 94.6% 94.5% Q30 Bases in UMI 94.3% 94.3% Cells Estimated Number of Cells 13,108 15,176 Fraction Reads Per Cell 70.0% 68.2% Mean Reads Per Cell 40,066 26,247 Median Genes Per Cell 1530 1478 Total Genes Detected 35,309 35,808 Figure 1 Integrated single-cell RNA sequencing (scRNA-seq) analysis of patients with adenomyosis identifies diverse cell types. ( A ) Workflow diagram showing the collection and processing of obtained endometrium samples for scRNA-seq. ( B ) t -distributed stochastic neighbor embedding ( t -SNE) plot visualizing all cells displayed with different colors for samples and clusters. ( C ) Proportion of all cell clusters in AM_EM and AM_EC. ( D ) Expression of special genes in all clusters is shown in the heatmap. Each column indicates an individual cluster; yellow represents the maximum, and pink represents the minimum expression values. ( E ) Number of differentially expressed genes in each cluster. Summary of Single-Cell Sequencing by 10 × Genomics Integrated single-cell RNA sequencing (scRNA-seq) analysis of patients with adenomyosis identifies diverse cell types. ( A ) Workflow diagram showing the collection and processing of obtained endometrium samples for scRNA-seq. ( B ) t -distributed stochastic neighbor embedding ( t -SNE) plot visualizing all cells displayed with different colors for samples and clusters. ( C ) Proportion of all cell clusters in AM_EM and AM_EC. ( D ) Expression of special genes in all clusters is shown in the heatmap. Each column indicates an individual cluster; yellow represents the maximum, and pink represents the minimum expression values. ( E ) Number of differentially expressed genes in each cluster. In total, 19 clusters were classified into 9 cell types in EC and EM samples based on markers (representative marker reference from http://biocc.hrbmu.edu.cn/CellMarker/ and https://panglaodb.se/search.html ): cluster 0, endothelial cells, epithelial cells, fibroblasts, myoepithelial cells, smooth muscle cells, lymphocytes, mast cells and macrophages ( Figure 2A –, Figure S3 ). The three most abundant cell types in the AM_EM group were endothelial cells (28%), myoepithelial cells (21%) and smooth muscle cells (19%); however, endothelial cells (28%), unknown cells (21%) and epithelial cells (20%) were most abundant in the AM_EC sample ( Figure 2D ). Considering the specificity of cluster 0, it was named unknown cells. The marker gene expression for epithelial cells ( EPCAM ) and endothelial cells ( PECAM1 ) was evaluated, and we found that the proportion of merged cells was higher in the AM_EC sample ( Figure 2E and Figure 2 ). Therefore, the unknown cell type was the same as the cluster that was been deeply analyzed in our previous study. 19 Endothelial cells, epithelial cells and unknown cells will not be emphatically discussed in this report. Figure 2 Identification of endometrial cell types with scRNA-seq. Representative cell markers were used to label cell types as presented by t -SNE ( A ) and violin plots ( B ). ( C ) t -SNE plots of cells from AM_EM and AM_EC, each color indicates the associated cell types. ( D ) The pie chart shows the proportion of each cell population in the AM_EM and AM_EC samples. ( E ) EPCAM-, PECAM1- and double positive cells are shown in t -SNE plots, and the proportion of EPCAM- and PECAM1-positive cells is shown in F. Identification of endometrial cell types with scRNA-seq. Representative cell markers were used to label cell types as presented by t -SNE ( A ) and violin plots ( B ). ( C ) t -SNE plots of cells from AM_EM and AM_EC, each color indicates the associated cell types. ( D ) The pie chart shows the proportion of each cell population in the AM_EM and AM_EC samples. ( E ) EPCAM-, PECAM1- and double positive cells are shown in t -SNE plots, and the proportion of EPCAM- and PECAM1-positive cells is shown in F. The complete tables of DEGs in other eight cell types were shown in Table S1 . Some of the significant DEGs between the AM_EC and AM_EM groups are shown in Figure S4 . Notably, the number of upregulated genes was obviously higher than that of downregulated genes in all cell types ( Table 4 ). To further evaluate the biological function of these genes, we performed functional enrichment analysis, which indicated that upregulated DEGs were mainly associated with the terms angiogenesis, focal adhesion and cadherin binding. Downregulated genes were enriched in protein folding-associated terms ( Figure S5 ). Table 4 The Number of Differentially Expressed Genes in Each Cell Type Up-Regulated Genes Down-Regulated Genes Endothelial cells 424 217 Epithelial cells 522 253 Fibroblasts 366 120 Smooth muscle cells 264 98 Myoepithelial cells 324 121 Lymphocytes 369 110 Macrophages 167 105 Mast cells 252 46 The Number of Differentially Expressed Genes in Each Cell Type A previous report showed the crucial role of fibroblast-to-myofibroblast transdifferentiation (FMT) in fibrosis and AM development and the fibroblast cell type was considered the most important target cell. Thus, fibroblasts within AM_EM and AM_EC were compared, and various DEGs were screened out ( Figure 3A ). We found that COL4A1, C7, SOD2, COL14A1, THBS4 and CXCL2 were overexpressed in AM_EC sample, but CYRAB was obviously downregulated. The GO functions were enriched in ECM-related terms ( Figure 3B ), while the KEGG pathways were enriched in focal adhesion, ECM–receptor interaction, PI3K-Akt signaling pathway and MAPK signaling pathway ( Figure 3C ). To identify the potential subtypes of the entire fibroblast population, we performed unsupervised clustering of all fibroblasts to examine their heterogeneity ( Figure 3D ). Functional enrichment analysis for nine subclusters showed that ECM-related terms, including extracellular space, extracellular region, focal adhesion, collagen-containing ECM and ECM organization were mostly enriched in the majority of subclusters ( Figure 3E ). All these functions have been characterized by the expression of genes that may promote cell migration and fibrogenesis. The above results may indicate a vital role of fibroblasts in AM progression. Figure 3 Distinct characteristics of the fibroblast subpopulation. ( A ) Volcano plot represents differentially expressed genes in fibroblasts. Functional enrichment analyses with GO ( B ) and KEGG pathway ( C ) of fibroblasts. ( D ) Nine subclusters were identified after fibroblasts were combined and clustered. ( E ) Functional enrichment analyses were compared in characteristic subpopulations from AM_EM and AM_EC. Distinct characteristics of the fibroblast subpopulation. ( A ) Volcano plot represents differentially expressed genes in fibroblasts. Functional enrichment analyses with GO ( B ) and KEGG pathway ( C ) of fibroblasts. ( D ) Nine subclusters were identified after fibroblasts were combined and clustered. ( E ) Functional enrichment analyses were compared in characteristic subpopulations from AM_EM and AM_EC. Given the critical role of collagen and ECM in fibrosis, genes encoding collagen-related proteins were enriched in all cell types ( Figure 4A ). Most of the genes were highly expressed in fibroblasts, especially FBLN1 and LUM . FBLN1 is an important ECM protein involved in matrix organization regulating fibrosis-related disease. 20 LUM, known as leucine-rich proteoglycans, is involved in the regulation of collagen fibril assembly and ECM structural constituents. 21 To predict the differentiation trajectory of fibroblast subclusters, we performed a pseudotime analysis ( Figure 4B ). We also identified genes differentially expressed by states that were altered in the two branches versus the prebranch. We found that the cells undergoing fate 1 expressed higher levels of genes including GSTO1, ANXA5 and VIM , but cells undergoing fate 2 highly expressed genes such as ATF3, JUN and IGFBP5 ( Figure 4C ). Gene expression dynamics showed that the relative expression level of VIM , a marker for fibroblasts, decreased as pseudotime progressed, but that of a myofibroblast marker ( ACTA2 ) increased ( Figure 4D ), indicating the differentiation of fibroblasts into myofibroblasts. Moreover, GO functional enrichment analysis showed that fibroblasts undergoing fate 2 were obviously related to the unfolded protein response (UPR), while fate 1 fibroblasts were associated with homeostatic molecules ( Figure 4E ). UPR is considered a common pathway involved in the stress response that plays a critical pathogenic role in fibrosis. 22 Additionally, KEGG enrichment showed that fibroblasts undergoing fate 2 were more related to the IL-17 signaling pathway, protein processing in the endoplasmic reticulum, estrogen signaling pathway and MAPK signaling pathway ( Figure 4F ). Figure 4 Pseudotime trajectory of fibroblast clusters in AM_EM and AM_EC. ( A ) Heatmap of extracellular matrix-associated genes in different cell types. ( B ) Monocle analyses showing the development of fibroblasts. ( C ) Heatmap of different blocks of the top 50 differentially expressed genes along the pseudotime trajectory. ( D ) Gene expression dynamics of representative genes of myofibroblast and fibroblast markers. GO ( E ) and KEGG pathway ( F ) enrichment analyses with enriched genes in cells through fate 1 and fate 2. Pseudotime trajectory of fibroblast clusters in AM_EM and AM_EC. ( A ) Heatmap of extracellular matrix-associated genes in different cell types. ( B ) Monocle analyses showing the development of fibroblasts. ( C ) Heatmap of different blocks of the top 50 differentially expressed genes along the pseudotime trajectory. ( D ) Gene expression dynamics of representative genes of myofibroblast and fibroblast markers. GO ( E ) and KEGG pathway ( F ) enrichment analyses with enriched genes in cells through fate 1 and fate 2. Emerging evidence suggests that smooth muscle metaplasia plays a role in AM and fibrosis in endometriosis. 23 , 24 To study the function of smooth muscle cells in AM, we assessed the average expression of the DEGs ( Figure 5A ). COL4A1 , which is related to ECM, and CXCL8 and IL24 , which are related to inflammation, were significantly upregulated in AM_EC, but the levels of HSPA1A, HSPA6 and CRYAB were decreased. The results of functional enrichment analysis showed that most DEGs were associated with focal adhesion, extracellular exosome and RNA binding ( Figure 5B ). Pathways associated with cancer were enriched in genes differentially expressed AM_EC ( Figure 5B ). Then, subclustering of these cells was performed to better characterize changes between AM_EC and AM_EM. This analysis resulted in the identification of nine distinct subclusters ( Figure 5C ). The expression of some genes in each subcluster is shown in Figure 5D . Moreover, functional analysis of each subcluster showed that subclusters 0, 4, and 8 were associated with focal adhesion, suggesting their roles in smooth muscle metaplasia ( Figure 5E ). These data indicate that the heterogeneity of smooth muscle cells contributes to AM progression. Figure 5 Distinct characteristics of smooth muscle cell subclusters. ( A ) Volcano plot of differentially expressed genes in smooth muscle cells. ( B ) Functional enrichment analyses with GO terms and pathways. ( C ) Subclusters of smooth muscle cells identified from AM_EM and AM_EC. ( D ) Heatmap representing the top ten positive marker genes within each cluster. ( E ) GO functional enrichment analysis of subclusters from two samples. Distinct characteristics of smooth muscle cell subclusters. ( A ) Volcano plot of differentially expressed genes in smooth muscle cells. ( B ) Functional enrichment analyses with GO terms and pathways. ( C ) Subclusters of smooth muscle cells identified from AM_EM and AM_EC. ( D ) Heatmap representing the top ten positive marker genes within each cluster. ( E ) GO functional enrichment analysis of subclusters from two samples. Given that adenomyotic lesions are wounds undergoing repeated injury and repair, the microenvironment plays critical roles in AM development. We sought to further characterize the microenvironment of AM lesions to analyze the associated immune cells using scRNA-seq data. First, we compared gene expression in lymphocytes, macrophages and mast cells between AM_EM and AM_EC samples ( Figure 6 ). CRYAB, TPM2 and TAGLN were obviously downregulated in all three cell types in the AM_EC sample compared with the AM_EM sample. Additionally, compared with those in the AM_EM sample, PTPRC, CXCL8 and MMP1 were upregulated in lymphocytes ( Figure 6A ); the expression of CREG1, THBS1 and COL3A1 was significantly increased in macrophages ( Figure 6D ); CPA3, KIT, and IL1RL1 were overexpressed in mast cells ( Figure 6F ). Figure 6 Distinct characteristics of immune cells. ( A ) Volcano plot indicates gene expression differences in lymphocytes. Bubble charts for results of GO functional ( B ) and pathway ( C ) enrichment analysis of dysregulated genes from lymphocytes between AM_EC and AM_EM. ( D ) Volcano plot indicates gene expression differences in macrophages. ( E ) Enriched GO and pathway terms for macrophages. ( F ) Volcano plot representing differentially expressed genes in mast cells. ( G ) Enriched GO and pathway terms for mast cells. Distinct characteristics of immune cells. ( A ) Volcano plot indicates gene expression differences in lymphocytes. Bubble charts for results of GO functional ( B ) and pathway ( C ) enrichment analysis of dysregulated genes from lymphocytes between AM_EC and AM_EM. ( D ) Volcano plot indicates gene expression differences in macrophages. ( E ) Enriched GO and pathway terms for macrophages. ( F ) Volcano plot representing differentially expressed genes in mast cells. ( G ) Enriched GO and pathway terms for mast cells. GO and KEGG enrichment analyses revealed that lymphocytes were enriched for genes associated with protein binding, cytoplasm, focal adhesion, estrogen signaling pathway and cellular senescence ( Figure 6B and Figure 6 ). Both macrophages and mast cells showed focal adhesion, cadherin binding, the PI3K-Akt signaling pathway and the MAPK signaling pathway ( Figure 6E and Figure 6 ). The above results suggest that these genes in lymphocytes, macrophages and mast cells may be associated with AM progression. To systematically assess the associated complex cellular responses and better understand cellular behavior as well as response to neighboring cells in AM, we mapped ligand–receptor interaction pairs to our scRNA-seq data. The expression levels of ligands and receptors were considered in each cell type and molecular interactions between cells were predicted in AM_EM ( Figure 7A ) and AM_EC samples ( Figure 7B ). The results showed that endothelial cells and epithelial cells interacted with other cell types more closely in both AM_EM and AM_EC, indicating mild pathological changes in the AM_EM sample. Furthermore, higher levels of interactions between fibroblasts and other cell types, including lymphocytes, macrophages, mast cells and smooth muscle cells, were found in the AM_EC sample than in the AM_EM sample. Functional enrichment of ligands showed that they were mainly related to locomotion and biological adhesion ( Figure 7C ). Moreover, the crosstalk among some cell types was studied by evaluating the expression patterns of ligand-receptor. Notably, fibroblasts in AM_EC expressed relatively higher levels of collagen molecules, while the corresponding receptors are widely expressed in immune cells and smooth muscle cells, which could enhance the ECM secretion and adhesion of cells ( Figure S6 ). Considering the ECM-related functions enriched in these genes, crosstalk between fibroblasts, macrophages and other cell types plays important roles in AM progression. Figure 7 Heatmaps showing the interaction events in AM_EM ( A ) and AM_EC ( B ). Color bars from blue to red indicate interaction events calculated between two cell types. ( C ) GO functional enrichment analysis of ligands in two samples. Heatmaps showing the interaction events in AM_EM ( A ) and AM_EC ( B ). Color bars from blue to red indicate interaction events calculated between two cell types. ( C ) GO functional enrichment analysis of ligands in two samples. To verify the results of scRNA-seq, we selected 10 DEGs among the top 10 up- and downregulated genes for RT-qPCR ( Figure 8 ). The levels of STC1 (P = 0.0075), MMP1 (P = 0.0291), CXCL8 (P = 0.0385) and COL4A1 (P = 0.034) were higher in AM_EC samples than in AM_EM samples, which was consistent with the scRNA-seq results. However, there was no significant difference in ESM1 (P = 0.2542), which may be due to our limited AM samples for verification. In contrast, levels of CRYAB (P = 0.0162), HSPA6 (P < 0.0001), TAGLN (P = 0.0077), HSPA1A (P = 0.0036) and TPM2 (P = 0.0354) were lower in AM_EC samples. Figure 8 Relative expression of genes was evaluated with Reverse Transcription Real-Time PCR. *P < 0.05; **P < 0.01; ***P < 0.001. Relative expression of genes was evaluated with Reverse Transcription Real-Time PCR. *P < 0.05; **P < 0.01; ***P < 0.001.

This study was approved by the Ethical Committee of The Affiliated Hospital of Shandong University of Traditional Chinese Medicine (No. 2019-061-KY) and done according to the principles of the Helsinki Declaration. One 46-year-old woman (sample used for scRNA-seq), with excessive menstruation, pain and pathological findings showing invasion of endometrium into the myometrium diagnosed with AM provided signed informed consent for the collection of uterine tissue. Fresh endometrial tissues from the AM patient mentioned above were obtained through operation under the supervision of experienced gynaecological surgery experts. AM_EM and AM_EC samples were also collected from women with AM and verified by hematoxylin–eosin (HE) staining. Three other AM patients were used for key result verification using the same methods, and information about all four patients used in the study is listed in Table 1 . All patients in the study presented regular menstrual cycles but did not receive any hormone treatment within the last three months. Table 1 Clinical Characteristics of AM Patient Samples Used in This Study ScRNA-seq Gender Age Infertility Indication for Surgery Operative Type AM (AM_EM & AM_EC) Female 46 No Excessive menstruation, pain Laparoscope Validation Gender Age Infertility Indication for Surgery Operative Type AM (AM_EM & AM_EC) Female 52 No Pain Laparoscope AM (AM_EM & AM_EC) Female 47 No Excessive menstruation Laparoscope AM (AM_EM & AM_EC) Female 50 No Pain Laparoscope Clinical Characteristics of AM Patient Samples Used in This Study Under sterile conditions, samples of the eutopic endometrium and ectopic lesions from that AM patient were collected and immediately saved in cold saline. Each sample was minced on ice to less than 1 mm 3 pieces using eye scissors, followed by collagenase IV (Gibco, ThermoFisher Scientific, USA) digestion in a 37°C water bath with shaking for 20 min at 100 rpm. Samples were then centrifuged at 300 relative centrifugal force (RCF) for 30s at room temperature. After removal of the supernatant without disturbing the cell pellet, samples were maintained in 1 × PBS (Gibco, ThermoFisher Scientific, USA, calcium, and magnesium free) containing 0.04% weight/volume BSA (Gibco, ThermoFisher Scientific, USA, 400 μg/mL) and then centrifuged at 300 RCF for 5 min at room temperature. The cell pellet was resuspended in 1 mL red blood cell lysis buffer and incubated for 10 min at 4°C. After red blood cell lysis, samples were resuspended in 1 mL PBS containing 0.04% BSA and then filtered over Scienceware Flowmi 40-μm cell strainers to remove large particles, erythrocytes and dead cells. Finally, cell concentration and cell viability were determined by hemocytometer and Trypan Blue staining. ScRNA-seq libraries were constructed using Chromium Single Cell 3’ Reagent Kits, version 3, according to the manufacturer’s protocol. Cellular suspensions were placed on a Chromium Single Cell Controller Instrument (GCG-SR-1, 10× Genomics, CA, USA) to obtain single-cell gel beads in emulsion (GEM). Briefly, single cells were suspended in 1 × PBS containing 0.04% BSA free of magnesium and calcium, and then added to each channel for the generation of GEMs. Reverse transcription reactions were performed to generate barcoded full-length cDNA, followed by the disruption of droplets and cDNA was purified with DynaBeads Myone Silane Beads (Thermo Fisher Scientific, MA, USA). The amplified cDNA was then subjected to fragmentation, end repair, A-tailing, adapter ligation and library amplification. Finally, these libraries were sequenced on an Illumina HiSeq XTEN platform (Thermo Fisher Scientific, MA, USA). All sequencing data were deposited in the BioProject of the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession numbers SRR12791872 (AM_EM) and SRR12791871 (AM_EC). The Cell Ranger software pipeline (version 4.0.0) was applied to conduct sample demultiplexing, barcode processing and mapping reads to the reference genome (human GRCh38). STAR aligner was used in FASTQ alignment. 15 The unique molecular identifier (UMI) count matrix from all samples was read into the Seurat R package (version 3.1.1, https://github.com/satijalab/seurat ). 16 Before clustering, data filtering was conducted by removing cells that expressed fewer than 500 genes, and had mitochondrial expression exceeding 25%. After filtering, we normalized the expression values for each cell by the total expression using the Seurat Normalize Data function. Then, dimensionality reduction was performed on the scaled data by computing the significant principal components (PCs), and the first PCs were obtained for downstream analysis. Principal component analysis (PCA) was run to reduce the number of dimensions of the top variable genes on the log transformed gene-barcode matrices. Based on the Euclidean distance in PCA space, a k-nearest-neighbor graph was constructed using the “FindNeighbors” function. This graph is transformed to encode Jaccard network similarities to represent phenotypic similarity. Shared nearest-neighbor graph-based clustering was also performed using the “FindClusters” function of the Seurat package. Cells were clustered and visualized in a t -distributed neighborhood embedding ( t -SNE). The top variable genes in single cells were identified, and the average expression for each gene was calculated. Specific marker genes in each cluster were identified by the “FindAllMarkers” function, and conventional markers presented in a previous study were used to cluster every cell into a known cell type. Seurat software in R was applied to construct violin plots and perform differential gene expression analysis. Genes that were differentially expressed between cell types and subtypes were screened out using Seurat software. Bimod likelihood-ratio test on all genes was used to identify 0.26-fold (log scale) enriched genes detected in at least 10% of cells in the cluster with a P value less than 0.01. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed via hypergeometric distribution using the ClusterProfiler package, and bubble charts were plotted using the ggplot2 package in R (P < 0.05). For exploration of the enriched receptor–ligand interactions between cell types, CellPhoneDB, 17 which is a publicly available repository of receptors, ligands and their interactions, was used with the custom argument “-result-precision=4”. Ten significantly differentially expressed genes (DEGs) related to fibrosis were selected for RT-qPCR in AM patients from a separate cohort of women (n=3): STC1, MMP1, ESM1, CRYAB, HSPA6, HSPA1A, COL4A1, CXCL8, TPM2 and TAGLN . The sequences for primers are listed in Table 2 . Total RNA from the endometrium was isolated using the RNeasy Mini Kit (BioFlux, China) according to the manufacturer’s protocol. Then, cDNA was synthesized using the PrimScript RT Reagent Kit (TOYOBO, Japan). The levels of GAPDH were used as an internal control for mRNA quantification. RT-qPCR was assayed in a LightCycler 480 system (Roche, Switzerland). All data output from experiments were analyzed using the 2 −ΔΔCT method. Student’s t -test was used to evaluate the differences between samples after normality testing. Table 2 Primer Sequences Gene Name Forward Primer Reverse Primer COL4A1 GGACTACCTGGAACAAAAGGG GCCAAGTATCTCACCTGGATCA CYRAB AGGTGTTGGGAGATGTGATTGA GGATGAAGTAATGGTGAGAGGGT STC1 CACGAGCTGACTTCAACAGGA GGATGTGCGTTTGATGTGGG CXCL8 TTTTGCCAAGGAGTGCTAAAGA AACCCTCTGCACCCAGTTTTC TPM2 GGTGGCCGAGAGTAAATGTGG TTTGGTGGAATACTTGTCCGC TAGLN AGTGCAGTCCAAAATCGAGAAG CTTGCTCAGAATCACGCCAT MMP1 CTCTGGAGTAATGTCACACCTCT TGTTGGTCCACCTTTCATCTTC ESM1 TGGTGAAGAGTTTGGTATCTGC TTTTCCCGTCCCCCTGTCA HSPA6 CAAGGTGCGCGTATGCTAC GCTCATTGATGATCCGCAACAC HSPA1A AGCTGGAGCAGGTGTGTAAC CAGCAATCTTGGAAAGGCCC Primer Sequences Trajectory analysis was performed for fibroblasts using the Monocle 2 package (version 2.8.0). 18 Differential gene expression analysis of fibroblasts was then conducted using the differential gene test function to identify significant genes (q value <0.01). The dimension reduction function was applied for cell ordering in an unsupervised manner, and then trajectory construction was performed with default parameters. Endometrium samples from an AM patient used for scRNA-seq were fixed with 4% paraformaldehyde (PFA, Servicebio, China) for 48h, embedded in paraffin and sliced into 4-μm sections with a microtome. After deparaffinization and rehydration, the pathological section was stained with hematoxylin and eosin. The slides were sealed with neutral gum, and stained sections were captured under a light microscope (Nikon Eclipse E100, Japan). GraphPad Prism (Version 5.0.1, GraphPad Software Inc., San Diego, USA) was applied to analyze mRNA expression levels. The data were collected and shown as mean ± standard error of the mean (SEM). Statistical analysis was performed by t -test. P < 0.05 was considered to be significant difference. Each assay was replicated three times.

AM is a benign gynecologic disorder that is characterized by infiltration of the basal endometrium into the underlying myometrium; moreover, extensive fibrosis can be secondary to the pathological symptoms mentioned above in women with AM. 25 Although increasing evidence indicates cell heterogeneity between AM_EM and AM_EC, our understanding of the cell-type complexity underlying AM development remains limited. In this study, we presented scRNA-seq survey of various cell types in the endometrium elucidating cell type-specific gene expression signatures. Notably, we identified nine cell types: unknown cell types, myoepithelial cells, epithelial cells, endothelial cells, smooth muscle cells, fibroblasts, lymphocytes, macrophages and mast cells. Importantly, we defined DEGs for each cell type and discovered that ECM-related terms and inflammation as well as focal adhesion were dysregulated in most cell types, indicating the contributions of endometrial cell types to AM. We also identified alterations in cell–cell communication in AM and found that cell interactions were increased in the AM_EC group. Here, we obtained endometrial samples from one AM patient and successfully isolated a subcluster (cluster 0) with high CNV levels and tumor-related characteristics for the first time. 19 Moreover, specific markers for epithelial and endothelial cells were colocalized in cluster 0, suggesting an important role of these two cell types in AM progression. The detailed function of cluster 0, epithelial and endothelial cells has been described in a previous study; therefore, other cell types were focused on in the present work. Notably, fibroblasts, as the main cell population, play a key role in the AM process. The contribution of fibroblasts that differentiate from fibrocytes to tissue fibrosis has been discussed for several years as well. scRNA-seq was applied to distinguish fibroblasts and identify AM-associated fibroblast subtypes. Significant genetic differences in fibroblasts between AM_EC and AM_EM may explain the occurrence and progression of AM. A previous study also demonstrated the significance of fibrosis in AM. 1 We identified multiple upregulated genes related to AM in fibroblasts, such as STC1 and COL4A1 , which have been demonstrated to be involved in tissue fibrosis. STC1 is a hypocalcemic glycoprotein that regulates calcium homeostasis, early wound healing and macrophage functions to maintain endothelial and epithelial homeostasis. 26 Moreover, STC1 is a downstream effector molecule of VEGF that facilitates the binding of VEGF and VEGFR2 to enhance neovascularization. 27 Previous reports revealed roles of STC1 in fibrosis of the kidney, 28 lung 29 and liver. 30 Additionally, the MMP family is devoted to maintaining ECM homeostasis and collagen formation, which is important for endometriosis fibrosis, suggesting a vital role of MMP1 in AM development. 31 We also found that a variety of collagen associated genes were highly expressed in fibroblasts, including FBLN1 and LUM . FBLN1 is an important ECM protein as a secreted glycoprotein that facilitates the stabilization and binding of other EMC components during collagen deposition. 32 , 33 LUM, a small leucine-rich ECM proteoglycan, has been implicated in collagen homeostasis and fibrogenesis. 34 , 35 Furthermore, our data identified an altered process in the fibroblast trajectory, suggesting that FMT may exist in this process. DEGs from fibroblasts were mainly involved in pathways of focal adhesion, ECM–receptor interaction and PI3K-Akt signaling. Cell adhesion can stabilize the interaction between ECM and endometrial cells. 36 Importantly, PI3K/Akt, activated by TGF-β1, has been reported to play a key role in fibrosis, 37 suggesting that pathogenic changes in fibroblasts may directly mediate progressive AM. In addition, fibrosis may be associated with wounds undergoing repeated tissue injury and repair, which contribute to processes of epithelial–mesenchymal transition (EMT), FMT and smooth muscle metaplasia. 38 In AM patients, EC tissue induces abnormal smooth muscle cells, including hyperplasia and hypertrophy. 25 We found increased levels of COL4A1 and CXCL8 in EC group compared with the AM_EM group. Fibrosis is defined by excessive ECM deposition, and collagen is the most abundant component of ECM. COL4A1 is deposited in early fibrotic lesions that may be implicated in the migration of lung lesion fibroblasts, suggesting its vital role in fibrosis. 39 Furthermore, most genes in smooth muscle cells were involved in antigen processing and presentation pathways, suggesting that inflammation of smooth muscle may play a role in AM progression. CXCL8 can recruit myeloid cells to the tumor microenvironment that directly or indirectly affects the function of almost all immune cell types and thus influence tumor development. 40 The results revealed the interactions between smooth muscle cells and immune cells in AM. As the most frequent cell type in fibrotic tissues, activated fibroblasts produce cytokines, growth factors and reactive oxygen species maintaining fibrogenesis and attracting immune cells to promote chronic inflammation. 41 A previous study has shown that immune microenvironment changes related to regulatory T cells, B cells, macrophages and mast cells have been linked to the pathogenesis of fibrosis. 42 , 43 M1- and M2-like macrophages have been found to promote lung fibrosis by secreting pro-inflammatory factors and contributing to the immune response. 44 Moreover, a previous study found increased M2 macrophages in endometriosis lesions that could decrease the immune response and induce angiogenesis and tissue repair. 45 , 46 Nevertheless, function of macrophages in AM development has rarely been studied. Our results showed that DEGs in macrophages and mast cells, including downregulated CRYAB, HSPA6, HSPA1A, TAGLN and TPM2 , were mainly involved in focal adhesion, the PI3K-Akt signaling pathway, the MAPK signaling pathway and VEGF signaling pathway, suggesting the important role of the immune system in modulating the transition of inflammation to fibrosis in AM patients. Importantly, CRYAB, as an LBH (limb-bud and heart)-interacting protein, has been proven to be involved in TGF-β1-induced fibrosis by regulating FMT and EMT-like processes. 47 TAGLN is considered an F-actin binding protein that regulates the organization of the actin cytoskeleton, motility and cellular contractility. Furthermore, this molecule is also found in the EMT process, induced by TGF-β1 signaling and enhanced by IL-1β signaling. 48 Importantly, TPM2 can stabilize F-actin filaments and is regulated by TGF-β1 to induce epithelial to myofibroblastic transition in fibrosis. 49 However, the levels of CRYAB, TAGLN and TPM2 were decreased in AM_EC compared with AM_EM samples, and thus, their functions will be further studied. Communication among cells in AM lesions plays a critical role in the progression of this disease. Notably, intercellular communication was increased in AM, such as interactions between fibroblasts, smooth muscle cells and macrophages. In addition, we observed high expression of collagen-associated ligands in AM_EC fibroblasts, suggesting an imbalanced ECM environment in this disease. In addition, high expression of fibroblast receptors associated with angiogenesis and adhesion, such as VEGFR, CD44 and CD74, was activated by ligands, indicating that the altered microenvironment played an important role in promoting AM development. Our study has some limitations. First, due to the limited sample size and a heterogeneous clinical presentation of AM patients, we could not statistically analyze samples or fully reflect the characteristics of AM lesions. Second, computational analyses for scRNA-seq data could not reflect the true biological effect related to disease conditions. With regard to gene dysregulation in this study, the levels of some genes were assessed using RT-qPCR. However, future functional validations including immunohistochemistry and RNA in situ hybridization would provide more data on AM progression. In conclusion, our comprehensive characterization of cells at the single-cell level from endometrium of AM patients revealed the cell composition and gene expression pattern in both AM_EM and AM_EC tissues. The single-cell transcriptome atlas presented in our study revealed cell type-specific gene alterations in fibroblasts involved in ECM and focal adhesions contributing to fibrosis in AM progression. Moreover, altered cell–cell communication may be indicative of disorders of the immune microenvironment. Therefore, we speculate that FMT-associated processes and the immune response may exert critical functions in fibrogenesis in AM. Under the combined action of different cells, the microenvironment in the endometrium changes. Thus, these findings indicate the potential role of fibrosis and microenvironment changes in AM progression and broaden our understanding of the critical molecular mechanism for AM as well as improve current therapeutic strategies for this disease.