Computational analysis for identification of the extracellular matrix molecules involved in endometrial cancer progression.

Recurrence and poorly differentiated (grade 3 and above) and atypical cell type endometrial cancer (EC) have poor prognosis outcome. The mechanisms and characteristics of recurrence and distal metastasis of EC remain unclear. The extracellular matrix (ECM) of the reproductive tract in women undergoes extensive structural remodelling changes every month. Altered ECMs surrounding cells were believed to play crucial roles in a cancer progression. To decipher the associations between ECM and EC development, we generated a PAN-ECM Data list of 1516 genes including ECM molecules (ECMs), synthetic and degradation enzymes for ECMs, ECM receptors, and soluble molecules that regulate ECM and used RNA-Seq data from The Cancer Genome Atlas (TCGA) for the studies. The alterations of PAN-ECM genes by comparing the RNA-Seq expressions profiles of EC samples which have been grouped as tumorigenesis and metastasis group based on their pathological grading were identified. Differential analyses including functional enrichment, co-expression network, and molecular network analysis were carried out to identify the specific PAN-ECM genes that may involve in the progression of EC. Eight hundred and thirty-one and 241 PAN-ECM genes were significantly involved in tumorigenesis (p-value <1.571e-15) and metastasis (p-value <2.2e-16), respectively, whereas 140 genes were in the intersection of tumorigenesis and metastasis. Interestingly, 92 of the 140 intersecting PAN-ECM genes showed contrasting fold changes between the tumorigenesis and metastasis datasets. Enrichment analysis for the contrast PAN-ECM genes indicated pathways such as GP6 signaling, ILK signaling, and interleukin (IL)-8 signaling pathways were activated in metastasis but inhibited in tumorigenesis. The significantly activated ECM and ECM associated genes in GP6 signaling, ILK signaling, and interleukin (IL)-8 signaling pathways may play crucial roles in metastasis of EC. Our study provides a better understanding of the etiology and the progression of EC.

Introduction

Extracellular matrix (ECM), a critical component of tissues, is an interlocking mesh of fibrous proteins and polysaccharides including collagen, elastin/nidogen, laminins, and proteoglycans that form a well-organized structure in multicellular organisms [1]. The basic structure of the ECM complex, as described in Fig 1A, provides a supporting structure for cells to anchor [2]. In addition, soluble molecules, such as glycopeptidases and growth factors, are present in spaces among cells and the ECM structure, may regulate cell proliferation, differentiation, and the movement [3], as well as many secreted signaling molecule such as semaphorin III family members genes (SEMA3F) also plays there role in Endometrial cancer [4]. The ECM networks and various cytokines and growth factors create a microenvironment for cell growth and development. ECM receptors at the cell surfaces mediate 2-way signals for cells to sense the microenvironment and react to the stimuli [5]. In cancer patients, tumor cells do not grow along, and their behaviours are regulated by their own biology and also by interacting with the surrounding microenvironment [6]. Many components of ECM, like fibrous collagen and elastin create tissue stiffness [7]. Increases in collagen linearization and tissue stiffness in precancerous tissues of breast cancer were reported, and such changes increase the tumor incidence and progression [7]. Proteolysis is a significant factor that leads to metabolic reprogramming of ECM dynamics, and disruption of these control mechanisms causes ECM dysregulation and disorganization, which may contribute to malignant processes or cancer [8, 9]. The ECM paves the way for cancer cell development by creating a favourable microenvironment, more-intense stiffness, etc. [10].

Fig 1

Extracellular matrix (ECM) components and their alterations during tumor development and progression.

(a) ECM components, including proteoglycans, collagens, integrins, and many other substances, form a network structure, and (b) alterations of ECM networks may allow tumor cells to migrate from the original site to distant organs and develop into tumors.

Endometrial cancer (EC) is a commonly occurred malignancy of the female reproductive tract that arises from the uterus lining [11]. While the occurrence of the disease varies widely among countries, EC has become the most common female cancer in areas like North America, Europe, and middle-income developing countries such as South Africa and India [12]. In most cases, women with more-aggressive EC have higher-grade tumors, and the disease spreads from the uterus within 1 year [13]. About 75% of EC patients are diagnosed at an early stage (International Federation of Gynecology and Obstetrics (FIGO) stages I and II) and have an excellent prognosis with a 5-year survival rate of 80% [14, 15]. However, survival rates decrease in severe cases, at 52% and 27% for patients at stages III and IV, respectively [15]. EC is classified into two groups, type I and II endometrioid tumors [16]. Type I is estrogen-dependent [17], obesity is the major risk factor [18], and it has a favourable prognosis; in contrast, type II tumors occur in elderly, non-obese women, are estrogen-independent and exhibit worse outcomes [16, 19]. The most common treatment for patients with EC is surgery combined with adjuvant radiotherapy for early-stage type I EC, and combined with chemotherapy for type II EC [11]. Many studies showed that women with a family history of EC have a higher risk of developing EC in their lifetime [20]. The pathological classification of EC is based on the FIGO system [21, 22]. Stage I cancer is described as the tumor only being found in the uterus or womb and not having spread to other parts of the body. Stage II or beyond is described when a tumor has spread out of the uterus, depending on the migration site and invasion of lymph nodes.

ECM plays a critical role in the female endometrium, being involved in major remodelling changes every month. The ECM of the reproductive tract in women undergoes extensive structural remodelling for decidualization, implantation, and renewal of the endometrial layer each month, while abnormal ECM changes can contribute to processes such as endometriosis, infertility, and cancer development and metastasis [23]. The invasive property is considered as cancer cells intrude into surrounding normal tissues and metastasize. Alterations of ECM, ECM-associated molecules, and the local ECM microenvironment that can promote tumor progression are illustrated in Fig 1. Some of the ECM molecules such as aggrecan, nidogen, collagen type VIII chain α1 and type XI chain α2, were overexpressed in stage III EC [24], matrix metalloproteinase (MMP)-9, an ECM degradation enzyme, might play a vital role in uterine cancer progression [25]. Previous studies by Grabarek B et.al [26] also shows normal human dermal fibroblast cells along with JAK/STAT signaling pathways regulates the numerous miR’s expression under the influence of adalimumab therapy, indicating ECM molecules can be targeted for treatment effectiveness. However, the regulatory mechanisms of ECM and ECM-associated molecules in EC are still unclear.

The ECM and the ECM associated molecules involved in the microenvironment surrounding the cells seem to play important roles in tumor development; however, a comprehensive analysis of the regulatory mechanisms of ECM during EC development is still lacking. In this study, we attempted to establish a PAN-ECM Data List by collecting ECM and ECM-related genes through a literature survey and incorporation of related databases. Subsequently, we retrieved the gene expression profiles of EC samples at different stages (normal, early, and late) from The Cancer Genome Atlas-Uterine Corpus Endometrial Carcinoma (TCGA-UCEC) data, and used different bioinformatic analysis tools including DESeq2, FunRich, Ingenuity Pathway Analysis (IPA), and Weighted correlation network analysis (WGCNA) to identify and verify the key ECMs and their associates that may participate EC progression. A flow chart and an overview of the computational analysis used in this study was shown in Fig 2.

Fig 2

An overview and flowchart of the study.

This study includes three main parts, (a) construction of the PAN-ECM Data List, (b) Identification of the ECMs, ECM-associateds molecules, the involved signaling pathways participating in EC development, and (c) Networking analyses for verification of the key ECM genes involved in EC progression.

Materials and methods

Construction of the PAN-ECM Data List

Through a literature survey (S1 Table) and incorporation of data from related databases (MatrisomeDB; http://matrisomeproject.mit.edu/) [27], a complete PAN-ECM Data List was generated for further analysis. The Matrisome project provides information on ECM proteins which are collected by computational identification methods, and a proteomics pipeline was developed to characterize compositions of ECM tissues and tumors [28]. However, several essential molecules were lacking, such as ECM-associated regulatory factors (e.g., SMC3), ECM synthetic/degradation enzymes (e.g., UXS1, XYLT1, and XYLT2), cell adhesion molecules (e.g., CDH1 and CD44), ECM microenvironment-associated molecules (e.g., SMAD1 and SMAD7), etc. that were previously identified as playing critical roles in ECM regulation [29-35].

Therefore, we herein focused on developing a complete dataset for ECM research and called it the PAN-ECM Data List. The construction pipeline is depicted in Fig 2A. ECM-associated genes were collected from a literature search using such keywords as cell adhesion molecules, extracellular matrix protein, proteoglycan, collagen, glycoprotein, ECM regulator, ECM micro-environment associated, proteoglycans and secreted factors, etc. Additionally, ECM genes in MatrisomeDB were also integrated into the PAN-ECM Data List, and their functions were annotated from functional annotation databases including Gencard.

EC data collection and analysis of differentially expressed genes (DEGs)

Inclusion and exclusion criteria for UCEC sample and the classification of DEGs

The RNA-Seq level 3-gene expression data for normal and solid tumors from UCEC patients were collected from TCGA [36] (Broad GDAC Firehose database: http://firebrowse.org/?cohort=UCEC). In total, 35 normal samples, and 476 tumor samples were used for this study. Based on FIGO staging for endometrial carcinoma, the tumor samples at stage IA was classified as early-stage samples, samples obtained from patients at stage II and above were classified as late-stage samples. In this study, 159 samples were late-stage and 317 samples were early-stage samples. Based on the description of the pathology in FIGO staging, EC tumors at stage IB has ≥ 50% invasion of the myometrium. To avoid the interferences, the data from the samples at stage IB patient was excluded in this study.

Both raw read counts and normalized fragments per kilobase of transcripts per million mapped reads (FPKM) values (FPKM-UQ) at the gene level were used for further analysis. Clinical XML files were also downloaded to retrieve stage and related clinical information, which was used to recognize tumorigenesis and metastasis samples, as explained in FIGO guidelines [21, 22]. Tumorigenesis DEGs were identified by comparison the DEGs from early-stage (stage IA) to normal samples, and metastatic DEGs were identified by comparing the DEGs from late-stage (stage II or above) to early-stage EC samples (stage IA).

Statistical analysis

Expression data were pre-processed and analyzed using DESeq2, an R package [37], and log2 fold change cutoffs of >1.5 and adjusted p values of <0.01 were used to determine significant DEGs. With thousands of genes tested in an RNA-Seq, multiple comparison adjustments were necessary. The Bonferroni method was applied for filtering DEGs; this controls the mean number of false positives, that can be used for multiplicity adjustment [38]. Samples at different stages were taken for a comparative analysis to identify tumorigenesis- and metastasis-associated DEGs during EC development, as illustrated in Fig 2B. A Chi-squared test was applied to test whether ECM genes were significantly associated with EC tumorigenesis or metastasis.

Analysis of ECM regulatory mechanisms during EC development

To identify potential functions and regulatory roles of the ECM during EC development, as shown in Fig 2C, three different network analytical approaches were applied: a canonical pathway network analysis, molecular interaction network analysis, and co-expression network analysis. The canonical pathway network analysis was performed by using an Ingenuity Pathway Analysis (IPA) [39], that uses the popular activation z-score analytical method. It was proposed by Kraèmer et al. in 2014 [39], and it measures activation states, (either increased or decreased) of pathways affected by DEGs. We used a statistical approach to define a quantitative z-score, which determines whether a biological function has significantly more “increased” predictions than “decreased” predictions (z-score >0) or vice versa (z-score <0). In general practice, an absolute z-score of >2 or <-2 may be significant. A molecular interaction network analysis was performed using FunRich [40] to generate networks based on predicted or experimentally validated interactions, such as protein-protein interactions and gene-gene regulation to identify critical ECM-involved networks or critical ECM hub genes. The coexpression network analysis was performed using WGCNA, an R package [41], for EC samples. The BlockwiseModules function was applied to construct gene coexpression hierarchical clustering, and dynamic tree cut methods were used to identify co-expression networks.

Results

PAN-ECM Data List

Based on a literature survey, a list of ECM and ECM-associated molecules including 1516 ECM-related genes, called as PAN-ECM Data List was generated (S1 Table). The PAN-ECM genes were then classified into 15 categories according to their biological functions and localizations: cell surface receptors, cytoskeleton, cell adhesion molecules, collagens, ECM receptor cofactors, ECM receptors, ECM, ECM-associated regulatory factors, ECM synthetic/degradation enzymes, ECM-affiliated proteins, ECM glycoproteins, ECM regulators, ECM microenvironment-associated, proteoglycans, and secreted factors. Numbers of genes in different categories of the PAN-ECM Data List are given in Fig 3.

Fig 3

PAN-ECM Data List description generated by a literature search.

Bar chart indicated the distribution of total 1516 PAN-ECM genes in the 15 classified categories. Numbers of the PAN-ECM genes in each category are shown at the top of each bar.

Tumorigenesis and metastasis DEG analysis

Total of 20,531 DEGs from EC patients was retrieved from TCGA. Background distribution of 1516 PAN-ECM genes in a total of 20,531 DEGs is 7.38%. Based on FIGO staging for endometrial carcinoma development, we investigated the expression level of the DEGs and classified the DEGs from normal endometrium as the normal group, the DEGs from EC samples at stage IA as an early-stage group, DEGs from EC samples at stage II and the above as late-stage group. By comparing the expression level, 7880 DEGs were considered to be tumorigenesis-associated (early vs. normal) and 1677 DEGs were metastasis-associated (late vs. early), respectively. In those tumorigeneses- and metastasis-associated DEGs, 831 (S2 Table) and 241 (S3 Table) DEGs were PAN-ECM genes, respectively. Interestingly, the 831 tumorigenesis-associated PAN-ECM genes of 7880 DEGs (10.54%, Chi-squared p<1.571e-15) and the 241 metastasis-associated PAN-ECM genes of 1677 DEGs (14.37%, Chi-square p<2.2e-16) were significant, as compared to the background distribution of 7.38% in the total of 20,531 DEGs. Detailed information about fold change values in correlation with tumorigenesis and metastasis and the list of 140 intersecting PAN-ECM DEGs were shown in Table 1. These PAN-ECM genes were separated into four different categories according to their fold changes as up/up (37 genes), up/down (35 genes), down/up (57 genes), and down/down (11 genes) in tumorigenesis and metastasis, respectively. It indicates that 48 PAN-ECM genes were either up/up or down/down throughout all stages during EC development, whereas 92 PAN-ECM genes were expressed in contrasting direction of fold change (up/down or down/up) at tumorigenesis and metastasis stages.

Table 1

The expression of 140 intersecting PAN-ECM DEGs in EC at tumorigenesis and metastasis stages.

Up/Up Genes (37#) (FC_T*, FC_M†)	Up/Down Genes (35#) (FC_T*, FC_M†)	Down/Up Genes (57#) (FC_T*, FC_M†)		Down/Down Genes (11#) (FC_T*, FC_M†)
ATP1A3 (1.30, 1.17)	ADAM28 (1.16, -0.72)	ADAM33 (-2.54, 1.05)	LOXL4 (-1.91, 0.59)	ADAMTSL3 (-1.97, -0.69)
C1QL1 (1.62, 1.27)	ADAMTS6 (1.41, -1.09)	ADAMTS1 (-2.36, 0.61)	MAG (-4.30, 1.48)	CLEC3B (-3.98, -0.62)
C1QL4 (4.78, 0.90)	ADAMTSL2 (1.41, -0.70)	ADAMTS3 (-2.22, 1.41)	MEGF10 (-1.87, 1.03)	DPT (-6.28, -1.04)
CAMK2N2 (1.76, 0.77)	BMPR1B (1.04, -0.59)	ADAMTSL4 (-2.36, 0.59)	NCAM1 (-4.28, 1.21)	MASP1 (-4.75, -0.62)
CCL7 (1.53, 0.84)	CCL24 (2.69, -1.37)	AGTR1 (-4.45, 1.49)	NELL2 (-1.89, 0.92)	MATN2 (-1.70, -0.68)
COL10A1 (1.42, 1.42)	CELA3A (3.42, -1.07)	ANGPT1 (-2.62, 1.06)	NTF4 (-3.87, 1.82)	NDP (-1.47, -0.68)
COL11A1 (2.44, 0.99)	CELA3B (2.05, -0.88)	ANGPTL7 (-2.52, 1.30)	PAPLN (-0.84, 0.85)	OGN (-5.52, -0.98)
COL11A2 (1.46, 0.98)	CST1 (6.09, -1.15)	CADM3 (-3.20, 1.09)	PCOLCE2 (-2.53, 0.83)	OVGP1 (-2.43, -1.81)
COL26A1 (4.72, 0.88)	CST4 (6.22, -1.96)	CCL11 (-1.74, 0.82)	POSTN (-2.72, 0.60)	PLG (-1.70, -1.64)
COL9A1 (3.51, 1.90)	CTSV (5.15, -0.77)	CHRD (-1.70, 0.94)	PRL (-1.52, 0.96)	SFRP4 (-2.05, -0.94)
COMP (1.85, 2.00)	DMBT1 (3.05, -1.43)	CHSY3 (-1.18, 0.63)	RELN (-1.98, 1.13)	VWA3B (-1.52, -0.96)
CRLF1 (1.32, 0.63)	FBLN1 (1.13, -0.68)	COL19A1 (-1.18, 0.62)	SCG2 (-1.13, 0.94)
DCSTAMP (1.28, 0.60)	GAD1 (5.24, -0.89)	COL20A1 (-2.10, 1.51)	SEMA6D (-2.26, 0.73)
EGFL6 (1.92, 0.66)	GDF5 (3.84, -1.05)	COL4A4 (-1.88, 1.67)	SRPX (-3.80, 0.82)
ERBB2 (0.60, 0.73)	HHIP (1.68, -0.62)	COL6A3 (-2.20, 0.80)	SYT1 (-3.47, 0.77)
EREG (2.84, 0.84)	IHH (2.89, -1.52)	COL6A6 (-2.89, 1.26)	TGM1 (-1.48, 0.72)
F2 (2.89, 0.77)	IL19 (4.92, -0.89)	COL8A2 (-1.13, 0.91)	THBS2 (-2.14, 0.65)
FGFR4 (1.83, 0.61)	ITIH2 (2.00, -0.86)	CRHBP (-4.00, 0.64)	TLL1 (-2.59, 0.63)
FZD9 (1.49, 0.81)	LMAN1L (3.55, -1.15)	DNM3 (-1.10, 0.67)	TPO (-2.49, 0.94)
IL11 (2.52, 1.21)	MATN1 (1.04, -0.65)	EYS (-0.79, 0.67)	WFIKKN2 (-2.37, 0.92)
INHBE (1.30, 0.62)	MMP26 (2.47, -3.07)	FAP (-2.45, 0.69)
ITGB6 (2.12, 0.92)	MUC13 (2.62, -0.88)	FN1 (-1.05, 0.79)
KCP (1.75, 0.70)	MUC5AC (6.62, -0.90)	FREM1 (-2.34, 0.62)
LIF (1.46, 1.24)	MUC5B (4.57, -1.13)	FSTL3 (-1.43, 0.60)
MATN4 (1.92, 0.62)	PF4V1 (4.66, -1.07)	GDF6 (-2.56, 1.05)
MMP1 (4.57, 1.05)	PLA2G10 (2.46, -0.72)	GREM1 (-2.92, 0.77)
MMP10 (3.71, 1.14)	BDNF (1.63, -1.76)	HGF (-2.38, 0.60)
MMP13 (4.07, 1.12)	S100A3 (2.23, -0.61)	HPSE (-1.08, 0.66)
PRSS1 (3.87, 1.34)	S100A7 (5.31, -1.33)	IGSF10 (-2.65, 0.69)
PRSS3 (2.75, 1.01)	SEMA3E (2.18, -0.84)	IL5 (-1.53, 1.00)
RSPO4 (1.73, 1.28)	SERPINA11 (4.81, -1.87)	IL6 (-2.88, 1.51)
SFTPB (2.76, 1.08)	NTF3 (1.43, -1.18)	IMPG2 (-3.78, 0.66)
SLIT1 (1.25, 0.59)	TNFSF14 (1.29, -0.79)	ISM2 (-1.86, 0.78)
SST (3.72, 2.93)	TPH1 (1.79, -0.84)	ITGB3 (-1.54, 1.02)
TMPRSS15 (2.33, 0.99)	WIF1 (2.75, -1.06)	ITLN1 (-1.70, 0.87)
TNF (2.55, 0.61)		LAMA2 (-3.10, 0.68)
WNT7A (3.64, 0.97)		LGI2 (-3.55, 1.07)

Total of 140 intersecting PAN-ECM genes were identified and grouped as Up/Up, Up/Down, Down/Up, Down/Down expression in EC at tumorigenesis and metastasis stages respectively.

#Number of genes.

*The fold change of DEGs in tumorigenesis (FC_T).

†the fold change of DEGs in metastasis (FC_M).

Pathway analysis of tumorigenesis- and metastasis-related PAN-ECM genes

IPA was applied for pathway analysis for the identified EC-associated PAN-ECM genes (Fig 4A). Numerous signaling pathways associated with the140 intersecting PAN-ECM DEGs were identified (Fig 4B). The identified pathways included GP6 signaling, IL-8 signaling, ILK signaling, neuroinflammation signaling pathway, colorectal cancer metastasis signaling, leukocyte extravasation signaling, and osteoarthritis pathway.

Fig 4

Diagram of the intersecting PAN-ECM genes and their involved canonical pathways in association with tumorigenesis and metastasis of EC by IPA analysis.

(a) The diagram shows the PAN-ECM genes in the intersection between tumorigenesis and metastasis. (b) The heat map represents the intersecting PAN-ECM gene involved canonical pathways and their expressions in tumorigenesis and metastasis of EC. The orange color indicates activation of functions/pathways, blue color indicates inhibition of functions/pathways. Intensive color indicates higher absolute z-score.

Interestingly, several signalling pathways were expressed in opposite fold changes at tumorigenesis and metastasis stages. GP6 signaling pathway was inhibited with the lowest z-score (-0.832) in tumorigenesis but activated in metastasis with the highest z-score (3.605). IL-8 signaling (z-score = -1), ILK signaling (z-score = -1), and neuroinflammation signaling pathways (z-score = -0.447) were also inhibited at tumorigenesis, but activated in metastasis. Moreover, pathways such as colorectal cancer metastasis signaling (z-score = 1.133 in tumorigenesis and 1.889 in metastasis), leukocyte extravasation signaling (z-score = 1 and 1), and osteoarthritis (z-scores = 0.447 and 1.341) were activated at both tumorigenesis and metastasis stages.

Pathway analyses were also carried out for the identified tumorigenesis and metastasis-associated PAN-ECM genes, detailed canonical pathways generated from 831 tumorigenesis- and 241 metastasis-associated PAN-ECM genes were shown in S1 Fig. Notably, GP6 signaling pathway had the highest z-score (4.14) in the metastasis-associated PAN-ECM genes and the lowest z-score (-1.414) in the tumorigenesis-associated PAN-ECM genes, that resembles the results from the analyses for the intersecting PAN-ECM genes. The PAN-ECM genes involved in the GP6 signaling pathway were listed in Fig 5. The PAN-ECM genes including COL4A4, COL6A6, ITGB3, PI3K, and αIIbβ3 in GP6 signaling pathway were upregulated in metastasis and downregulated in tumorigenesis.

Fig 5

Tumorigenesis- and metastasis-associated PAN-ECM genes and their expressions inGP6 signaling pathway which were generated using QIAGEN’s Ingenuity® Pathway Analysis (IPA®, QIAGEN Redwood City) Software (27).

Typical EC related tumorigenesis- and metastasis-associated PAN-ECM genes were identified. Most of them were collage-related genes and were downregulated in tumorigenesis and upregulated in metastasis of EC. The PAN-ECM genes in volved in GP6 signaling pathway with opposite expression in tumorigenesis and metastasis of EC were colored. Green colour genes indicates an inhibition z-score of -1.414. Red colour genes indicates an activation z-score of 4.146; The color is more intensive, the absolute z-score is higher.

Functional enrichment analysis of tumorigenesis- and metastasis-associated DEGs

To depict the possible roles of the intersecting PAN-ECM genes in tumorigenesis and metastasis of EC, the disease-related and bio functions of the DEGs from IPA were enumerated in Table 2. Interestingly, the most correlated bio functions, cell movement, colony formation, cell proliferation, cell migration, adhesion of immune cells, and invasion of cells were activated with higher z-scores (>2) in metastasis and inhibited in tumorigenesis with lower z-scores (<-2). The bio functions, such as cell death and survival, cell-to-cell signaling, proliferation of epithelial cells, cellular growth, and cell cycle mitogenesis, were activated in both metastasis and tumorigenesis of EC. It suggested that the differential activation and inhibition of the bio functions of the PAN-ECM genes may lead to EC progression.

Table 2

Diseases and bio-functions of the intersecting PAN-ECM genes associated with differential expressions in tumorigenesis and metastasis of EC.

Diseases and bio-functions	Tumorigenesis (z-score)	Metastasis (z-score)
Cardiovascular system development and function: cell movement of endothelial cells	-2.116	0.761
Cellular growth and proliferation: colony formation of cells	-1.377	2.021
Cell death and survival, hematological system development and function: cell viability of myeloid cells	-1.217	2.177
Cellular movement, hair and skin development and function: cell movement of epithelial cell lines	-1.05	2.418
Cellular movement: migration of cells	-0.755	3.864
Cell-to-cell signaling and interaction, hematological system development and function, immune cell trafficking: adhesion of immune cells	-0.487	2.381
Cellular development, cellular growth and proliferation: proliferation of myeloma cell lines	-0.478	2.146
Free radical scavenging, small molecule biochemistry: biosynthesis of hydrogen peroxide	-0.462	2.202
Cellular development, cellular growth and proliferation: proliferation of prostate cancer cell lines	-0.349	2.195
Post-translational modification: phosphorylation of proteins	-0.34	3.06
Cell-to-cell signaling and interaction, cellular growth and proliferation: induction of cells	-0.232	2.201
Proliferation of connective tissue cells	-0.194	2.391
Cellular growth and proliferation, connective tissue development and function, tissue development	-0.128	2.402
Lipid metabolism, small molecule biochemistry: synthesis of leukotriene	-0.111	2.352
Cellular movement: invasion of cells	0	2.138
Cell-to-cell signaling and interaction, hematological system development and function: binding of lymphocytes	0.068	3.174
Cellular movement: invasion of tumor cell lines	0.113	2.72
Cell-to-cell signaling and interaction, hematological system development and function: binding of mononuclear leukocytes	0.119	2.584
Cell-to-cell signaling and interaction: binding of myeloid cells	0.221	2.697
Cellular growth and proliferation, tissue development: proliferation of epithelial cells	0.528	4.022
Cellular movement: migration of ovarian cancer cell lines	0.577	2.194
Cell signaling, molecular transport, small molecule biochemistry, vitamin and mineral metabolism: release of Ca2+	0.603	3.046
Cellular development, cellular growth and proliferation: proliferation of stem cells	0.672	2.167
Cell-to-cell signaling and interaction: binding of lymphoma cell lines	0.757	2.366
Cell signaling, molecular transport, vitamin and mineral metabolism: quantity of Ca2+	1	2.387
Cell death and survival: apoptosis of gonadal cells	1.076	2.19
Cellular growth and proliferation: expansion of cells	1.258	2.01

Activated pathways are coloured orange (z-score >2) and inhibited pathways are coloured blue (z-score <-2). The absolute z-score is higher, the color is more intensive.

Network analysis of intersecting PAN-ECM genes between tumorigenesis and metastasis DEGs

To decipher the possible regulatory mechanisms of the ECM involved microenvironment during EC development, a molecular interaction network was generated by applying the 140 intersecting PAN-ECM genes using FunRich analysis. The generated network is shown in Fig 6. Each network consists of connected PAN-ECM genes with a hub gene. As indicated, hub genes with the most intensive connections were FN1, ITGB3, ERBB2, and BMPR1B. In addition, the hub genes with the contrasting fold change in tumorigenesis and metastasis included 100A7, FBLN1, SERPINAS, TNFSF14, BMPR18, and GDF5, which were upregulated in tumorigenesis and downregulated in metastasis, and FSTL3, ADAMTS1. COL4A4, HGF, MAG, FN1, ITGB3, NCAM1, AGTR1, and DNM3, which were downregulated and upregulated in tumorigenesis and metastasis, respectively.

Fig 6

Molecular network of PAN-ECM genes that intersected in tumorigenesis and metastasis.

The molecular interaction network of the intersecting PAN-ECM genes are shown. Upregulated and downregulated genes are shown in red and green, respectively. Some genes, such as S100A7, FBLN1, SERPINAS, TNFSF14,BMPR18, and GDF5 were upregulated in tumorigenesis and downregulated in metastasis. In contrast, FSTL3, ADAMTS1. COL4A4, HGF, MAG, FN1, ITGB3, NCAM1, AGTR1, and DNM3 were downregulated an upregulated in tumorigenesis and metastasis respectively.

Gene coexpression network analysis

The gene coexpression network was constructed using the retrieved 20,531 genes from EC tumor samples using a systematic unsupervised method, WGCNA. As shown in Fig 7A, numerous modules were generated, whereas four modules consisting of the highest proportions of PAN-ECM genes, ranging 11%~52%, were selected for further analyses Fig 7B. Significantly involved pathways among the modules are enumerated in Table 3.

Fig 7

Gene coexpression network during endometrial cancer (EC) development.

(a) A WGCNA study for gene coexpression network analysis between late and early stages of EC disease was carried out to determine the relationships between PAN-ECM genes and the other co-expressed differential exressed genes (DEGs) during EC development was varried out. (b) Percentage of PAN-ECM genes in each generated WGCNA module was indicated. (c) The identified hub genes and their molecular interaction of the molecular network of the genes in module 1 by FunRich analysis.

Table 3

Gene clustering of the metastasis-associated PAN-ECM differentially expressed genes (DEGs) by WGCNA module analysis.

Canonical pathway (IPA Knowledgebase database)	Metastasis PAN-ECM (241+) (z-score)	Module 1 (77) (z-score+)	Module 2 (74) (z-score+)	Module 3 (164) (z-score+)	Module 4 (253) (z-score+)
GP6 signaling pathway	4.15	2.11	1.41	0.00	0.00
Dendritic cell maturation	2.45	2.24	2.45	0.00	0.90
Osteoarthritis pathway	2.31	1.34	0.45
Colorectal cancer metastasis signaling	2.89		2.14
Wnt/Î2-catenin signaling	2.45		-2.00		0.00
Neuroinflammation signaling pathway	2.45	0.00	1.73	0.00	0.30
ILK signaling	2.45		1.89
Role of pattern recognition receptors in recognition of bacteria and viruses	2.24	0.00	2.45	0.00	0.38
Neuregulin signaling	2.24	0.00	2.24	0.00
Signaling by Rho family GTPases	2.00	0.00	2.53		0.00
IL-6 signaling	2.00		3.00
ErbB signaling	2.00	0.00	2.71	0.00
Growth hormone signaling	2.00	0.00	2.65
Glioblastoma multiform signaling	1.89		2.71
Leukocyte extravasation signaling	1.63		2.53		-1.67
HMGB1 signaling	1.41	0.00	2.83
IL-8 signaling	1.34	0.00	3.46	-1.00	-0.71
Role of NFAT in cardiac hypertrophy	1.34		2.71
eNOS signaling	1.00	0.00	2.65
FGF signaling	0.82	0.00	2.65	0.00
Acute-phase response signaling	0.45		2.45
Gα12/13 signaling	-0.45	0.00	3.00
RhoGDI signaling	-1.00		-1.34
LXR/RXR activation	-1.00	0.00	0.00	0.00	0.45

The canonical pathways involved in each modules was listed and sorted by z-score.

Number of genes.

As described previously, high grade and re-occurred EC have poor prognosis outcome. Therefore, we attempted to dig the possible mechanisms and the characteristics of the ECM microenvironment during high-grade EC at metastasis stage. The pathway analyses for metastasis-associated PAN-ECM genes listed in Table 3 showed that module 1 and 2 possessed the highest proportions (52% and 30% respectively) of PAN-ECM genes, where those PAN-ECM genes are involved in three pathways, GP6 signaling, dendritic cell maturation, and osteoarthritis pathways, suggesting that those PAN-ECM genes and their involved pathways may play crucial roles during EC metastasis.

In addition, colorectal cancer metastasis signalling, IL-6 signaling, and IL-8 signaling also were also predicted to be activated during EC metastasis, whereas RhoGDI signaling pathway was predicted to be inhibited in metastatic EC. The molecular network of coexpression module 1 containing the highest proportion of PAN-ECM genes including COL1A1, BGN, NID1, ADAM12, COL1A2, INHBA, THBS1, SPARC, FBN1, and CSPG4 was shown in Fig 7C. The typical metastasis-associated hub genes were identified.

Discussion

In this study, PAN-ECM Data List was generated, that covers not only ECMs, but also synthetic and degradation enzyme for ECMs and ECM associated molecules and the regulators, also present in the extracellular microenvironment surrounding cells. The PAN-ECM Data List has a total of 1516 ECMs and ECM-associated genes. Comparison to the Matrisome Project [28], there were 489 additional genes that contribute in 9 categories including cell surface receptor, cytoskeleton, cell adhesion molecules, ECM receptor cofactors, ECM receptors, ECM membrane-associated, ECM-associated regulatory factors, ECM synthetic/degradation enzymes, and ECM microenvironment associated. Some of the genes were reported to be involved in many critical biological processes, for instance, integrin subunit alpha 2b (ITGA) and intercellular adhesion molecule 1 (ICAM1) genes, which is an ECM receptor and a molecule that can bind to integrin, respectively, both play crucial roles in tumor development [42]. In addition, the cell surface receptor category includes transforming growth factor-beta receptor genes (TGFBR1 and TGFBR3), and that the elevated TGFBR levels could induce collagen expression in metastatic breast cancer cells [43]. Several genes in the cell adhesion molecule category, such as CD44 and CADM1, were involved in cancer cell migration and the regulation of ECM adhesion [44, 45]. ITGB3 and BMPR2 genes were key regulators of cell migration and invasion of cancer cells [46, 47]. XYLT1 in ECM synthetic/degradation enzymes play important roles in the initiation step of the biosynthesis of glycosaminoglycan [48]. Differential expression of XDH gene modulated the migration of human breast cancer cells [49].

One forty ECM and ECM-associated genes have been identified, that may participate in tumorigenesis and metastasis of EC. Numerous identified genes have been reported to participate in the development of other cancers, for instance, COL11A1, MMP1, COL10A1, HHIP, and COL6A6 could be prognostic and predictive indicators of early-stage non-small cell lung cancer (NSCLC) [50]; upregulation of ITBG3 and COL4A3 genes could respectively promote colorectal cancer and breast cancer development [51, 52]. Alterations of EGFL6 associated with metastasis of ovarian cancer [53].

Interestingly, several PAN-ECM genes involved signaling pathways, GP6 signaling pathway, ILK signaling pathway, and IL-8 signaling pathway were significantly inhibited in tumorigenesis and activated in metastasis of EC, suggesting that those specific PAN-ECM genes may participate in EC metastasis and progression. To our knowledge, PAN-ECM gene involved GP6 signaling in association with EC progression has not been reported. GP6, a glycoprotein, is mainly expressed in platelets. GP6 signaling pathway was primarily thought to participate in platelets and their precursor megakaryocyte activation [54]. Aggregation of platelets may shield circulating tumor cells from immunosurveillance and to promote tumor progression and metastasis [55]. Another report suggested that GP6 served as the primary signaling receptor for collagen, which induces platelet activation and thrombus formation [54]. Defect of GP6 could cause malfunction of platelet reactivity to collagens in blood and lung cancer metastasis [56, 57]. An elevated level of serum platelets lowered the chemotherapy response rate and increased the risk of recurrence in EC [58]. Our current study identified several collagen genes that were upregulated in metastatic EC, suggesting that the upregulated collagens may activate GP6 signaling and protect the EC cells from immunosurveillance and promote tumor metastasis.

ILK signaling pathway was shown to represent an essential kinase activity and also regulated the decidualization of endometrial stromal cells [59]. It also played a role in hormone-dependent cancer progression and regulated a variety of cellular functions including cell survival, migration, and angiogenesis, its gain or loss of function resulted in oncogenic transformation and progression of invasive and metastatic phenotypes [60]. IL-8 signaling belongs to CXC chemokine family. Many studies confirmed high levels of IL-8 in HER2-enriched and basal-like (ER–) primary breast cancer, whereas patients with a low level of IL8 had better prognosis [61].

The hub gene they interact with the other PAN-ECM genes may play crucial roles in regulating the EC development and the progression. From the molecular network analysis (Fig 6), we identified several hub genes, such as FN1, ITGB3, and BMPR1B, interacted with the other ECM molecules intensively. These genes were associated with many critical biological functions related to cancer cell metastasis. FN1 was previously reported to be associated with adverse or poor prognosis of breast and ovarian cancer [62]. Increased expression of ITGB3 could promote and control the metastasis of breast cancer and EC [63, 64]. BMPR1B was correlated with cell proliferation, suggesting that it may play a role in protecting against endometriosis progression [65]. A previous study by Frida et al.'s team explained the transformation of normal cells to malignant and metastatic tumor cells. VCAN, an ECM gene was initially downregulated but upregulated in metastasizing cells in the prostate and colorectal cancers [66]. Our results showed that the metastatic EC associated ECM genes were highly expressed at a late stage and might participate in EC progression.

From the coexpression network analysis identified several modules possessing higher proportions of ECM genes. Interestingly, the GP6 signaling pathway along with the dendritic cell maturation and osteoarthritis pathways were activated in modules 1 and 2 and in metastatic DEGs (Table 3). Notably, GP6 signaling was identified not only by functional enrichment analysis but by coexpression network analysis. PAN-ECM genes, such as COL1A1, COL1A2, BGN, NID1, ADAM12, INHBA, THBS1, SPARC, and CSPG4, were identified to act as hub genes. Increased expressions of COL1A1 and COL1A2 genes were associated with lower survival of ovarian cancer patients [67]. BGN gene could enhance migration and invasion of EC [68]. NID1 gene was found to be a new therapeutic target biomarker in breast cancer and EC invasion [69]. ADAM12 was a tumor marker and has negative prognostic value for overall survival in ovarian cancer [70]. Strong colocalization of INHBA was highly associated with malignant endometrial tissues [71]. THBS1 expression might associate with the survival of EC patients [72]. Overexpression of SPARC gene could promote migration of EC [73]. Increased expression of CSPG4 gene was correlated with disease reappearance in patients with breast cancer [74].

Increased matrix stiffness has a profound effect on tumor development and progression [75]. Upregulation of collagen and integrins is known to participate in stiffness [76, 77]. Increased matrix stiffness resulted in increases in cell outgrowth and sprouting, which were correlated with cancer metastasis [78]. Collagen expressions, mainly types I, III, and VI, have been studied in the past few years, and higher expressions of those ECM genes were positively correlated with matrix stiffness [27, 79–81]. In this study, many collagen-associated genes, such as COL4A4, COL19A1, and COL6A6 (Table 1), showed decreased expressions at early stages and increased expressions at the late stage of EC. These findings suggested that those ECM genes may participate in EC progression.

Conclusions

In this work, we compiled a comprehensive PAN-ECM Data List including ECM and ECM-associated genes. The ECM genes including collagens and integrins and the signaling pathway, GP6, ILK signaling, and IL-8 signaling pathways, were highly correlated with metastatic EC. The identified PAN-ECM genes and their participated signaling pathways may lead to the development of novel biomarkers and new therapeutic strategies for the re-occurred, poor prognosis of high-grade EC.

Ingenuity Pathway Analysis (IPA) of tumorigenesis- and metastasis-associated PAN-ECM gene involved canonical pathways.

(a) List of tumorigenesis-associated PAN-ECM gene involved canonical pathways (b) List of metastasis-associated PAN-ECM gene involved canonical pathways. Orange color indicates activated pathways with positive z-scores. Blue color indicates inhibited pathways with negative z-scores. The color is darker indicate the absolute z-score is higher.

(DOCX)

Click here for additional data file.

Complete PAN-ECM Data List with PUBMED ID (1516).

(DOCX)

Click here for additional data file.

Tumorigenesis PAN-ECM DEGs (831 genes).

(DOCX)

Click here for additional data file.

Metastasis PAN-ECM DEGs (241).

(DOCX)

Click here for additional data file.