Identification of endometriosis-associated genes and pathways based on bioinformatic analysis.

ABSTRACT: Endometriosis is associated with dysmenorrhea, chronic pelvic pain, and infertility. The specific mechanism of endometriosis remains unclear. The aim of this study was to apply a bioinformatics approach to reveal related pathways or genes involved in the development of endometriosis.The gene expression profiles of GSE25628, GSE5108, and GSE7305 were downloaded from the gene expression omnibus (GEO) database. Differentially expressed gene (DEG) analysis was performed using GEO2R. The database for annotation, visualization, and integrated discovery (DAVID) was utilized to analyze the functional enrichment, gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) pathway of the differentially expressed genes. A protein-protein interaction (PPI) network was constructed and module analysis was performed using search tool for the retrieval of interacting genes and cytoscape.A total of 119 common differentially expressed genes were extracted, consisting of 51 downregulated genes and 68 upregulated genes. The enriched functions and pathways of the DEGs and hub genes include DNA strand separation, cellular proliferation, degradation of the extracellular matrix, encoding of smooth muscle myosin as a major contractile protein, exiting the proliferative cycle and entering quiescence, growth regulation, and implication in a wide variety of biological processes.A bioinformatics approach combined with cell experiments in this study revealed that identifying DEGs and hub genes leads to better understanding of the molecular mechanisms underlying the progression of endometriosis, and efficient biomarkers underlying this pathway need to be further investigated.

Introduction

Endometriosis is a debilitating disease with features of chronic inflammation and is defined as the presence of functional endometrial glands and stroma outside the uterine cavity, the most common locations for the ectopic endometrial implants being the ovaries, the fossa ovarica, the uterosacral ligaments, and the posterior cul-de-sac.[ Endometriosis currently affects approximately 5.5 million reproductive-aged women in the United States.[ Worldwide, it represents a significant cause of morbidity in approximately 10% to 15% of women in their reproductive years.[ Endometriosis is one of the major causes of economic burden and compromised quality of life in a very large percentage of Asian women.[ While it is perceived as a benign condition, recent research has shown that it may be a significant cause of infertility and metastatic cancer.[

Although the cause of endometriosis remains unclear, genetic,[ hormonal, and immunological factors[ as well as endometrial progenitor cells have been implicated in the development of lesions.[ Endometriosis-associated genes and pathways still remain unclear. This present study aimed to identify critical genes and pathways contributing to endometriosis.

Materials and methods

Microarray data

The gene expression omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) is a public functional genomics data repository of high-throughput gene expression data, chips, and microarrays.[ Three gene expression datasets [GSE25628 (Platforms: GPL571, Affymetrix Human Genome U133A 2.0 Array),[ GSE5108 (Platforms: GPL2895, GE Healthcare/Amersham Biosciences CodeLink Human Whole Genome Bioarray),[ and GSE7305 (Platforms: GPL570, Affymetrix Human Genome U133 Plus 2.0 Array)][ were downloaded from GEO. The GSE25628, GSE5108, and GSE7305 datasets contained 22, 22, and 20 samples, respectively.

Identification of differentially expressed genes

Differentially expressed gene (DEG) analysis for the 3 datasets were carried out using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/). GEO2R is an interactive web tool that allows users to compare 2 or more datasets in a GEO series in order to identify DEGs across experimental conditions. For the 3 gene expression datasets, an adjusted P value of <.01 and a ≥ fourfold change were set as the cutoff to be considered a statistically significant differentially expressed gene. The common differentially expressed genes of the 3 datasets were selected using the Draw Venn Diagram web tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Kyoto encyclopedia of genes and genomes and gene ontology enrichment analyzes of differentially expressed genes

The database for annotation, visualization, and integrated discovery (DAVID; http://david.ncifcrf.gov) (version 6.8) is an online biological information database that integrates biological data and analysis tools, and provides a comprehensive set of functional annotation information of genes and proteins for users to extract biological information.[ KEGG is a database resource for understanding high-level functions and biological systems from large-scale molecular datasets generated by high-throughput experimental technologies.[ GO is a major bioinformatics tool to annotate genes and analyze the biological processes of these genes.[ In order to analyze the function of DEGs, biological analyzes were performed using the DAVID online database. P < .0001 for GO functional enrichment analysis and P < .05 for KEGG pathway enrichment analysis were considered statistically significant.

Protein-protein interaction network construction and module analysis

The PPI network was predicted using the search tool for the retrieval of interacting genes (STRING; http://string-db.org) (version 11.0) online database.[ Analyzing the functional interactions between proteins may provide insights into the mechanisms of generation or development of diseases. In the present study, a PPI network of DEGs was constructed using the search tool for the retrieval of interacting genes database, and an interaction with a combined score of >0.4 was considered statistically significant. Cytoscape (version 3.8.0) is an open-source bioinformatics software platform for visualizing molecular interaction networks.[ The molecular complex detection (MCODE) (version 1.6.1) plug-in of cytoscape is an application for clustering a given network based on topology to discover densely connected regions.[ The PPI networks were drawn using cytoscape and the most significant module in the PPI networks was identified using MCODE. The criteria for selection were as follows: MCODE score >5, degree of cut-off=2, node score cut-off=0.2, Max depth=100, and k-score=2. Subsequently, the KEGG and GO analyzes for genes in this module were performed using DAVID.

Hub gene selection and analysis

Genes with a connectivity degree ≥10 were selected as hub genes. Biological process analysis of hub genes was performed and visualized using the Biological Networks Gene Oncology tool (version 3.0.4) plugin of cytoscape.[

The Institutional Review Board or Ethics Committee approval was not needed.

Results

Overview of the gene expression omnibus microarray data and identification of differentially expressed genes

For GSE25628, GSE5108, and GSE7305, an adjusted P value of <.01 and a ≥ fourfold change were set as the cutoff to be considered a statistically significant differentially expressed gene. A total of 2514 elements in GSE25628, 2767 elements in GSE5108, and 2403 elements in GSE7305 were selected. Using Draw Venn Diagram, 2001, 2131, and 1678 unique elements were identified in GSE25628, GSE5108, and GSE7305, respectively. A total of 119 common DEGs were extracted from the 3 groups after comparing (Fig. 1), and consisted of 51 downregulated genes and 68 upregulated genes.

Figure 1

Venn diagram. Differentially expressed genes with an adjusted P value of <.01 and a ≥ four-fold change were selected from among the gene expression profiling sets GSE25628, GSE5108, and GSE7305. The 3 datasets showed an overlap of 119 genes.

Gene ontology functional enrichment analysis

Three categories of GO functional annotation analysis were performed on these potential target genes mentioned above, including biological process, cellular component, and molecular function. As shown in Table 1, the GO analysis results for the common DEGs indicated circulatory system development and regulation of the meiotic cell cycle in the biological process category; extracellular region, bicellular, and apical junction in the cellular component category; and protein binding, metallocarboxypeptidase activity, and carboxypeptidase activity in the molecular function category.

Table 1

GO analysis of the common differentially expressed genes (P < .0001).

Category	Term	Description	Count	Fold enrichment
Downregulated
GO MF	GO:0030674	Protein binding, bridging	2	33.48245614
GO CC	GO:0005923	Bicellular tight junction	4	17.74306445
	GO:0070160	Occluding junction	4	17.02375102
	GO:0043296	Apical junction complex	4	14.82067736
	GO:0072686	Mitotic spindle	3	20.9959596
GO BP	GO:0007049	Cell cycle	11	3.825756266
	GO:0000278	Mitotic cell cycle	8	5.095283318
	GO:1903047	Mitotic cell cycle process	7	5.183658447
	GO:0051726	Regulation of cell cycle	8	4.168868169
	GO:0000280	Nuclear division	6	5.914747673
	GO:0007067	Mitotic nuclear division	5	7.85630613
	GO:0048285	Organelle fission	6	5.398926655
	GO:0051445	regulation of meiotic cell cycle	3	28.13986014
Upregulated
GO MF	GO:1901681	Sulfur compound binding	8	12.56706192
	GO:0005539	Glycosaminoglycan binding	7	12.13978182
	GO:0008201	Heparin binding	6	14.13794466
	GO:0004181	Metallocarboxypeptidase activity	4	37.70118577
	GO:0004180	Carboxypeptidase activity	4	25.50374332
	GO:0008238	Exopeptidase activity	5	13.21840355
GO CC	GO:0031012	Extracellular matrix	17	10.92711213
	GO:0044421	Extracellular region part	32	2.171676543
	GO:0005576	Extracellular region	32	2.008260584
	GO:0005578	Proteinaceous extracellular matrix	8	7.509960509
	GO:0070062	Extracellular exosome	23	2.071438795
	GO:1903561	Extracellular vesicle	23	2.06223649
	GO:0043230	Extracellular organelle	23	2.061320755
	GO:0005615	Extracellular space	14	2.961868597
GO BP	GO:0072359	Circulatory system development	12	3.726522906
	GO:0072358	Cardiovascular system development	12	3.726522906
	GO:0001568	Blood vessel development	9	4.553181156

GO = gene ontology, BP = biological process, CC = cellular component, MF = molecular function.

Kyoto encyclopedia of genes and genomes pathway enrichment analysis

To further analyze the enriched pathways of these DEGs, KEGG pathway enrichment analysis was subsequently conducted. As shown in Table 2, KEGG pathway analysis revealed that the DEGs were mainly enriched in leukocyte transendothelial migration, cellular junction, vascular smooth muscle contraction, focal adhesion, malaria, phagosome, and arachidonic acid metabolism signaling pathways.

Table 2

KEGG pathway enrichment analysis of the common differentially expressed genes (P < .05).

KEGG pathway	Term	Description	Count	Fold enrichment
Downregulated
	hsa05200	Pathways in cancer	5	4.326923077
	hsa04270	Vascular smooth muscle contraction	3	8.728448276
Upregulated
	hsa04510	Focal adhesion	6	5.590062112
	hsa04512	ECM-receptor interaction	4	8.866995074
	hsa04530	Tight junction	4	8.766233766
	hsa04670	Leukocyte transendothelial migration	4	6.650246305
	hsa05144	Malaria	3	11.80758017
	hsa04145	Phagosome	4	5.672268908
	hsa00590	Arachidonic acid metabolism	3	9.484777518

KEGG = kyoto encyclopedia of genes and genomes, ECM-receptor = extracellular matrix-receptor.

A PPI network of DEGs was constructed (Fig. 2) and the most significant module was obtained using cytoscape (Fig. 3). The functional analyzes of genes involved in this module were performed using DAVID. Results showed that genes in this module were mainly enriched in glycosaminoglycan binding, sulfur compound binding, metallocarboxypeptidase activity, heparin binding, extracellular matrix-receptor interaction, vascular smooth muscle contraction, tight junction, and focal adhesion (Table 3).

Figure 2

A protein-protein interaction network of differentially expressed genes was constructed using cytoscape. Upregulated genes are marked in light red; downregulated genes are marked in light blue.

Figure 3

The most significant module from the protein-protein interaction network. (Upregulated genes are marked in light red; downregulated genes are marked in light blue).

Table 3

GO and KEGG analysis of DEGs in the most significant module (P < .01).

Pathway ID	Pathway description	Count	FDR
GO:0005539	Glycosaminoglycan binding	8	0.0675376
GO:1901681	Sulfur compound binding	8	0.1265086
GO:0004181	Metallocarboxypeptidase activity	4	0.9614553
GO:0008201	Heparin binding	6	0.9705039
hsa04512	ECM-receptor interaction	6	0.8194419
hsa04270	Vascular smooth muscle contraction	6	3.4249874
hsa04530	Tight junction	5	6.3134673
hsa04510	Focal adhesion	7	7.225584

GO = gene ontology, KEGG = kyoto encyclopedia of genes and genomes, DEGs = differentially expressed genes, FDR = false discovery rate, ECM-receptor = extracellular matrix-receptor.

A total of 5 genes with a connectivity degree ≥10 were selected as hub genes. The names, abbreviations, and functions for these hub genes are shown in Table 4, and the biological process analysis of the hub genes is shown in Figure 4.

Table 4

Functional roles of 5 hub genes.

Number	Gene symbol	Full name	Function
1	HELLS	Helicase, lymphoid specific	DNA strand separation, cellular proliferation
2	TIMP2	TIMP metallopeptidase inhibitor 2	Degradation of the extracellular matrix
3	MYH11	Myosin heavy chain 11	Encode smooth muscle myosin as a major contractile protein
4	QSOX1	Quiescin sulfhydryl oxidase 1	Exit the proliferative cycle and enter quiescence, growth regulation
5	LAMA4	Laminin subunit alpha 4	Implicate in a wide variety of biological processes including cell adhesion, differentiation, migration, signaling, neurite outgrowth and metastasis

Figure 4

The biological process analysis of hub genes. The color depth of the nodes refers to the corrected P value of ontologies. The size of the nodes refers to the number of genes that are involved in the ontologies. P < .05 was considered statistically significant.

Discussion

Endometriosis is a chronic inflammatory hormonal, immune, systemic, and heterogeneous disease defined as the presence of endometrial glands and stroma-like lesions outside of the uterus, often associated with inflammation, severe and chronic pain, and infertility.[ Diagnosis of endometriosis should be based on patient interviews, examination, and imaging.[ Lesions identified during laparoscopy are categorized as superficial peritoneal lesions, endometriomas, or deep infiltrating nodules, with a high degree of individual variability in lesion color, size, and morphology.[ Histopathological analysis requires the presence of at least 2 features - endometrial epithelium, endometrial glands, endometrial stroma, or hemosiderin-filled macrophages - for a diagnosis of endometriosis.[

The most well-accepted pathophysiological hypothesis for endometriosis is based on retrograde menstruation.[ Other hypotheses proposed include Müllerian metaplasia, lymphovascular emboli of endometrial cells, and proliferation of endometrial stem cells or bone marrow progenitors.[

It is understood that several factors are involved in the pathogenesis and progression of endometriosis, including inflammation, angiogenesis, cytokine/chemokine expression, and endocrine alterations such as steroid and steroid receptor expression.[

Angiogenesis is the formation of new blood vessels, and subsequently, is a key process in forming functional blood vessels to ectopic menstrual tissue for the establishment and maintenance of endometriotic lesions. The vascular endothelial growth factor protein family is well known for its roles in angiogenesis. A variety of rodent endometriosis models have shown that vascular endothelial growth factor levels increase in endometriosis-like lesions.[ Matrix metalloproteinases are proteases required for reorganizing existing blood vessels during budding angiogenesis.[ They play a known role in endometriosis.[

Cytokines and chemokines are emerging as key players in endometriosis pathobiology. Altered levels of a large number of cytokines and chemokines have been found in cyst fluid removed from endometriomas and chocolate cysts.[

Endometriosis is intimately associated with steroid metabolism and associated pathways, corresponding to the paramount roles estrogen receptors and progesterone receptors play in uterine biology. Many studies have shown that endometriosis is estrogen dependent and is regulated through alpha and beta estrogen receptors (ESR1 and ESR2).[

Although, endometriosis is intimately associated with interaction between inflammation and the endocrine system, which was considered the major mechanism of endometriosis, the exact etiology and pathophysiological mechanisms of endometriosis still remain unclear.

In the present study, a bioinformatics approach was applied to reveal the possible pathways and critical genes related to the development of endometriosis. As some significant biological functions were considered to be common, either in normal endometrium or the development of endometriosis, 3 gene expression datasets (GSE25628, GSE5108, and GSE7305) were analyzed to obtain DEGs. A total of 119 DEGs were identified among the 3 datasets, including 51 downregulated genes and 68 upregulated genes. GO and KEGG enrichment analyzes were performed to explore interactions between the DEGs. The genes were mainly enriched in protein binding, bridging, bicellular tight junction, cell cycle, extracellular region, metallocarboxypeptidase activity, and circulatory and cardiovascular system development.

GO and KEGG enrichment analyzes revealed that changes in the most significant modules were mainly enriched in glycosaminoglycan binding, sulfur compound binding, metallocarboxypeptidase activity, heparin binding, extracellular matrix-receptor interaction, vascular smooth muscle contraction, tight junction, and focal adhesion.

Five genes with a connectivity degree ≥10 were selected as hub genes. Among these hub genes, helicase, lymphoid-specific (HELLS) showed node degrees with 28 of them.

The gene TIMP metallopeptidase inhibitor 2 is a member of the TIMP gene family. The proteins encoded by this gene family are natural inhibitors of the matrix metalloproteinases, a group of peptidases involved in degradation of the extracellular matrix. In addition to an inhibitory role against metalloproteinases, the encoded protein has a unique role among TIMP family members in its ability to directly suppress the proliferation of endothelial cells. As a result, the encoded protein may be critical to the maintenance of tissue homeostasis by suppressing the proliferation of quiescent tissues in response to angiogenic factors, and by inhibiting protease activity in tissues undergoing remodeling of the extracellular matrix.[ Increased gene TIMP metallopeptidase inhibitor 2 expression has been reported in endometriosis patients.[

Literature retrieval results showed that a connection between endometriosis and the hub genes HELLS, myosin heavy chain 11 (MYH11), quiescin sulfhydryl oxidase 1 (QSOX1), and LAMA4 has not been widely reported.

The gene HELLS, a protein coding gene, encodes a lymphoid-specific helicase. Other helicases function in processes involving DNA strand separation, including replication, repair, recombination, and transcription.[ This protein is thought to be involved in cellular proliferation, associated with the occurrence of cancer. Cheuk-Ting Law found that HELLS, an SWI2/SNF2 chromatin remodeling enzyme, was remarkably overexpressed in hepatocellular carcinoma.[ A study by Xi Liu revealed that HELLS was significantly upregulated in colorectal cancer.[ Aside from cancers, immunodeficiency–centromeric instability–facial anomalies syndrome is caused by ATPase-defective point mutations in HELLS.[

MYH11, a protein coding gene that belongs to the myosin heavy chain family, encodes smooth muscle myosin. The gene product is a hexameric protein subunit that consists of 2 heavy chain subunits and 2 pairs of non-identical light chain subunits. It functions as a major contractile protein, converting chemical energy into mechanical energy through adenosine-triphosphate hydrolysis. MYH11 has been traditionally thought of as a specific and exclusive marker for vascular smooth muscle cells and pericytes.[ Bruce A. Corliss identified MYH11 as a marker of a subset of corneal endothelial cells.[ Recent studies have identified homozygous or compound heterozygous variants in MYH11 as a candidate gene for megacystis-microcolon-intestinal hypoperistalsis syndrome, a rare and severe disorder characterized by functional obstruction in the urinary and gastrointestinal tract.[

QSOX1, a protein coding gene that is a member of 2 long-standing gene families, encodes a protein that contains domains of thioredoxin and resolvin E1 receptor.[ Gene expression is induced as fibroblasts begin to exit the proliferative cycle and enter quiescence, suggesting that this gene plays an important role in growth regulation.[ Two transcript variants encoding 2 different isoforms have been found for this gene. Amber L Fifield concluded that overexpressed QSOX1 is a potential novel anticancer agent in tumors.[

Laminin subunit alpha 4 is a protein-encoding gene. Laminins, a family of extracellular matrix glycoproteins, are the major noncollagenous constituent of basement membranes. They have been implicated in a wide variety of biological processes, including cell adhesion, differentiation, migration, signaling, neurite outgrowth, and metastasis. Down-regulating LAMA4 expression inhibits the proliferation and migration of breast cancer,[ renal cell carcinoma,[ gastric cancer,[ and ovarian cancer.[

The biological process analysis the specific signaling pathways involved in the key genes and explore the potential molecular mechanisms by which the key genes influence endometriosis progression. High expression of hub genes mainly enriched maintenance of DNA methylation, protein thiol-disulfide exchange, methylation-dependent chromatin silencing, skeletal myofibril assembly, (skeletal/striated muscle) myosin thick filament assembly, system development, cardiac muscle fiber development, extracellular matrix assembly, centromere complex assembly, elastic fiber assembly, (centromeric) heterochromatin formation, heterochromatin organization, chromatin remodeling at centromere and other signaling pathways. (Fig. 4). Kuei-Yang Hsiao's study found that Epigenetic modifications, including DNA methylation, histone modification, and microRNA expression, are involved in the pathogenesis of endometriosis[ Mohamed G Ibrahim’ study got the conclusion that the nuclear membrane irregularities are evidence for ultramicro-trauma in adenomyosis.[ But the connections between endometriosis and the other enriched signaling pathways has not been widely reported. Further studies are needed to confirm the role of these unproved signaling pathways in endometriosis.

However, there are still some limitations in this study:

A lack of research on detailed molecular mechanisms that hub genes regulate endometriosis progression;

elated animal studies are deficient in this study.

Conclusion

In conclusion, the present study was designed to identify DEGs that may be involved in the progression of endometriosis. A total of 119 DEGs and 5 hub genes were identified. However, further studies are needed to elucidate the biological function of these genes in endometriosis.

Author contributions

Conceptualization: Ting Wang.

Data curation: Ruoan Jiang, Linhua Qian.

Formal analysis: Ruoan Jiang, Yingsha Yao.

Methodology: Ting Wang, Ruoan Jiang.

Project administration: Xiufeng Huang.

Software: Yingsha Yao, Yu Zhao.

Validation: Ting Wang.

Writing – original draft: Yingsha Yao.

Writing – review & editing: Ruoan Jiang.