Identification and Validation of MSX1 as a Key Candidate for Progestin Resistance in Endometrial Cancer.

PURPOSE: Progestin resistance is a critical obstacle for endometrial conservative therapy. Therefore, studies to acquire a more comprehensive understanding of the mechanisms are urgent. However, the pivotal molecules are still unexplored.
MATERIALS AND METHODS: We downloaded GSE121367 from the GEO database. The "limma" R language package was applied to identify differentially expressed genes (DEGs). We conducted Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA). Protein-protein interaction was constructed by STRING and visualized in Cytoscape. The tumor immune microenvironment was explored by the TISIDB database. Methylation validation and overall survival analysis were conducted by the TCGA database. In addition, the upstream modulators of hub genes were predicted by miRTarBase and Network Analyst databases. The expression levels of candidate genes were validated by quantitative real-time PCR (qRT-PCR), Western blot, and immunohistochemical assay (IHC). Cell growth, clone formation, migration, invasion, and wound healing assays were studied to explore the role of MSX1 in progestin resistance in vitro.
RESULTS: A total of 3,282 DEGs were identified and they were mostly enriched in the cell adhesion pathway. We screened out ten hub genes whose genomic alteration rates were low based on the current endometrial carcinoma sample sets. Has-miR-335-5p, has-miR-124-3p, MAZ, and TFDP1 were the most prominent upstream regulators. The methylation status of CDH1, JAG1, EPCAM, and MSX1 was decreased, corresponding to their high protein expression, which also predicted better overall survival. The homeobox protein of MSX1 showed significant tissue specificity and better prognostic value and its knockdown inhibited epithelial-mesenchymal transitions (EMT) and enhanced progesterone efficacy.
CONCLUSION: Our study identified that the gene of MSX1 promised to be the specific indicator and therapeutic target for progestin resistance. This would shed new light on the underlying biological mechanism to overcome progestin resistance of endometrial cancer.

Introduction

Endometrial carcinoma (EC), which results from aberrant regeneration in terms of excessive growth of endometrial glands,1 accounts for 4.4% of carcinoma cases among women in 20182 with more than 60,000 cases estimated in the United States in 2019.3 As for endometrial precancerous lesions including atypical hyperplasia or endometrial intraepithelial neoplasia and well-differentiated cancer, hysterectomy would not be a feasible and effective optimal choice for them and conservative treatment to preserve fertility for young patients is becoming significantly essential. While a progestin remedy is commonly applied, approximately 30% of such patients do not respond to the therapy, which causes a poor effect for fertility preservation.4 Up until now, there is no effective solution to detect or predict which group of patients may respond to the progestin treatment.

A more comprehensive exploration of the precise molecular targets of progestin resistance would facilitate further improvements in disease diagnosis and would probe new biomarkers, continuous research, including ours, has been carried out in the last decade to address the problem.5–7 At present, microarray technology and bioinformatics processing are considered to be promising tools for genomic analysis and could be well applied to identify genetic or epigenetic alterations in carcinogenesis and drug resistance, which become a necessary complement to experimental research.8 Considering recent developments in open-access datasets like the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA), the exploration of key genes and detection of functional pathways have been implemented in EC.9 The GEO and TCGA database contains thousands of clinical information and gene sequencing data, allowing for well-rounded analysis of various cancers. The gene expression profiling interactive analysis (GEPIA) acts as a web server for gene expression profiling, survival analysis, and correlative analyses on the basis of different tumor characteristics such as grades or stages.10 Epigenetics, which covers fields of aberrant DNA methylation, dysregulated noncoding RNA, and altered post-translational histone modification, refers to heritable changes in gene expression which are not associated with an alteration in DNA sequence but play an essential role in carcinogenesis and resistance detection.11 Aberrant DNA methylation is most widely explored and may become an effective detection indicator.12 Up to now, there was no relevant analysis of bioinformatics focused on progestin resistance of EC and the exploration of methylation marker of resistant genes was needed.

In this research, bioinformatics analysis was applied to reveal the differential expressed genes (DEGs) that lead to progestin resistance based on microarray datasets from GEO databases and screen out significant hub genes. Gene-related microRNAs (miRNAs), transcription factors (TFs), methylation status, and survival analysis, as well as biological functions and pathways were also integrated to explore the mechanisms and potential therapeutic value of these DEGs in resistance by constructing networks. Tumor Immune Estimation Resource (TIMER), Gene Set Enrichment Analysis (GSEA), and Gene Set Variation Analysis (GSVA) were utilized to detect underlying biological mechanisms, and in vitro experiments were utilized to validate the role of MSX1 in the presence of medroxyprogesterone (MPA). Our results may help understand the pathogenesis of progestin resistance. Moreover, it may provide insight regarding the novel treatment for EC.

Materials and Methods

Microarray Data and Data Procession

The Gene Expression Omnibus (GEO) is a public repository for data storage. In the present study, the gene expression profiling data sets (GSE121367) were obtained from the GEO database. It included endometrial cancer cell line Ishikawa, which was purchased commercially from the American Type Culture Collection (Manassas, USA), and IshikawaPR, which was established from the Ishikawa cell as its acquired a medroxyprogesterone acetate (MPA) resistant subline. Normalized data of GSE121367 was downloaded from the GEO database and further processed by the “limma” R language package to identify differentially expressed genes (DEGs) between IshikawaPR and Ishikawa cell lines. A P-value<0.05 and |log fold change (FC)|>2 were set as criteria to screen DEGs. The Cancer Genome Atlas (TCGA) database of EC was used to verify the expression status and survival function of hub genes. Subsequently, Web-based software OmicShare () and Heml () were used to draw a volcano plot and heatmap, respectively.

Functional and Pathway Enrichment Analysis

The Database for Annotation, Visualization, and Integrated Discovery (DAVID, ) is an online program offering systematic and integrative functional annotation tools for researchers to explore biological meaning behind a large list of genes.13 In this study, the DAVID database and Metascape () were introduced to perform both Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the top 250 DEGs.14 P<0.05 was set as the cut-off criterion.

Data Analysis of Gene Set Enrichment Analysis (GSEA) and Gene Set Variation Analysis (GSVA)

In order to explore biological pathways of different groups, GSEA software () was used. The annotated gene sets of c5.all.v7.0.symbols.gmt and h.all.v7.0.symbols.gmt were downloaded from the website and considered as the reference gene sets. The number of permutations was 1,000. Other parameters were set to default. A significant difference at P-value<0.05 was defined as the cutoff criteria. Normalized enrichment score (NES) and false discovery rate (FDR) were applied to determine the statistical differences. The differential results were visualized by Enrichment Map plug-in of Cytoscape.15 Furthermore, the “GSVA” R package was utilized to explore the pathways most associated with hub genes.16 On the basis of the median expression of hub gene, 91 EC samples were divided into two groups (high expression and low expression). P<0.01 was defined as statistically significant.

Protein–Protein Interaction (PPI) Network Construction and Hub Genes Screening

Firstly, an online database Search Tool for the Retrieval of Interacting Genes (STRING, ) was employed to explore the functional interactions between DEG-encoded proteins and build the PPI network.17 PPI pairs with a combined score ≥0.4 were considered as the threshold value. Subsequently, the PPI network was visualized by Cytoscape software,18 and the degree of connectivity was also analyzed. Then the network relationship file was downloaded and the top 10 hub genes were identified by the analysis tool of its plug-in (degrees ranking of cytoHubba).19

Validation of the Hub Genes in Datasets

Gene Expression Profiling Interactive Analysis (GEPIA, ) is a web-based server for Hub genes expression analysis, correlation analysis, and patient survival analysis.10 Survival analyses of hub genes were conducted by log‐rank tests and Kaplan–Meier survival curves were plotted. Then the mutation and DNA copy-number alterations of hub genes were investigated in cBioPortal (), the methylation status of hub genes was validated in Ualcan (), which were based on TCGA analysis. Furthermore, on the basis of the TCGA database, the RNA expression level of hub genes in different carcinoma tissues were detected and the difference in protein levels was displayed by immunohistochemistry (IHC) on the basis of the Human Protein Atlas database (HPA, ).20

Prediction of Relevant MicroRNAs and Transcriptional Factors (TFs) of Hub Genes

For the hub DEGs identified from the PPI network, the related miRNAs were predicted by miRTarBase (), which is the database aiming to provide hundreds of published experimentally validated miRNA–gene interactions.21 The Network Analyst () is designed to support integrative analysis of gene expression data through statistical, visual, and network-based approaches.22 In this study, Network Analyst was introduced to predict hub TFs. A list of the hub genes were enrolled into the input area and proceeded step by step, finally, the gene-related TFs as well as TFs–gene interactions pairs were presented. Then the results were visualized using the Cytoscape software, and the correlations were also evaluated based on GEPIA. All databases that have been used in this study are listed in Table 1.

Table 1

The Related Databases and Their Functions in This Study

Databases	Functions
Gene Expression Omnibus (GEO)	Extraction of gene expression microarray data
The Cancer Genome Atlas (TCGA) database	Extraction of gene expression data and survival data
The Database for Annotation, Visualization, and Integrated Discovery (DAVID)	Gene functional annotation
Metascape	Pathway enrichment analysis
Search Tool for the Retrieval of Interacting Genes (STRING)	Functional interactions between proteins
Gene Expression Profiling Interactive Analysis (GEPIA)	Hub genes expression analysis and survival analysis
cBioPortal for Cancer Genomics	Mutation and DNA copy-number alterations analysis
UALCAN	Evaluate gene expression and epigenetic regulation
Human Protein Atlas database (HPA)	Protein expression detection
TISIDB database	Explore immune microenvironment
miRTarBase database	Detect related miRNAs
Network Analyst	Predict transcriptional factors

Cell Culture and Cell Transfection

Human EC cell line Ishikawa and normal endometrial epithelial cells EEC were obtained commercially from the American Type Culture Collection (Manassas, USA) and conserved in our team. IshikawaPR as an acquired medroxyprogesterone acetate (MPA) resistant subline was established from Ishikawa cells.23–27 Briefly, Ishikawa cells were maintained in DMEM/F12 medium with 10% FBS and the synthetic progestin MPA (MedChemExpress, USA) for 10 months. The concentration of MPA was gradually increased by 2.5 µM every 4 weeks until it reached 20 µM. The medium with MPA was changed every 2 days. When the remaining surviving cells grew to a more suitable density but were less than confluent, they were passaged by 0.02% EDTA and 0.25% trypsin resolved in Hank’s balanced salt solution. After being screened for up to 10 months, cells that proliferated in 20 μM MPA were considered to be progestin-resistant Ishikawa cells (IshikawaPR). Subsequently they were kept in 20 μM MPA for further experiments. The plasmid of MSX1 was constructed in Obio Technology (Shanghai, China). A siRNA targeting human MSX1 mRNA was designed by GenePharma Biotechnology (Shanghai, China). The sequences of siMSX1 were 5ʹ-GCAUUUAGAUCUACACUCUTT-3ʹ (sense) and 5ʹ-AGAGUGUAGAUCUAAAUGCTA-3ʹ (anti-sense). The pcDNA3.1-MSX1 and siRNAs (20 μmol/L) were transfected into the sensitive and resistant subline IshikawaPR using Lipofectamine 3000 (Invitrogen, USA), respectively, when cells reached 60–70% confluency. The effect of interference was determined by Western blot and real time PCR analysis.

Western Blot Analysis

Protein extraction and Western blotting were conducted, as previously described.28 Briefly, total protein was resolved by RIPA buffer containing PMSF. Then, after electrophoresis, 50 µg protein was loaded to SDS-polyacrylamide gel and transferred to the membranes, which were incubated overnight with primary antibodies including MSX1 (bs-8512R, Bioss). After being incubated with secondary antibodies, the protein was detected using the chemiluminescence detection system. Each experiment was repeated three times.

Semi-Quantitative Reverse Transcription-PCR (RT-PCR) Analysis

Total RNA from cells was extracted using the TRIzol reagent (Thermo Fisher Scientific, USA). High-capacity cDNA reverse transcription kits (TransGen Biotech Co., China) were used to synthesize cDNA from RNA. Real-time PCR was performed using the SYBR Green PCR kit (Thermo Fisher Scientific) according to the manufacturer’s protocol. The primers used in our study are listed in Table 2.

Table 2

Sequences of Primers Used for Amplification of Target Genes

Gene	Primer Nucleotide Sequence
CDH1	Forward: 5ʹ-GAACGCATTGCCACATACAC-3’
	Reverse: 5ʹ-TGGTGTAAGCGATGGCGGCA-3
EPCAM	Forward: 5′-AATCGTCAATGCCAGTGTACTT-3′
	Reverse: 5′-TCTCATCGCAGTCAGGATCATAA-3′
MSX1	Forward: 5′-CCTCTTTGCTCCCTGAGTTCA-3′
	Reverse: 5′-GGGACTCTTCCAGCCACTTTTT-3′
KRT19	Forward: 5ʹ-GCGAGCTAGAGGTGAAGATC-3’
	Reverse: 5ʹ-CGGAAGTCATCTGCAGCCA-3’
OAS1	Forward: 5′-CGTGTTTCCGCATGCAAAT-3′
	Reverse: 5′-ACCTCGGAAGCACCTTTCCT-3′
JAG1	Forward: 5ʹ-ATCGTGCTGCCTTTCAGTTT-3′
	Reverse: 5ʹ-GATCATGCCCGAGTGAGAA-3′
TBX1	Forward: 5ʹ-CTGACCAATAACCTGCTGGATGA-3′
	Reverse: 5ʹ-GGCTGATATCTGTGCATGGAGTT-3′
DAB2	Forward: 5ʹ-GTAGAAACAAGTGCAACCAATGG-3’
	Reverse: 5ʹ-GCCTTTGAACCTTGCTAAGAGA-3’
E-cadherin	Forward: 5′-CGAGAGCTACACGTTCACGG-3′
	Reverse: 5′-GGGTGTCGAGGGAAAAATAGG-3′
N-cadherin	Forward: 5′-TGCGGTACAGTGTAACTGGG-3′
	Reverse: 5′-GAAACCGGGCTATCTGCTCG-3′
Vimentin	Forward: 5′-TGCCGTTGAAGCTGCTAACTA-3′
	Reverse: 5′-CCAGAGGGAGTGAATCCAGATTA-3′
Snail	Forward: 5′-ACTGCAACAAGGAATACCTCAG-3′
	Reverse: 5′-GCACTGGTACTTCTTGACATCTG-3′
Twist	Forward: 5′-ATTCAAAGAAACAGGGCGTGG-3′
	Reverse: 5′-CCTTTCAGTGGCTGATTGGC-3′
MMP9	Forward: 5′-TTGACAGCGACAAGAAGTGG-3′
	Reverse: 5′-GCCATTCACGTCGTCCTTAT-3′
MMP2	Forward: 5′-TCTCCTGACATTGACCTTGGC-3′
	Reverse: 5′-CAAGGTGCTGGCTGAGTAGATC-3′
CCND1	Forward:5ʹ-AAACAGATCATCCGCAAACAC-3’
	Reverse:5ʹ-GTTGGGGCTCCTCAGGTTC-3’
MYC	Forward:5ʹ-CCTGGTGCTCCATGAGGAGA-3’
	Reverse:5ʹ-TCCAGCAGAAGGTGATCCAGAC-3’
GAPDH	Forward: 5ʹ-ACCCAGAAGACTGTGGATGG-3’
	Reverse: 5ʹ-TCAGCTCAGGGATGACCTTG-3’

Immunohistochemical (IHC) Assay

Twelve endometrial cancer and paracancer tissue samples were collected from the International Peace Maternity and Child Health Hospital (IPMCH). The collection complied with the regulations of the Medical School, Shanghai Jiao Tong University, and was in accordance with the declaration of Helsinki. The sections were incubated with a MSX1 antibody overnight followed by a secondary antibody. Then, the images were captured by microscopy (Leica, Germany). The score was determined by both the intensity and the percentage of positive-stained tumor cells.28

Cell Proliferation and Colony Formation Assay

The viability of IshikawaPR cells transfected with siMSX1 or siControl and treated with or without MPA was detected by CCK-8 kit (Beyotime Biotechnology, Shanghai, China). Briefly, after the resistant cells were transfected for 0, 24, 48, and 72 hours, 5×103 cells/well in cultured medium (90 μL/well) were incubated with CCK-8 (10 μL/well) at 37°C for 2 hours. The OD value was recorded at 450 nm wavelength. Meanwhile, A total 500 transfected cells/well were plated and incubated at 37°C for 2 weeks. After cells were fixed with methanol and stained with 0.25% crystal violet (Beyotime, Shanghai, China), colonies were counted and photographed by using a microscope. All experiments were carried out in triplicate.

In vitro Migration and Invasion Assays

The 24-well chamber with polycarbonate filters (8-µm pore) (Corning, USA) was utilized and 1x105 IshikawaPR cells suspended in 100 µL serum-free medium were seeded into the upper chamber and 500 µL complete medium were added to the lower chamber. After 24 hours, the crystal violet-stained cells were counted at x200 magnification.

Cell Scratch-Wound Assay

After cells reached 90% confluence in 6-well plates, the monolayers were wounded by scratching the well with a 10 µl pipette tip. After washing twice slowly with Phosphate Buffer Saline (PBS), fresh serum-free medium was added into the well. Then, images of the scratch were captured at 0 and 48 hours by phrase-contrast microscope and calculated by using Image J software.

Statistical Analyses

Statistical analyses were performed using GraphPad Prism 7. Statistical analyses between groups were conducted with the Student’s t-test. Data are presented as mean±SD. P-values<0.05 were considered statistically significant.

Results

Identification of Aberrantly Expressed Genes

After data preprocessing and quality evaluation, we obtained the expression matrices from the cell samples in the set GSE121367 and further processed with GEO2R tool. Results showed that a total number of 3,282 common DEGs were screened out from the dataset GSE121367 with 1,819 up-regulated genes (log FC>2) and 1,463 down-regulated genes (log FC<2). The volcano plot was shown in Figure 1A. The red plots represent the up-regulated genes, green plots represent down-regulated genes, and black plots represent genes without differentialexpression based on the cut-off criteria of P-value<0.05 and |log fold change (FC)|>2. The top 250 of DEGs with significant fold change were also drawn with a heatmap (Figure 1B). The red pane shows the high expression level and green pane shows the low expression level.

Figure 1

Identification of DEGs from the GSE121367 dataset and functional enrichment analysis. (A) Volcano plot of the DEGs. Red dots and green dots represent the upregulated and downregulated genes, respectively; black represents genes with no differential expression based on the threshold of P-value <0.05 and |log FC| >2.0. (B) Heatmap of top 250 DEGs. Gene expression levels were shown by color bar. Red color denotes high level and green color denotes low level. (C) Biological process (BP). (D) Cellular component (CC). (E) Molecular function (MF). (F) KEGG pathways. (G) Boxplot of enriched terms across DEGs, colored by P-values. (H) Network of enriched terms, colored by cluster ID, where nodes that share the same cluster ID are typically close to each other.

Gene Functional Enrichment Analysis

To illustrate the biological classification of DEGs, GO term enrichment analysis was carried out using DAVID. Three categories of GO terms including biological process (BP), cellular component (CC), and molecular function (MF) results are presented in Figure 1C-E. Results suggested that changes in BP of the top 250 key genes were significantly enriched in “negative regulation of DNA binding”, “type I interferon signaling pathway”, and “neuron migration” (Figure 1C), As for the CC term, these genes showed enrichment in “anchored component of membrane”, “extracellular space”, and “multivesicular body” (Figure 1D). Besides, the MF term indicated enrichment predominantly at “protein dimerization activity”, “Notch binding”, and “transporter activity” (Figure 1E). To further analyze the DEG-enriched pathways, KEGG pathway analysis was subsequently conducted (Table 3). As shown in Figure 1F, it covered the “Cell adhesion molecules pathway” and “Endocytosis pathway”. This functional investigation identified that these DEGs had close associations with changes of DNA binding and the cell adhesion pathway. Furthermore, we also analyzed the pathway of differential genes by the website of Metascape (Figure 1G-H), which revealed that differential genes were enriched in “mesenchymal cell differential”, “cell-cell adhesion mediated by cadherin”, and “negative regulation of DNA binding and cell proliferation”.

Table 3

GO and Pathway Enrichment Analysis of DEGs

Category	Term	P-value
GOTERM_BP_DIRECT	GO:0021983~pituitary gland development	2.52E-04
GOTERM_BP_DIRECT	GO:0043392~negative regulation of DNA binding	2.52E-04
GOTERM_BP_DIRECT	GO:0060337~type I interferon signaling pathway	7.45E-04
GOTERM_BP_DIRECT	GO:0010628~positive regulation of gene expression	7.56E-04
GOTERM_BP_DIRECT	GO:0001764~neuron migration	0.001169
GOTERM_CC_DIRECT	GO:0031225~anchored component of membrane	0.000235
GOTERM_CC_DIRECT	GO:0005615~extracellular space	0.000690
GOTERM_CC_DIRECT	GO:0005771~multivesicular body	0.001670
GOTERM_CC_DIRECT	GO:0005576~extracellular region	0.002465
GOTERM_CC_DIRECT	GO:0016324~apical plasma membrane	0.004640
GOTERM_MF_DIRECT	GO:0046983~protein dimerization activity	0.006223
GOTERM_MF_DIRECT	GO:0008201~heparin binding	0.008454
GOTERM_MF_DIRECT	GO:0005112~Notch binding	0.016286
GOTERM_MF_DIRECT	GO:0046872~metal ion binding	0.016772
GOTERM_MF_DIRECT	GO:0005215~transporter activity	0.023963
KEGG_PATHWAY	hsa04514:Cell adhesion molecules	0.0046
KEGG_PATHWAY	hsa04612:Antigen processing and presentation	0.0096
KEGG_PATHWAY	hsa05168:Herpes simplex infection	0.0153
KEGG_PATHWAY	hsa04144:Endocytosis	0.0166
KEGG_PATHWAY	hsa04145:Phagosome	0.0246

Data Processing and Gene Set Enrichment Analysis (GSEA)

Although differential expression of individual genes could play a critical role in mechanistic aspects of cellular regulation, many compounds and genes are regulated complicatedly. For the sake of categorizing such modules of cellular regulation, bioinformatics approaches for “gene set enrichment” (GSEA) statistics have been developed.29 The consequences of GSEA analysis revealed that 428 gene sets were upregulated in the IshikawaPR cell line with a P-value<0.05, among which 145 gene sets were significantly enriched at a nominal P-value<0.01. A total of 116 gene sets were upregulated in the Ishikawa cell line with a P-value<0.05, among which 34 gene sets were significantly enriched at a nominal P-value<0.01. As shown in Table 4, pathways including interferon gamma response, TNF-a signaling via NF-KB, epithelial mesenchymal transition, interleukin1 beta production, and negative regulation of response to drug were significantly enriched in the IshikawaPR cell line (Figure 2B). While in the Ishikawa cell line (Table 5), the consequences showed that pathways about mesenchymal to epithelial transition, negative regulation of insulin secretion, apical junction assembly, and plasma membrane receptor complex compounds were highly enriched (Figure 2C). All the differential results of gene sets were visualized by the Enrichment Map plug-in of Cytoscape (Figure 2A). Altogether, these data suggested that when Ishikawa cells were stimulated and selected by MPA for almost 10 months, the functions of cell signal transduction such as nuclear receptor activity and cytokine biosynthetic process including interferon gamma, interleukin1 production, and epithelial cell polarity were dramatically changed.

Table 4

Upregulated Gene Sets in the IshikawaPR Cell Line

Gene Sets	Size	ES	NES	NOM P-value	FDR q-value
Interferon gamma response	196	0.55	2.07	0.00	0.00
TNF-a signaling via NF-KB	196	0.46	1.72	0.00	0.01
Epithelial mesenchymal transition	195	0.43	1.62	0.00	0.01
Hypoxia	191	0.39	1.47	0.00	0.04
Complement	193	0.44	1.67	0.00	0.01
Negative regulation of regulated secretory pathway	23	0.68	1.73	0.00	0.27
Chronic inflammatory response	18	0.69	1.72	0.01	0.27
Interleukin1 production	90	0.49	1.69	0.00	0.30
Interferon gamma mediated signaling pathway	87	0.51	1.74	0.00	0.30
Negative regulation of response to drug	25	0.63	1.68	0.00	0.30

Figure 2

The results of GSEA analysis. (A) The pathway network of the group. The red dots represent upregulated pathways. (B) Significantly enriched gene sets in IshikawaPR cell line. (C) Significantly enriched gene sets in Ishikawa cell line. MPA-R represents cell line of IshikawaPR; MPA-S represents cell line of Ishikawa; NES, normalized enrichment score.

Table 5

Upregulated Gene Sets in the Ishikawa Cell Line

Gene Sets	Size	ES	NES	NOM P-value	FDR q-value
ATP dependent microtubule motor activity plus end directed	26	−0.65	−1.76	0.00	0.62
Mesenchymal to epithelial transition	20	−0.69	−1.75	0.00	0.58
Phospholipid catabolic process	38	−0.59	−1.75	0.00	0.55
Respiratory chain complex IV	15	−0.72	−1.73	0.01	0.56
Transcytosis	18	−0.69	−1.72	0.01	0.52
Positive regulation of protein localization to cell periphery	60	−0.54	−1.69	0.00	0.59
Negative regulation of insulin secretion	36	−0.55	−1.59	0.01	0.66
Apical junction assembly	58	−0.48	−1.52	0.01	0.88
Cell–cell adhesion via plasma membrane adhesion molecules	253	−0.36	−1.39	0.01	0.94
Plasma membrane receptor complex	184	−0.36	−1.37	0.01	0.94

PPI Network Construction of DEGs and Verification of Hub Genes

The top 250 DEGs (P<0.05) of GSE121367 were used to construct a PPI network by the database of STRING and visualized in Cytoscape software (Figure 3A). The red color of a node reflects the upregulated gene and the green means the downregulated gene. The size of the node indicates the connectivity degree and the width of edge displays the combined score. PPI network analysis had been studied by using the threshold value of confidence >0.4 and connectivity degree ≥10. In this network, it contained 159 nodes and 244 edges. The plug-in of cytoHubba in Cytoscape was used to screen hub genes, then a significant submodule was obtained (Figure 3B), from which we chose the hub genes with high scores. Finally, 10 common hub genes (CDH1, JAG1, PTGES, EPCAM, CNTNAP2, TBX1, MSX1, KRT19, OAS1, and DAB2) were identified in the subnetwork (Tables 6 and 7). Next, we observed the mutation and DNA copy number alterations of 10 key genes (Figure 3C). As is shown in the OncoPrint tab, it demonstrated a visual summary of the different alterations of 10 hub genes across all sets of uterine corpus endometrial carcinoma samples based on a query of the 10 genes. Each row represents a gene, and each column represents a tumor sample. Red bars indicate gene amplifications, blue bars are deep deletions, grey bars are no alterations, and green squares are missense mutation. The genomic alteration rates of hub genes were <10% in all enrolled endometrial cancer cases. Furthermore, 10 hub genes were validated in the database TCGA to compare gene expression between endometrial carcinoma samples and normal samples (Figure 3D). The genes of CDH1, EPCAM, MSX1, KRT19, and OAS1 were overexpressed in tumor tissues, while genes including JAG1, TBX1, and DAB2 were downregulated in cancer tissues. There were no differences in the expression of PTGES and CNTNAP2 in cancer and normal samples (results not shown). In addition, we explored the expression profiles of ten hub genes in other cancers by the GEPIA database. As shown in Figure 3E, more than 10 genes were all expressed in other tissues, however, they were significantly expressed in paired endometrial tissues.

Figure 3

Identification and verification of hub genes in datasets. (A) The PPI network of top 250 DEGs. (B) The construction of submodule by the plug-in of CytoHubba in Cytoscape. (C) The OncoPrint tab showed a visual summary of the different alterations of 10 hub genes by the website of cBioPortal. (D) The protein expression of hub genes in GEPIA. *P<0.05 compared with normal endometrial tissues. (E) The expression heatmap of 10 hub genes in human cancers.

Table 6

Identification of Hub Genes by cytoHubba

Name	Betweenness	BottleNeck	Closeness	ClusteringCoefficient	Degree	DMNC	EcCentricity	EPC	MCC	MNC	Radiality	Stress
CDH1	7,440.58	76	66.299	0.085	27	0.227	0.100	43.793	96	17	11.806	17,516
JAG1	1,707.27	14	52.316	0.321	8	0.428	0.100	39.736	26	6	11.457	3,744
PTGES	510.00	1	38.399	1.000	2	0.308	0.082	21.129	2	2	10.558	500
EPCAM	423.05	6	52.267	0.378	10	0.339	0.090	41.186	58	10	11.312	1,750
CNTNAP2	2,582.66	10	38.916	0.000	5	0.000	0.100	15.668	5	1	10.609	4,672
TBX1	1,156.20	12	52.016	0.333	7	0.454	0.100	39.344	20	5	11.438	3,108
MSX1	944.50	4	43.892	0.167	4	0.308	0.100	25.883	4	2	11.039	1,976
KRT19	405.98	5	49.634	0.389	9	0.334	0.090	39.797	46	9	11.198	1,452
OAS1	81.05	1	39.275	0.500	8	0.408	0.100	29.613	58	8	10.419	450
DAB2	751.752	2	45.945	0.267	6	0.379	0.112	32.237	10	4	11.191	1,378

Table 7

The Information of Ten Hub Genes

Gene Name	logFC	P-value
CDH1(cadherin 1)	−5.41	2.16E-10
JAG1(jagged 1)	−5.44	2.35E-10
PTGES(prostaglandin E synthase)	7.77	2.42E-11
EPCAM(epithelial cell adhesion molecule)	−5.39	2.66E-10
CNTNAP2(contactin associated protein-like 2)	−7.10	1.76E-11
TBX1(T-box 1)	7.61	2.71E-10
MSX1(msh homeobox 1)	−5.45	5.79E-10
KRT19(keratin 19)	−6.60	4.00E-10
OAS1(2ʹ-5ʹ-oligoadenylate synthetase 1)	5.83	1.24E-10
DAB2(disabled homolog 2, mitogen-responsive phosphoprotein)	9.00	3.02E-10

Gene‑Associated MicroRNAs Network Analysis

To explore the potential upstream regulator of hub genes, predicted miRNAs of hub gene were analyzed by the miRTarBase database. Main miRNAs with interactions of more than two genes are listed in Table 8. Moreover, Cytoscape was used to construct the hub gene-relevant miRNAs network (Figure 4A). There were a total of eight genes, 118 miRNAs, and 128 gene-miRNA pairs contained in the network. Some miRNAs were found to play a critical role in regulating essential genes. Has-miR-335-5p was predicted to regulate PTGES, OAS1, KRT19, and DAB2, has-miR-26b-5p may regulate DAB2, CDH1, and JAG1. Furthermore, genes of PTGES and JAG1 were regulated by has-miR-124-3p, whose high expression may be associated with worse survival, suggesting that it may be involved in tumor resistance and may become a prognostic indicator for endometrial cancer.

Table 8

The Main Related MicroRNAs of Hub Genes

MicroRNAs	Genes	Count
has-miR-335-5p	PTGES, OAS1, KRT19, DAB2	4
has-miR-26b-5p	DAB2, CDH1, JAG1	3
has-miR-9500	MSX1, PTGES	2
has-miR-124-3p	PTGES, JAG1	2
has-miR-129-5p	DAB2, CDH1	2
has-miR-199a-5p	CDH1, JAG1	2
has-miR-193b-3p	CDH1, KRT19	2

Figure 4

Hub gene‑relevant MicroRNAs and Transcriptional Factors network analysis. (A) Hub gene-relevant miRNA network. Red nodes stand for hub genes, blue nodes stand for relevant miRNA, and yellow nodes stands mainly relevant miRNA. (B) Hub gene-transcription factors (TFs) regulatory network. Red nodes represent hub genes, blue nodes represent TFs, and yellow nodes represent major TFs. (C) Correlation analysis between hub genes and MAZ. (D) Correlation analysis between hub genes and TFDP1.

Core Transcriptional Factors Mediation Network Analysis of Hub Genes

To identify the transcriptional regulation of the hub genes and assess the effect of TFs on the expression of the hub genes, the gene-TFs regulation network was performed by using a Network Analyst network-based service. Totally, 143 TFs were included in the network, constructing 203 gene–TFs interaction pairs (Figure 4B). In this network, MAZ was considered as the key TF to regulate five hub genes: CDH1, EPCAM, KRT19, MSX1, and TBX1. In addition, TFDP1 plays a second important role in regulating CDH1, CNTNAP2, KRT19, and MSX1 (Table 9). Furthermore, we explored the correlation of hub genes and core TFs of MAZ and TFDP1 in endometrial carcinoma using TCGA datasets, respectively. From these results, we found that MAZ and TFDP1 had positive correlations with CDH1, EPCAM, MSX1, KRT19, and OAS1, and had negative correlations with JAG1, TBX1, and DAB2 (Figure 4C-D). Additionally, we found that MAZ was most positively related with EPCAM and negatively related to DAB2. Meanwhile, TFDP1 was positively associated with gene CDH1 and negatively associated with gene JAG1.

Table 9

The Main Related TFs of Hub Genes

TFs	Genes	Count
MAZ	CDH1, EPCAM, KRT19, MSX1, TBX1	5
TFDP1	CDH1, CNTNAP2, KRT19, MSX1	4
PPARG	CDH1, EPCAM, KRT19	3
NR2F1	CDH1, KRT19, TBX1	3
DMAP1	CDH1, KRT19, OAS1	3
EZH2	CDH1, MSX1, TBX1	3
ELF3	CDH1, EPCAM, KRT19	3
CHD1	CDH1, KRT19, PTGES	3
BCOR	EPCAM, PTGES, TBX1	3
E2F5	KRT19, PTGES, TBX1	3
JUND	EPCAM, KRT19, TBX1	3
SMARCA5	KRT19, PTGES, TBX1	3

Methylation Status and Expression Validation of Hub Genes in HPA

The initiation of cancer resistance was controlled by both genetic and epigenetic events. Epigenetic changes also make an important impact on the occurrence of drug resistance. Therefore, we decided to detect the methylation status of hub genes. As is performed in Figure 5A, the genes of CDH1, JAG1, EPCAM, and MSX1 were aberrantly methylated, which was inconsistent with their protein expression, on the basis of the Ualcan website. In addition, Immunohistochemistry (IHC) staining obtained from the HPA database showed the dysregulation of the expression of hub genes (Figure 5B), among which CDH1, JAG1, EPCAM, and MSX1 were all upregulated in endometrial carcinoma samples and the patient’s age, ID number, and staining status were listed in detail.

Figure 5

Gene methylation status and overall survival evaluation of hub genes. (A) The methylation status and protein expression of hub genes in Ualcan website. (B) The demonstration of Immunohistochemistry (IHC) staining of genes by HPA website. (C) Survival analysis of hub genes in endometrial carcinoma by GEPIA. (D) Basic expression of hub genes in different human cancer organs based on TCGA. (E) MSX1 staining of endometrial cancer and paracancer samples and the results of statistical analysis. *P<0.05 compared with normal endometrial tissues; ***P<0.001 compared with normal tissues.

Prognostic Value Evaluation and Tissue Specificity Analysis of Hub Genes

Meanwhile, we explored the prognostic values of the four essential genes further, which were obtained in GEPIA and are displayed in Figure 5C. Overall survival for endometrial cancer patients was analyzed in correspondence with the low or high expression of each gene. As is shown, high mRNA expression of CDH1 (P=0.01) was associated with better overall survival for endometrial cancer patients, along with EPCAM (P=0.045), JAG1 (P=0.02), and MSX1 (P=0.001). To further evaluate the expression level of hub genes in different human carcinoma tissues, we explored the HPA website. From Figure 5D, it indicated that the selected four hub genes had various RNA expression levels in different cancer tissues including glioma, thyroid cancer, lung cancer, colorectal cancer, head and neck cancer, stomach cancer, liver cancer, pancreatic cancer, renal cancer, urothelial cancer, prostate cancer, testis cancer, breast cancer, cervical cancer, endometrial cancer, ovarian cancer, and melanoma. Moreover, CDH1 and EPCAM displayed moderate expression levels in carcinoma of the endometrium, while JAG1 had a relatively low level. Nevertheless, MSX1 displayed the highest expression level in endometrial cancer, which showed high tissue specificity. Therefore, we examined its expression level using an anti-MSX1 antibody by IHC assay (Figure 5E), revealing that MSX1 positive staining was present at a higher level in samples with endometrial cancer than in samples of paracancer endometrium.

MSX1 is Significantly Upregulated in Human EC Resistant Cells and Its Suppression Increased Progesterone Sensitivity by in vitro Experiments

To identify whether the expression of the hub genes could be recapitulated in vitro experiments, we performed qRT-PCR assays in both the progesterone sensitive Ishikawa cell line and progesterone resistant IshikawaPR cell line. High levels of MSX1 mRNA were observed in IshikawaPR cells compared with Ishikawa cells (Figure 6A), which were the most statistically significant and may play a crucial role in the formation of progestin resistance. Therefore, the gene MSX1 was selected for the following validation. We transfected IshikawaPR cells with specific siRNA for MSX1 at a concentration of 50 nM to investigate the role that MSX1 played in progesterone resistance. As illustrated in Figure 6B, the effect of gene silencing was obvious. Knockdown of MSX1 repressed the expression of proliferation-related gene and markers associated with epithelial-mesenchymal transition. When validated in protein level, the gene of MSX1 in IshikawaPR cells was overexpressed in both Ishikawa and IshikawaPR cells, while it was weakly expressed in EEC cells. After being transfected with siMSX1, the protein levels were obviously downregulated in Ishikawa and IshikawaPR cells, respectively (Figure 6C). Then we detected the cell viability of Ishikawa and IshikawaPR cell lines when treated with the indicated doses of MPA after MXS1 overexpression and silencing. Results showed that down-regulation of MSX1 made cells more sensitive to the drug of MPA, while the up-regulation of MSX1 attenuated the inhibitory effect of the drug in both progesterone sensitive and resistant cell lines (Figure 6D). Moreover, the IC50 concentrations of MPA after transfection are listed in Table 10. Meanwhile, the knockdown of MSX1 significantly inhibited cell growth and clonogenic formation compared to an untargeted siControl and increased the sensitivity of cells to progesterone (Figure 6E-F). Besides, to evaluate whether MSX1 influenced the migratory function of IshikawaPR cells, we conducted cell migration, invasion, and wound healing assays, the downregulated MSX1 apparently decreased the migration ability of resistant cells and enhanced the effects of MPA (Figure 6G-H).

Figure 6

Validation of hub genes through in vitro experiments. (A) The RT-PCR analysis of candidate genes in Ishikawa and IshikawaPR cell lines. (B) Expression of MSX1 mRNA and other relevant molecules were evaluated by qRCR after transfection of siMSX1 for 24 hours. (C) The protein level of MSX1 in EEC, Ishikawa, and IshikawaPR cell lines, respectively. (D) Cell viability of Ishikawa and IshikawaPR cell lines when treated with the indicated doses of MPA after MXS1 overexpression and silencing. (E) Cell growth of siMSX1 transfected cells with or without MPA measured by CCK-8 assay. (F) The transfection of siMSX1 decreased clone number of IshikawaPR cells in vitro. (G) IshikawaPR-siMSX1 cells were subjected to transwell invasion and migration assays in the presence or absence of MPA. (H) Wound healing assays for IshikawaPR-siMSX1 cells and its relevant control cells. Data were shown as mean±SD; *P<0.05; **P<0.01; ***P<0.001.

Table 10

IC50 Concentrations of MPA in Sensitive and Resistant Cell Lines After MXS1 Overexpression and Silencing

Cell Lines	IC50 of MPA (µM)
Ishikawa	19.3±0.6
Ishikawa-siMSX1	15.5±0.9
Ishikawa-exMSX1	25.2±0.2
IshikawaPR	30.0±0.4
IshikawaPR-siMSX1	19.8±0.9
IshikawaPR-exMSX1	35.1±0.2

Validation of MSX1 Based on the TCGA Dataset

The relationship between MSX1 expression and different pathological grade was measured, which suggested that mRNA expression of MSX1 was significantly correlated with tumor grades (Figure 7A). To further investigate the potential functions of MSX1, we performed GSVA on the TCGA data. As shown in Figure 7B, genes in high expression groups of MSX1 were enriched in “positive regulation of intrinsic apoptotic signaling pathway by p53 class mediator” and “epithelial to mesenchymal transition” pathways. Furthermore, we enrolled survival and follow-up data from the TCGA cohort. It is suggested that high expression of MSX1 was significantly associated with favorable prognosis in EC patients (Figure 7C). We also used the tumor-immune system interactions (TISIDB) online database to detect the expression and prognostic value of MSX1 in other types of tumors (Figure 7D). It showed better prediction function in endometrial cancer (Figure 7E). In addition, we explored the association of hub genes’ expression with immune infiltration, the results showed that no or weak associations were observed between MSX1 and infiltration of lymphocytes, while we found MSX1 was positively associated with immunostimulator NT5E (Figure 7F-G).

Figure 7

Validation of hub gene in the TCGA dataset. (A) Transcriptional expression of MSX1 was significantly correlated with pathological grades of EC. (B) GSVA-derived clustering heatmap of differentially expressed pathways for MSX1. (C) The survival value of MSX1 based on TCGA data. (D) The prognostic value of MSX1 in human tumors based on TCGA cohort. (E) High expression of MSX1 is related to the better prognosis of patients with EC. (F) Spearman correlations between expression of MSX1 and immunostimulators across human cancers. (G) Correlation between expression of MSX1 and immunostimulator NT5E. *p<0.05.

Discussion

The majority of women in their reproductive period with endometrial precancerous lesion and well-differentiated endometrial neoplasm have an intense desire to preserve fertility. However, when closely following up younger patients with the progestin conservative therapy, clinicians observed that more than 30% of them responded poorly to the treatment.4 The major account for such a high failure rate is that the potential molecular mechanisms of drug resistance remain unclear. Our team has been keen on this research theme for several years.5–7,30,31 On the basis of our previous work, we recently conducted bioinformatics research replying on the dataset from the GEO website. By comparing the endometrial adenocarcinoma cell line Ishikawa with its counterpart IshikawaPR (MPA resistant cell line), we systematically analyzed their significant different genes, relevant molecular pathways, as well as their DNA copy number and the status of methylation in order to identify the essential candidate genes that promoted the carcinogenesis and progestin resistance.

In our present research, among the 24,384 DEGs, we utilized the top 250 genes to conduct functional and pathway enrichment analysis and GSEA analysis, which may provide novel insights for clarifying pathogenesis of progestin resistant. As was exhibited in DAVID, the genes were enriched in biological processes (BP) of negative regulation of DNA binding, type I interferon signaling pathway, neuron migration, and axonogenesis involved in innervation, which made complementary remarks to previous points that EGF/EGFR32 and insulin33 signaling pathways may lead to progestin resistance. Meanwhile, the KEGG pathway showed that cell adhesion molecules pathway and endocytosis signals were involved. Furthermore, the results of GSEA reported that nuclear receptor activity, chronic inflammatory response, and endothelial cell proliferation pathways were enriched in MPA resistant cells, which were significantly different from MPA sensitive cells, suggesting that prolonged progestin treatment may result in changes of cell membranes and nuclear receptors activity, affect signaling transduction, and induce peripheral inflammation response and neurological development. As was said, intestinal bacteria can induce a chronic subclinical inflammatory process, leading to insulin resistance,34 confirming the association between inflammation and resistance, while how inflammatory cytokines evoked progestin resistance remained uncovered, and the role neurological factors played in drug resistance provided a new perspective for us to explore.

The PPI network of DEGs demonstrated the functional correlations, in which hub genes were screened out. Afterwards, we investigated the mutation status of essential genes, all mutation rates were less than 10%, suggesting that their expressions were regulated by other factors, other than genomic mutations. Then, the protein expressions were testified in the TCGA database, genes with statistically significant differences were enumerated, while PTGES and CNTNAP2 had no remarkable difference between normal and tumor samples, so the results were not shown. Moreover, the relevant microRNA and transcription factors were detected and has-miR-335-5p, has-miR-124-3p, MAZ, and TFDP1 played a vital role. As was reported, has-miR-335-5p inhibits invasion and metastasis of thyroid cancer cells,35 breast cancer cells,36 and non-small cell lung cancer.37 In endometrial stromal sarcomas (ESS), has-miR-335-5p was more highly expressed in patients with tumor metastasis and relapse.38 In this study, high expression of has-miR-124-3p predicted worse survival, which was consistent with its effect on hepatocellular carcinoma.39 MAZ was thought to act as a therapeutic target for aerobic glycolysis and the progression of neuroblastoma40 and prostate cancer bone metastasis.41 TFDP1 was involved in colon cancer stemness and cell cycle progression42 and in the endometrium of women with deep infiltrating endometriosis (DIE).43 The relevant microRNA and interactions between hub genes and core TFs MAZ and TFDP1 had already been verified, suggesting they may participate in the formation of progestin resistance, while the detailed regulated mechanisms between TFs and hub genes needed to be further confirmed both in vitro and vivo.

Epigenetic processes, such as DNA methylation, are known to regulate specific gene expression.28 Therefore, in the current study, we researched the methylation status of key genes and presented the significantly different genes such as CDH1, JAG1, EPCAM, and MSX1, which corresponded to their expressions based on TCGA datasets and HPA website. According to the RNA expressions of four hub genes in different carcinoma organs depending on the TCGA database, both CDH1 and EPCAM showed moderate expression levels, while the level of JAG1 was low. Meanwhile, MSX1 was reported to demonstrate the highest expression in endometrial neoplasm, showing its high tissue specificity. Additionally, existing researches reported that MSX1 inhibited the growth and metastasis of breast cancer cells and was frequently silenced by promoter methylation,44 and that MSX1 induced G0/G1 arrest and apoptosis in cervical cancer45 and that epigenetic regulation of MSX1 associated with platinum-resistant disease in high-grade serous epithelial ovarian cancer.46 However, the expression and function of MSX1 had not been explored in endometrial resistant lesions. Therefore, although all of the four genes were associated with patients’ overall survival, due to its high tissue specificity, MSX1 may become a promising tissue specific marker to predict the therapy effect of progestin. In in vitro experiments, we found MSX1 was the most significantly upregulated gene in progesterone resistant IshikawaPR cells than other candidates and its knockdown made progesterone treatment more effective.

We further verified the expression of MSX1 in a larger TCGA cohort. The results confirmed our findings that MSX1 had a significant prediction value in EC and may make an impact on cancer progression through the P53 pathway. We also explored associations of MSX1’s expression with the tumor microenvironment and found the interaction between MSX1 and immunostimulator NT5E.47 However, a series of experimental verifications about our found mechanism need to be investigated and its immune-reinforcing effect needs further detection.

In conclusion, a comprehensive analysis of the hub genes based on the GEO dataset will likely shed new light on progestin therapy. Our study highlighted a novel understanding of the potential biological mechanism in progestin resistance and identified the homeobox gene MSX1 as a biomarker to detect the sensitivity and efficacy of progestin treatment.

Conclusions

In summary, by using comprehensive bioinformatics analysis and in vitro experiment, we have identified DEGs and demonstrated for the first time that MSX1 is likely one of the main molecular indicators of progestin resistance in endometrial cancer. If validated in a larger cohort, MSX1 may become a useful target to detect the progestin therapy effect, which is beneficial for younger patients who want to preserve fertility.