Identification of differentially expressed genes and signaling pathways involved in endometriosis by integrated bioinformatics analysis

Endometriosis, a benign disease with malignant properties, is defined by the presence of active endometrial cells outside the uterus, including in the pelvis, abdominal cavity, intestines, rectovaginal septum, abdominal wall and myometrium (also known as adenomyosis) ( 1 , 2 ). Patients afflicted with endometriosis are associated with higher occurrences of anxiety and depression ( 3 ). In addition, women with endometriosis are more likely to experience non-menstrual pelvic pain (36.7% vs. 14.3%), dyspareunia (29.5% vs. 13.4%) and infertility (11.6% vs. 3.4%) compared with women without endometriosis ( 4 ). The current gold standard treatments for endometriosis are surgical resection and hormone suppression; however, none of these therapies are ideal due to their various side effects and high recurrence rates ( 5 ). In order to improve the treatment strategies, it is important to study the underlying mechanisms involved in disease development and progression. In addition, identifying reliable molecular markers can aid in the diagnosis and treatment of endometriosis. The upregulation and downregulation of genes associated with disease susceptibility serves an important role in the progression of endometriosis ( 6 ). With the development of next-generation sequencing, gene expression microarrays have been widely performed to identify the differentially expressed genes (DEGs) that may be involved in the development and progression of endometriosis ( 5 , 7 ). However, due to the different sample size, technology detection platforms and inconsistent data processing methods across different studies, the DEGs identified in previous studies are inconsistent or even contradictory. Thus, there are certain limitations in using a single gene expression profile. Integrated bioinformatics analysis has emerged as a promising tool for exploring the molecular markers and signaling pathways involved in diseases, and has previously been applied to study ovarian cancer ( 8 ), breast cancer ( 9 ) and non-small cell lung cancer ( 10 ). In the present study, four microarray expression datasets, namely GSE11691 ( 11 ), GSE23339 ( 12 ), GSE25628 ( 13 ) and GSE78851 ( 14 ) were downloaded from the Gene Expression Omnibus (GEO) database of the National Center for Biotechnology Information. A total of 58 samples, including 27 cases of endometriosis and 31 normal samples, were included in the present study. Firstly, the gene expression profiles were normalized, and the DEGs were then identified using the limma package function of R software. Subsequently, Gene Ontology (GO) enrichment analysis of DEGs was performed on Database for Annotation, Visualization and Integrated Discovery (DAVID), while the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were analyzed via the KOBAS online analysis database. Finally, a protein-protein interaction (PPI) network was constructed using the STRING online database. Cytoscape software was applied for further visualization. The current study identified key signaling pathways and potential candidate genes involved in the development of endometriosis, which may facilitate a better understanding of the underlying molecular mechanisms and provide effective targets for the diagnosis and treatment of this disease.

The gene expression profiles of the datasets GSE11691 , GSE23339 , GSE25628 and GSE78851 were normalized, as shown in Fig. 1 . According to the criteria of |log 2 FC|>1 and P<0.05, a total of 903 DEGs were identified in GSE11691 using the limma R package, including 575 upregulated and 328 downregulated genes. A total of 1,139 DEGs were identified from the GSE23339 dataset, including 608 upregulated and 531 downregulated genes. Additionally, 1,731 DEGs were identified from the GSE25628 dataset, consisting of 708 upregulated and 1,023 downregulated genes, while there was a total of 2,118 DEGs in the GSE78851 dataset, including 221 upregulated and 1,897 downregulated genes. Subsequently, the volcano plots for the identified DEGs and the cluster heatmaps of the top 200 DEGs in each dataset were constructed, and are presented in Figs. 2 and 3 , respectively. The RRA method assumes that each gene in each dataset is randomly arranged, which is widely used in integrated bioinformatics analysis ( 17 , 18 ). Through rank analysis (corrected P-value of <0.05), 275 integrated DEGs were identified. In order to obtain more reliable DEGs, genes with inconsistent upregulation and downregulation in the expression profiles were deleted. In total, 103 integrated genes were identified, including 47 upregulated and 56 downregulated genes ( Table II ). The top 20 upregulated and downregulated genes were represented on heatmaps using Sanger Box software, as shown in Fig. 4 . The GO functional analysis was divided into the BP, MF and CC categories. As displayed in Fig. 5 , the DEGs were mainly enriched in cell adhesion, cell migration, cell-cell junction and heparin binding in the GO function annotation. Furthermore, according to the KEGG pathway analysis, the DEGs were mainly involved in adherens junction and hypoxia-inducible factor (HIF)-1 signaling. A PPI network was subsequently constructed, which consisted of 54 nodes (proteins) and 62 edges (interactions), as shown in Fig. 6 . The genes showing the most significant interaction in the network were PIK3R1, ERBB3, MRPS31, HSPA5, ZFPM2, NUP88, SUZ12, MRPL39, HSPA4, GATA6, NUPL2 , and EP300 . To further investigate whether the expression of the identified genes in endometriosis tissues was consistent with the bioinformatic analysis results, the expression of certain genes in each pathway was randomly verified. HSPA5, ENO2 and TJP1 are associated with cell migration, adherens junction and the HIF-1 signaling pathway, respectively. As shown in Fig. 7 , the findings of immunohistochemical analysis verified that the expression levels of HSPA5 and TJP1 were evidently reduced in endometriosis as compared with that in normal tissues. However, ENO2 was significantly upregulated in endometriosis, which was consistent with the bioinformatics results.

It is estimated that there are 176 million women with endometriosis worldwide, and this condition seriously affects 10% of women of reproductive age ( 19 ). Chronic pelvic pain and infertility cause great physical pain and mental distress to women with endometriosis and their partners, greatly reducing the family happiness index and increasing the domestic burden ( 20 ). Accumulating evidence suggests that the endometrium of patients with endometriosis exhibits abnormal molecular expression, which gives the tissue the ability to implant, invade and develop into endometriosis lesions ( 21 , 22 ). In order to identify more stable and reliable molecular markers, the present study mapped out the genetic alterations that may be involved in the development of endometriosis by integrated bioinformatics analysis. Four gene expression profile datasets from different groups were integrated in the present study, and R software and bioinformatics analysis were used to analyze these datasets. A total of 103 DEGs were identified using the RRA analysis method, including 47 upregulated and 56 downregulated genes. Furthermore, through GO and KEGG analyses, these DEGs were found to be closely associated with cell migration, adherens junction and the HIF-1 signaling pathway. The results revealed that the DEGs associated with cell migration in endometriosis were PIK3R1, PODXL, HSPA5 and LRRC15 , while the genes IQGAP1, TJP1 and EP300 were involved in adherens junction. Notably, the DEGs EP300, ENO2 and PIK3R1 were mainly associated with the HIF-1 signaling pathway. The most widely accepted theory for the development of endometriosis is implantation and invasiveness. Accumulating studies have indicated that the degradation of extracellular matrix and the alteration of gene expression serve critical roles in the pathophysiological processes of endometriosis ( 23 , 24 ). In addition, PI3K/Akt signaling has been reported to be involved in these processes ( 25 ). Rai and Shivaji indicated that DJ-1 regulated cell proliferation, migration and invasion in endometriotic epithelial cells via the PI3K/Akt signaling pathway. In the present study, the findings demonstrated that the gene PIK3R1 was downregulated and may be involved in cell migration ( 26 ). In addition, PODXL, HSPA5 and LRRC15 may also have potential value in this process. Intercellular junctions (including tight junctions and adherens junctions) play a critical role in the endometrium. The development of endometriosis is accompanied with changes in cell-cell tight junctions ( 27 ). Extensive research has demonstrated that claudin-3, claudin-4, ZO-3, E-cadherin and α-catenin are downregulated in the ectopic endometrium as compared with their expression in the corresponding eutopic endometrium ( 28 – 30 ). In the present study, integrated bioinformatics analysis revealed that the expression levels of genes associated with the adherens junction pathway, namely IQGAP1, TJP and EP300 , were significantly reduced in endometriosis. In the last decade, researchers have indicated that the expression of HIF-1α was higher in ectopic endometriosis tissue as compared with that in eutopic tissue ( 31 , 32 ). Furthermore, hypoxia can induced the invasion of endometrial stromal cells and promoted the endometriosis-associated angiogenesis ( 33 , 34 ). Additionally, the expression of HIF-1α in the serum was reported to be proportional to the stage of endometriosis and the severity of pain ( 32 ). Indeed, bioinformatics analysis in the present study deonmonstrated that the expression of genes associated with HIF-1α, such as ENO2 , was upregulated in endometriosis. In conclusion, the present study revealed that cell migration, adherens junction and the HIF-1 signaling pathway may be involved in the development of endometriosis via integrated bioinformatics analysis. In addition, these identified DEGs may be of clinical significance for the diagnosis and treatment of the endometriosis. However, as the present study is solely based on data analysis and experimental verification, further studies with larger samples and clinical trials are required to confirm the function of the identified genes in endometriosis.

The keyword ‘endometriosis’ was used to search the GEO database ( http://www.ncbi.nlm.nih.gov/geo ), and the gene expression profiles of GSE11691 ( 11 ), GSE23339 ( 12 ), GSE25628 ( 13 ) and GSE78851 ( 14 ) were downloaded. The dataset GSE11691 , based on the platform GPL96 [HG-U133A] Affymetrix Human Genome U133A Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA), included 9 endometriosis tissues and 9 normal endometrial tissue samples. The platform for GSE23339 was GPL6102 (Illumina Human-6 v2.0 expression beadchip; Illumina, Inc., San Diego, CA, USA), including 10 endometrioma samples and 9 control endometrium specimens. GSE25628 included 7 samples of ectopic endometrioma and 6 samples of normal endometrial tissue, and its platform was GPL571 [HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc.). The platform for GSE78851 , consisting of 3 tissues from patients with adenomyosis and 5 normal tissues, was GPL6244 [HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version] (Affymetrix; Thermo Fisher Scientific, Inc.). The platform and series matrix files were downloaded as CSV files. The dataset information is displayed in Table I . The gene IDs within each gene expression profile was converted into a gene symbol, and then the data were log 2 transformed and normalized using R 5.3.1 ( https://www.r-project.org/ ). DEGs between endometriosis and non-endometriosis samples were screened out under the thresholds of |log 2 fold change (FC)|>1 and P<0.05 using the limma package in the Bioconductor 3.9 tool ( http://www.bioconductor.org/packages/release/bioc/html/limma.html ). The volcano map of the DEGs and the heatmap of the top 200 DEGs in each microarray datasets were obtained using R. SangerBox 1.0.8 ( http://sangerbox.com/ ) is a computerized and powerful software for biological information analysis, and is used as a visualization tool. The robust rank aggregation (RRA) method can be applied as a useful and general solution for gene list integration and meta-analysis in an unbiased manner, using a probabilistic model to make the algorithm parameter free and robust to outliers, noise and errors, and to assign a significance score to each gene ( 15 ). The RRA method can rank each item in each list and compare this ranking with the baseline case where all preference lists are randomly ordered. The P-value can represent the rank location, with a smaller P-value indicating a higher gene rank. In the present study, RRA in SangerBox was performed for comprehensive sorting of DEGs in the four gene expression profiles. P<0.05 was set as the threshold, and DEGs that were inconsistent across the four data sets were excluded. GO analysis ( 16 ), which is composed of biological process (BP), cellular compartment (CC) and molecular function (MF) terms, is a common method for large-scale genomic data function annotation. In order to better understand the mechanism of DEGs involved in the development of endometriosis, GO and KEGG pathway enrichment analyses were performed using the DAVID 6.8 ( https://david.ncifcrf.gov/ ) and the KOBAS 3.0 ( http://kobas.cbi.pku.edu.cn/ ) online analysis tool. P<0.05 was considered to indicate a statistically significant difference in these analyses. The STRING database ( http://string-db.org/ ) was used to identify the interacting protein pairs among DEGs with the criterion of combined score of ≥0.4. Upon removal of the isolated and partially connected nodes, a complex network of DEGs was constructed. The file of STRING interactions was downloaded and visualized with Cytoscape 3.7.0 ( https://cytoscape.org/ ). For immunohistochemical analysis, archival samples of normal endometrial and endometriosis specimens were used. The samples had been collected between May 2018 and December 2018 from patients that underwent surgery at Renmin Hospital of Wuhan University (Wuhan, China). The age of the females from which these samples were collected ranged between 20 and 40 years old. The present study was approved by the Ethics Committee of Renmin Hospital of Wuhan University, Patients and their families signed an informed consent form in advance. In short, six normal endometrial and six endometriosis specimens were confirmed by a pathologist. The tissue samples were cut into sections of 3 µm in thickness and 3 mm in diameter. Once the samples had been dewaxed, hydrated and treated with sodium citrate (pH=6), hydrogen peroxide was used to block any endogenous peroxidase activity. Immunohistochemical staining was conducted with a rabbit polyclonal primary antibody against HSPA5 (1:150; cat. no. ab108615; Abcam, Cambridge, MA, USA), TJP1 (1:150; cat. no. 21773-1-AP; Wuhan Sanying Biotechnology, Wuhan, China) and ENO2 (1:100; cat. no. ab79757; Abcam) at 4°C overnight. Subsequently, the samples were incubated with a horseradish peroxidase-conjugated goat anti-rabbit secondary antibody (1:200; cat. no. AS-1107; Aspen) at 37°C for 50 min, and a 3,3′-diaminobenzidine solution and hematoxylin were then used for staining and counterstaining at room temperate for 1 min. The integrated option density was analyzed using the ImageJ software (version 1.4.6; National Institutes of Health).