Transcriptome analysis of eutopic endometrial stromal cells in women with adenomyosis by RNA-sequencing

ABSTRACT This study aimed to identify differentially expressed genes (DEGs) and molecular pathways in eutopic endometrial stromal cells (EuESCs) from adenomyosis (AM) patients and to provide a new insight into the disease mechanisms. The gene expression profiles in adenomyotic EuESCs (A-EuESCs) and normal ESCs (N-ESCs) were analyzed by RNA-sequencing (RNA-Seq) and validated by quantitative reverse transcription-polymerase chain reaction (qRT-PCR). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses were performed to obtain insights into the functions of DEGs. The protein-protein interaction (PPI) network was constructed using the STRING database and visualized by Cytoscape software, and their hub genes were identified. A total of 458 up-/363 down-regulated genes were identified in A-EuESCs versus N-ESCs. The GO enrichment analysis showed that these genes were significantly enriched in calcium-dependent cell-cell adhesion. The most significant term of the KEGG pathway analysis was cytokine-cytokine receptor interaction. There were 145 nodes in the PPI network of the 157 DEGs, which were identified in significant enrichment pathway by the KEGG pathway analysis in N-ESCs and A-EuESCs. The PPI network revealed that IL-6 was a central hub gene. Besides, IL-6 was found as a central hub gene in the pro-inflammatory/chemotactic subnetwork, and EGF was noted as a central hub gene in the angiogenesis subnetwork. Our study indicated the alterations of transcriptomic profiles in A-EuESCs and provided new insights into the pathogenesis of AM. The A-EuESCs in women with AM have fundamental abnormalities that may predispose to pro-invasion/migration and angiogenesis.


Introduction
Adenomyosis (AM) is a common gynecological condition causing uterine enlargement, pelvic pain, menorrhagia and/or dysmenorrhea, in which and the condition is mainly refractory to drug treatments, and hysterectomy is sometimes essential for the complete alleviation of clinical symptoms [1,2]. Although several studies have concentrated on AM, but the etiology and pathogenesis of the AM have still remained elusive, which are worthy of further assessment [3,4].
As there is often a visualized histologic continuity between the ectopic and the basal endometrium in AM [5], it is widely accepted that AM results from the invagination of basalis endometrium into the myometrium through an altered or interrupted junctional zone, representing a highly specialized hormone-responsive structure located in the inner third of the myometrium [6][7][8].
Recent studies have suggested that a heritable or acquired alteration in the eutopic endometrium may play an essential role in the occurrence of AM [9][10][11][12]. Eutopic endometrium of the women with AM showed abnormal biological processes, including decreased apoptosis [8,13], increased proliferation [10] and angiogenesis [14], and impaired cytokine expression and local production of estrogens, which involved the pathogenesis of AM by enhancing the infiltration of the endometrium to the junctional zone and the growth of ectopic tissue [15]. Numerous AM-based studies have concentrated on the eutopic endometrium. Initially, gene expression profiles of endometrium from women with AM and age-matched healthy controls (HCs) had been explored with microarray platforms [16]. Subsequently, one study used RNA sequencing (RNA-Seq) to analyze differentially expressed genes (DEGs) that were in the eutopic endometrium of women with AM and HCs [17]. Recently, single-cell RNA sequencing (scRNA-seq) has been applied to identify the changes in gene expression patterns among ectopic lesions, eutopic endometrium, and normal endometrium at the single-cell level and to explore a potential novel pathogenesis of AM [18]. However, these studies were carried out without separating the endometrial stromal cells (ESCs) from glandular epithelial cells in eutopic endometrial tissue.
Invasion of abnormal adenomyotic eutopic endometrial stromal cells (A-EuESCs) has been reported in the etiology of AM [19,20]. The stromal cells may play a primary pathogenetic role in accelerating epithelial downgrowth [6]. In addition, both exogenous and endogenous interleukin-22  have enhanced the invasiveness of A-EuESCs in vitro [21]. Besides, A-EuESCs were proliferated more rapidly than normal endometrial stromal cells (N-ESCs), whether they were treated with or without estradiol (E2), medroxyprogesterone acetate (MPA), interleukin 6 (IL-6), lipopolysaccharide (LPS) and interferon γ (IFN-γ) [10,14]. A-EuECSs also expressed a higher level of AMP-activated protein kinase (AMPK) than N-ESCs [22]. As described above, studies on A-EuESCs were limited to the evaluation of expression levels of one or several particular genes [13,23,24], while few systematic studies have concentrated on the DEGs or the major pathways involved in AM from the perspective of A-EuESCs [25]. Hence, identifying distinct gene expression profiles in A-EuESCs and N-ESCs is of great significance to better understand the pathogenesis of AM.
As mentioned earlier, A-EuESCs exhibited dysregulation of pathways that globally predispose toward the development, invasion/migration, and survival of ectopic endometrial implants beyond the myometrial interface. It is reasonable to postulate the existence of intrinsic abnormalities in A-EuESCs from women with AM. In the present study, we performed the transcriptome analysis of A-EuESCs from women with or without AM by RNA-Seq, and aimed to identify DEGs and molecular pathways/networks in A-EuESCs and to provide new insights into underlying mechanisms of AM.

Patients and tissue sample collection
Female patients with AM and women without AM as age-matched HCs were enrolled in the present study. All patients had not received hormone therapy or used intrauterine contraceptive device for at least 6 months prior to surgery. Eutopic endometrial tissues were collected by hysterectomy from symptomatic women with pathologically confirmed diffuse AM, and also from female age-matched HCs by hysterectomy who had no endometrial disorders (e.g., intramural myomas) and pathologically confirmed to be free of AM or endometriosis.
All eutopic endometrial tissues were collected on the days of 20-23 of the menstrual cycle (middle-secretory phase) for isolation and culture of primary ESCs [26]. The study was approved by the Ethic Committee of Integrated Traditional Chinese and Western Medicine Hospital of Jiangsu Province (Nanjing, China; Approval No. 2019LWKYS-001). Informed consent was obtained from all participants prior to enrollment.

Primary cell culture
The isolation and culture of A-EuESCs and N-ESCs were carried out based on previously described procedures with slight modification [27]. Briefly, endometrial samples obtained by surgery were immediately placed in ice-cold sterile phosphate-buffered saline (PBS) and transferred to the laboratory. Tissues were thrice washed with sterile PBS and minced into small pieces, and then incubated with 0.1% (w/v) collagenase type II (Sigma-Aldrich, St. Louis, MO, USA) in a shaking water bath for 0.5 h at 37°C. The cell suspension was sequentially filtered through a 100-μm filter, and then, through a 40μm cell strainer (BD Falcon, Bedford, MA, USA), followed by removing the debris and epithelial cells, respectively. The cell suspension was collected and centrifuged at 200 × g for 5 min to obtain ESCs. The pellet was re-suspended in a Dulbecco's modified Eagle's medium (DMEM)/F12 (1:1) (Invitrogen, Carlsbad, CA, USA) containing 10% fetal bovine serum (FBS) (Gibico, New York, NY, USA) and 1% penicillin/ streptomycin (Invitrogen), 10 nmol/L 17estradiol (Sigma-Aldrich), and 1 umol/L medroxyprogesterone acetate (MPA) (Sigma-Aldrich) [28]. Cells were seeded at a density of 2 × 10 5 cells per T25 flask and incubated in a 5% CO 2 atmosphere at 37°C. The cultured ESCs were identified by immunocytochemical staining for vimentin and cytokeratin 8 as previously described (data were not shown) [29].

RNA-seq and data analysis
To preserve the biological properties of the ESCs, RNA-Seq was performed using primary cells (all passage 1) when the cell confluence reached 80% at 6-7 days after culture. Total RNA was extracted from cultured cells using TRIzol® reagent (Invitrogen) according to the manufacturer's instructions. RNA-Seq transcriptome library was prepared through a TruSeqTM RNA sample preparation kit (Illumina Inc., San Diego, CA, USA) using 1 μg of total RNA [30,31]. Then, mRNA sequencing was performed on an Illumina NovaSeq 6000 platform (Illumina Inc.) by Shanghai Majorbio Bio-pharm Technology Co.,Ltd (Shanghai, China). To identify DEGs in two different groups, the expression level of each transcript was calculated according to transcripts per million (TPM) reads. RSEM software was used to quantify gene abundances. DESeq2 R package (ver. 1.22.2) was utilized for annotation and differential expression analysis. DEGs were identified with adjusted P < 0.05 and absolute Log 2 fold-change (FC) ≥ 1. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses were performed to identify significant biological processes and pathways [32,33], using Goatools and KOBAS, respectively, with P-value < 0.05.
In order to explore the relationship between key genes and GO terms/KEGG pathways at a clearer glance, chord plots of GO terms and KEGG pathways were also drawn. The data were analyzed using Majorbio cloud platform (www.majorbio. com). The data discussed in the present study were deposited in the Gene Expression Omnibus (GEO) database (Accession No. GSE157718) [34].

Quantitative reverse transcription polymerase chain reaction (qRT-PCR)
Total RNA was extracted from the ESCs (all passage 2) using TRIzol® reagent as described previously. The cDNA was synthesized using the a HiScript®III Reverse Transcriptase kit (Vazyme Biotech Co., Ltd., Nanjing, China) for qRT-PCR. The thermocycling conditions of the reverse transcription were as follows: Removing genomic DNA by gDNA wiper at 42 °C for 2 min; synthesizing first-strand cDNA at 25 °C for 5 min, 37 °C for 45 sec, and at 85 °C for 5 sec. qPCR was subsequently performed using an AceQ qPCR SYBR Green Master Mix system (Vazyme Biotech Co., Ltd.) according to the manufacturer's instructions. Experiments were performed by a QuantStudio™ 5 Real-Time PCR system (Applied Biosystems, Waltham, MA, USA). The thermocycling conditions of the qPCR were as follows: pre-denaturation at 95 °C for 5 min; 40 cycles at 95 °C for 10 sec, and at 60 °C for 30 sec; and the final dissociation stage (at 95 °C for 15 sec, at 60 °C for 60 sec and at 95 °C for 15 sec) was performed at the end of the amplification procedure. The relative mRNA expression levels were normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH), as a reference gene, and calculated using the 2 -ΔΔCT method [35,36]. The qRT-PCR reactions were performed in triplicate in 96-well optical reaction plates. The GraphPad Prism 9.0 software package (GraphPad Software Inc., San Diego, CA, USA) was used to draw figures. The primers used for RT-qPCR are listed in Supplementary Table 1.

Construction of protein-protein interaction (PPI) network
The Search Tool for the Retrieval of Interacting Genes (STRING) database (ver. 11.0; https:// string-db.org/) was used to elucidate the interactive relationships of the DEGs identified in significantly enriched pathways (adjusted P < 0.05) by the KEGG pathway analysis [37]. The interacting pairs with a confidence score greater than 0.4 were considered significant and were retained. Subsequently, Cytoscape software (ver. 3.8.1) was used to establish the PPI network [38]. The network topology property indicators, including degree centrality, betweenness centrality, and closeness centrality were analyzed using CytoNCA in Cytoscape software. A node with a higher score of network topology property indicators indicated a more important role in that node in the PPI network, which was considered as a hub node [39].

Cytokine assay
N-ESCs (n = 15) and A-EuESCs (n = 15) at passage 2 were cultured at a density of 5 × 10 5 cells per a 60-mm dish. The conditioned medium from each dish was collected after 48 h of inoculation, then, centrifuged at 2,500 rpm for 5 min at room temperature, and the culture supernatants were filtered through a 0.22-μm pore-sized filter and stored at −80°C. The levels of interleukin-6 (IL-6) and epidermal growth factor (EGF) in the supernatants were determined using enzyme-linked immunosorbent assay (ELISA) kits (Elabscience Biotechnology Co., Ltd., Wuhan, China). The measurements were performed according to the manufacturer's instructions [40]. Each experiment was carried out in triplicate and repeated three times.

Statistical analysis
The statistical analysis was performed using SPSS 22.0 software (IBM Corp., Armonk, NY, USA). Continuous variables were presented as mean ± standard deviation (SD). Categorical data were expressed as number (percentage). The Chisquare and test or the Fisher's exact test, whatever appropriate, were used for data analysis. The independent two-sample t-test or the Mann-Whitney -U test, whatever appropriate, were utilized to compare the continuous variables between the two groups. P < 0.05 was considered statistically significant.

Clinical characteristics in the AM and control groups
A total of 18 female patients with AM and 18 agematched female HCs were enrolled in our study. All patients' endometrial tissues were collected by hysterectomy during the secretory phase. Besides, 18 patients with AM were multipara who aged 32 to 53 (median age, 41.6) years old; 18 age-matched female HCs were multipara who aged 34 to 55 (median age, 42.3) years old. None of the agematched female HCs suffered from AM. The clinical data of participants are presented in Table 1.

RNA-Seq analysis and identification of DEGs
To better understand the molecular mechanism of AM, we conducted a comparative transcriptomic analysis on ESCs from 3 AM and 3 age-matched female HCs. Totally, 29655 annotated mRNAs were identified, of which 821 mRNAs (458 upregulated and 363 downregulated) were significantly deregulated (Figure 1(a,b)). The top 10 genes that were significantly upregulated or downregulated are listed in Table 2. The principle component analysis (PCA) showed that the A-EuESCs exhibited distinct gene expression profiles compared with N-ESCs (Figure 1(c)). Hierarchical cluster analysis of the DEGs in ESCs from AM and control groups indicated that the gene expression patterns were clustered separately after unsupervised clustering (Figure 1(d)). The above-mentioned results suggested that the gene expression profiles in A-EuESCs in AM group were significantly altered compared with N-ESCs in control group.

Functional enrichment analysis of DEGs
To define the biological functions of the 821 DEGs, GO and KEGG pathway enrichment analyses were performed. The GO term enrichment analysis showed that the up-regulated DEGs were significantly enriched in regulation of angiogenesis (GO:0045765), tissue morphogenesis (GO:0048729) and regulation of vasculature development (GO:1901342) (Figure 4(a)). Meanwhile, the down-regulated DEGs were significantly enriched in regulation of chemotaxis (GO:0050920), cell-cell signaling (GO:0007267), and synaptic signaling (GO:0099536) (Figure 4(b)). The functions of all DEGs were significantly enriched in calcium-dependent cell-cell adhesion (GO:0016339), proteoglycan metabolic process (GO:0006029), and negative regulation of chemotaxis (GO:0050922) (Figure 4(c)).

PPI networks of the DEGs in N-ESCs and A-EuESCs
PPI networks of the 157 DEGs identified in significant enrichment pathway (adjusted P < 0.05) by the KEGG pathway analysis in N-ESCs and A-EuESCs were constructed using the STRING database and visualized by the Cytoscape software. In the PPI networks of 157 DEGs, there were 145 nodes and 793 edges (Figure 6(a)). According to the ranking of network topology property indicator of degree centrality, the top 10 nodes were separately identified as hub genes, and IL-6 was a central hub gene in the network, with the maximum number of degree (n = 61) (see Table 3). Besides, PPI subnetwork of the DEGs related to inflammatory cytokines and chemokines were  constructed as described above. In these PPI subnetworks, there were 71 nodes and 177 edges ( Figure 6(b)). As mentioned earlier, the top 10 nodes were separately identified as hub genes (Table 4), and IL-6 was also a central hub gene in the subnetworks (n = 61) (Table 4). Moreover, another subnetwork of the DEGs related to the regulation of angiogenesis was constructed, including 82 nodes and 240 edges  ( Figure 6(c)). Similarly, the top 10 nodes were also identified as hub genes (Table 5), and EGF was a central hub gene in this subnetwork (n = 55).

Quantitative analysis of ESCs-secreted protein levels of IL-6 and EGF
To further detect ESCs-secreted protein levels of IL-6 and EGF in the cell culture supernatants, two ELISA kits were used. The protein levels of IL-6 and EGF in supernatants were significantly higher in A-EuESCs than those in N-ESCs (Figure 7(a,b)).
These results further confirmed that IL-6 and EGF were key hub genes in the subnetworks of the DEGs related to the regulation of inflammatory cytokines/chemokines and angiogenesis.

Discussion
To date, no study has directly investigated the transcriptomic profile of EuESCs from women with clinically significant AM. Given the cyclical  changes in the uterine endometrium [41,42], the A-EuESCs and N-ESCs were collected in the midsecretory phase of the menstrual cycle, as confirmed by histopathological examination, for RNA-Seq. Hence, our study presented the first genome-wide view of the gene expression profiles in secretory EuESCs from female patients with AM.
In the present study, we compared the transcriptomic profiles between A-EuESCs and N-ESCs. The potential functions of dysregulated genes in A-EuESCs and N-ESCs, were predicted through the GO and KEGG pathway enrichment analyses. The GO enrichment analysis revealed that the DEGs in N-EuESCs and A-EuESCs were mainly enriched in cell-cell adhesion, angiogenesis, extracellular matrix organization, chemotaxis, etc. Meanwhile, the KEGG pathway enrichment analysis revealed that the abovementioned DEGs were mainly enriched in cytokine-cytokine receptor interaction, MAPK signaling pathway, PI3K-AKT signaling pathway, and chemokine signaling pathway.
The three tandem-arrayed protocadherin (PCDH) gene clusters, including Pcdh-α, Pcdh-β, and Pcdh-γ, play fundamental roles in the development of the vertebrate central nervous system [43]. A growing body of evidence suggested that PCDHs are widely involved in the pathogenesis and progression of multiple types of cancers by enhancing invasion and metastasis [44,45]. According to the results of chord plot analysis, several PCDH-β family member genes (PCDHB3, PCDHB4, PCDHB5, PCDHB9, PCDHB13 and PCDHB16), with a particularly pronounced downregulation in A-EuESCs were directly related to calcium-dependent cell-cell adhesion. Downregulation of PCDH-β genes may affect invasion and migration of A-EuESCs by regulating intercellular adhesion and cell spreading. Meanwhile, 2 known genes, including FOXC2 and SERPINE1, with a particularly pronounced    up-regulation in A-EuESCs were directly related to in extracellular matrix organization. Up-regulation of these genes may promote cell invasion and migration by degrading extracellular matrix. In addition, we found that the chemotactic genes, including IL-6, LPAR1, NTF3, CCL26, PTK2B, ADGRA2, ANGPTL4, and SERPINE1 were expressed through up-regulation in A-EuESCs, suggesting that A-EuESCs may have more movement tendency than N-EuESCs.
Previous studies have shown that inflammation was accumulated in the eutopic endometrium compared with in the control endometrium [46,47], which was supported by the results of the comparison between the A-EuESCs and N-ESCs in the present study. The PPI subnetworks showed that IL-6 acted as the most significant hub gene and interacted with several important inflammatory cytokines and chemokines. Recently, Xiang et al. has reported the increased mRNA expression of SERPINE1 in both eutopic and ectopic endometrium compared with that in controls during proliferative and secretory phase, while the altered expression of SERPINE1 in cellular components of endometrial tissues remained elusive [17].
Remarkably, our study further demonstrated that SERPINE1 was significantly upregulated in A-EuESCs compared with N-ESCs.
It has also been shown that angiogenesis participates in the pathophysiology of abnormal uterine bleeding and subfertility in AM [48]. Wang et al. demonstrated that A-EuESCs treated with βestradiol presented stronger pro-angiogenetic capacities, accompanied by the increased expression levels of VEGFB and ANGPTL4 proteins [14]. Furthermore, our chord plot analysis revealed that upregulated genes of FOXC2, TLR3, GATA6, HGF, IL-6, PTK2B, ANGPTL4, and HIF1A were involved in the regulation of angiogenesis ( Figure 5). Importantly, the PPI subnetworks showed that EGF acted as most significant hub gene and interacted with other angiogenic factors (HGF, ERBB2, FGF2 and SERPINE1) ( Figure 6). These results demonstrated that A-EuESCs have the characteristics of pro-angiogenic activity compared with N-ESCs. Besides, the recent scRNA-seq analysis has shown that the cell motility, cell proliferation, angiogenesis, and inflammation terms were enriched in eutopic endometrium versus normal endometrium. The DEGs were mainly functioned in angiogenesis and cell mobility-related cytoskeleton regulation and chemotaxis. However, there is no information related to stromal cell subpopulation [18]. Hence, the overlapping of our findings with the previous studies may reflect the universal mechanisms underlying the AM and the reliability of our data.
In summary, these findings, for the first time, revealed the overall characteristics of A-EuESCs from the perspective of transcriptomic profiles. Besides, there are several limitations in our study. First, endometrial samples analyzed herein were from the mid-secretory phase of menstrual cycle, and study of gene dysregulation during proliferative phase may increase the understanding of abnormalities of A-EuESCs in the pathogenesis of AM. Second, the sample size was small, and further study should be conducted with a larger sample size. Third, the dysregulated genes were only validated through qRT-PCR and ELISA in ESCs, and no validation was carried out by in situ biospecimens. Last but not least, the biological functions of several DEGs, including HOXC8, TCF21, CACNG7, FOXC2, SERPINE1, ADGRA2, NTF3, TLR3, GATA6, KCTD8 and DTX1 in the pathophysiology of AM remained unclear and should be further explored by in vitro and in vivo studies.

Conclusions
The present study provided an important basis to prior focal studies that revealed fundamental abnormalities in A-EuESCs in women with AM. The findings of the DEGs lay a foundation for further investigation to elucidate the mechanisms of AM, and the altered pathways in A-EuESCs may predispose to pro-invasion /migration and angiogenesis, which may be involved in the development of AM. The DEGs and pathways may assist scholars in the development of further efficacious therapies for AM. However, specific roles and mechanisms of the DEGs in A-EuESCs should be investigated and confirmed in the future study.

Disclosure statement
No potential conflict of interest was reported by the author(s).