HOX cluster and their cofactors showed an altered expression pattern in eutopic and ectopic endometriosis tissues

We evaluated 84 genes from the HOX clusters A-D and their cofactors in the control, eutopic, and ectopic tissue biopsies. Figure  1 A-D and supplementary figures S 1 -S 3 show the changes in expression patterns of HOX clusters A-D and their cofactors, respectively. Table S 2 lists the genes that showed significant differential expression (DEGs) ( P  < 0.05) in the ectopic and eutopic groups compared to the control group in addition to DEGs in the ectopic compared to the eutopic group. Among the genes we evaluated in this study, 38.1% (No: 32/84) were non-DEG and about 61.9% (52/84) were differentially expressed between the ectopic, eutopic, and control groups. There were 18 genes from the DEG that were common among the studied groups (eutopic, ectopic and control tissues) (Table S 2 , Fig. 2 A,B). Table S 2 lists the various cluster of genes that had different expression patterns in the studied groups. Fig. 1 Gene expression array analysis for HOX genes. There were three biological repeats in each group. A, B, C and D in each graph show the nine data that were analyzed with SPSS software. P  < 0.05 was considered to be significant Fig. 2 A Venn diagram for DEGs ( P  < 0.05). DEGs in the ectopic samples compared to normal samples (black ellipse), eutopic samples in comparison to normal samples (green ellipse), and ectopic samples compared to eutopic samples (red ellipse). Genes that were significantly upregulated are enclosed in black, green and red rectangles according to the above mentioned legend. In the multi-colored rectangles, each color represents upregulation in the relevant analysis in the corresponding agreement. B Heat map for all genes. Right column is logFC of ectopic samples compared to normal samples and the left column is the logFC of eutopic samples in comparison to normal samples. logFC: Logarithm fold change; UP: Significant upregulation; Down: Significant downregulation; None: Non-significant changes. C , D PCA analysis of gene expression in the nine sample groups. Triangle, circle, and square represent normal endometrial tissues of the control women, and ectopic endometrial tissues and eutopic endometrial tissues from patients, respectively. C: PCA of DEGs. D: PCA of non-DEGs Gene expression array analysis for HOX genes. There were three biological repeats in each group. A, B, C and D in each graph show the nine data that were analyzed with SPSS software. P  < 0.05 was considered to be significant A Venn diagram for DEGs ( P  < 0.05). DEGs in the ectopic samples compared to normal samples (black ellipse), eutopic samples in comparison to normal samples (green ellipse), and ectopic samples compared to eutopic samples (red ellipse). Genes that were significantly upregulated are enclosed in black, green and red rectangles according to the above mentioned legend. In the multi-colored rectangles, each color represents upregulation in the relevant analysis in the corresponding agreement. B Heat map for all genes. Right column is logFC of ectopic samples compared to normal samples and the left column is the logFC of eutopic samples in comparison to normal samples. logFC: Logarithm fold change; UP: Significant upregulation; Down: Significant downregulation; None: Non-significant changes. C , D PCA analysis of gene expression in the nine sample groups. Triangle, circle, and square represent normal endometrial tissues of the control women, and ectopic endometrial tissues and eutopic endometrial tissues from patients, respectively. C: PCA of DEGs. D: PCA of non-DEGs For further investigation, we performed PCA on the nine studied samples from the ectopic, eutopic, and control groups by separate analysis of the DEGs and non-DEGs (Fig. 2 C,D). Each point in a PCA graph represented one sample and all samples were properly classified into separate clusters based on the DEGs; however, the non-DEGs could not correctly classify the same samples. Then, we performed hierarchical clustering in the DEGs (Fig. S 4 -A) and non-DEGs (Fig. S 4 -B) in the eutopic, ectopic, and control samples . Hierarchical clustering analysis for the DEGs showed that control tissues were clustered together, and the eutopic and ectopic samples were classified in separate clusters (Supplementary figure S 4 ). We identified three co-expression networks based on pairwise correlation of the gene expression data and visualized them using Cytoscape (Fig.  3 A). Network A was mainly occupied by genes that had significant expression changes in eutopic tissue compared to control tissue (78.6%); however, most genes in networks B and C showed significant expression changes in ectopic tissues. A large number of genes from all of the networks had significant expression changes in both ectopic and eutopic tissues compared to control tissues (Fig.  3 B). Fig. 3 A Co-expression gene network. Blue, yellow, and green colors represent genes whose expressions significantly changed in the ectopic, eutopic and both tissues, respectively. B The number of genes with significant different expressions between sample tissues. All genes: The total number of nodes, Eutopic: Number of nodes with significant changes in gene expressions of the eutopic to control tissues, Ectopic: Number of nodes with significant changes in gene expression of the ectopic to control tissues, All DEGs: Number of nodes with significant changes in gene expressions of the patient to control tissues A Co-expression gene network. Blue, yellow, and green colors represent genes whose expressions significantly changed in the ectopic, eutopic and both tissues, respectively. B The number of genes with significant different expressions between sample tissues. All genes: The total number of nodes, Eutopic: Number of nodes with significant changes in gene expressions of the eutopic to control tissues, Ectopic: Number of nodes with significant changes in gene expression of the ectopic to control tissues, All DEGs: Number of nodes with significant changes in gene expressions of the patient to control tissues The KEGG and Reactome databases were searched to determine the signaling pathways with involvement of networks A, B, and C. The number of genes grouped in networks A, B, and C were involved in “activation of anterior HOX genes in hindbrain development during early embryogenesis”. Also, the number of genes grouped in network A were involved in “signaling pathways regulating pluripotency of stem cells”. Some genes from network B were involved in “transcriptional misregulation in cancer” (Table 1 ). Table 1 Pathways with genes in the constructed networks in the study Gene symbol ISL1, OTX1 HOXB2, HOXB3, HOXB4 HOXB1, HOXD1, HOXD3 SIX1, SIX4 Network name A B C B Pathway name S R R T R: Activation of anterior HOX genes in hindbrain development during early embryogenesis S: Signaling pathways regulating pluripotency of stem cells T: Transcriptional misregulation in cancer Pathways with genes in the constructed networks in the study R: Activation of anterior HOX genes in hindbrain development during early embryogenesis S: Signaling pathways regulating pluripotency of stem cells T: Transcriptional misregulation in cancer DEGs and the genes in networks A, B, and C were mapped to the DAVID database to investigate GO, interference with other diseases, and expressions in other tissues (Tables 2 , 3 , S 2 and S 4 ). Table 2 List of diseases that show interference by the DEGs and genes from the five networks Description Count Percent P -value DEGs in ectopic to control  Developmental 19 44.19 2.14E-11  Metabolic 17 39.53 0.038479  Neurological 12 27.91 0.024028  Cleft lip/Cleft palate 8 18.60 1.21E-05  Bone mineral density 7 16.28 1.64E-04  Vision 7 16.28 0.002492  Parkinson's disease 4 9.30 0.003099  Colobomatous microphthalmia 3 6.98 1.86E-04  Clubfoot 3 6.98 0.001869  Obstructive sleep apnea 3 6.98 0.00642  Hirschsprung's disease 2 4.65 0.034825 DEGs in eutopic to control  Developmental 16 43.24 8.90E-09  Neurological 10 27.03 0.084075  Psychological 9 24.32 0.027439  Cleft lip/Cleft palate 6 16.22 8.20E-04  Parkinson's disease 5 13.51 1.46E-04  Autism 5 13.51 0.004516  Bone mineral density 5 13.51 0.007412  Vision 5 13.51 0.041046 Genes in network A  Metabolic 13 46.43 0.03939  Developmental 12 42.86 1.51E-06  Psychological 11 39.29 1.60E-04  Autism 6 21.43 1.30E-04  Cleft lip/Cleft palate 6 21.43 2.09E-04  Bone mineral density 6 21.43 2.51E-04  Parkinson's disease 5 17.86 4.70E-05  Talipes equinovarus 2 7.14 0.005847  Neurodevelopmental psychiatric disorders 2 7.14 0.007304  SIDS/sudden infant death syndrome 2 7.14 0.038844 Genes in network B  Cleft lip/Cleft palate 3 25 0.012792  Cleft lip with and without cleft palate 2 16.67 0.013342  Bone mineral density 3 25 0.013805  Developmental 4 33.33 0.023658  Sleep apnea, obstructive 2 16.67 0.026536 Genes in network C  Developmental 7 46.67 1.47E-04  Clubfoot 3 20 1.97E-04  Bone mineral density 4 26.67 0.002301  Metabolic 9 60 0.004216 Table 3 Other tissues in which the DEGs and the genes in five networks are expressed Description Count Percent P-value DEGs in ectopic to control  Cerebellum 30 69.77 0.007552  Salivary gland 30 69.77 0.023191  Lymphoma Burkitts Raji 28 65.12 0.043332  Adrenal cortex 13 30.23 0.004518  Craniofacial 4 9.3 1.32E-06  Retina 4 9.3 0.036957 DEGs in eutopic to control  Salivary gland 26 70.27 0.005005  Lymphoma Burkitts Raji 25 67.57 0.005565  Cerebellum 25 67.57 0.008935  BM-CD33 + myeloid 24 64.86 0.010012  Pons 17 45.95 0.028449  Adrenal cortex 10 27.03 0.023199  Craniofacial 4 10.81 7.19E-07  Ovary neoplasia 3 8.11 0.033658 Genes in network A  Cerebellum 19 67.86 0.002102  Salivary gland 19 67.86 0.005265  Lymphoma Burkitts Raji 18 64.29 0.012292  Prostate 17 60.71 0.015867  BM-CD33 + myeloid 17 60.71 0.03038  Pons 14 50 0.007547  Eye_normal 13 46.43 0.032975  Craniofacial 2 7.14 0.012725 Genes in network B  Uterine tumor_disease 7 58.33 0.003695  Embryo_development 8 66.67 0.01666  Cingulate cortex 5 41.67 0.025441  Heart 6 50 0.047369 Genes in network C  Subthalamicnucleus 6 40 0.020087  Skin_normal 7 46.67 0.027182  Embryo_development 10 66.67 0.032412  Germ cell tumor_disease 8 53.33 0.038866 List of diseases that show interference by the DEGs and genes from the five networks Other tissues in which the DEGs and the genes in five networks are expressed Many DEGs in ectopic tissues play a role in developmental, neurological, and metabolic diseases. Numerous DEGs in eutopic tissues are also involved in developmental and neurological disease classes such as Cleft lip, Cleft palate, and Parkinson's disease. As shown in Table 2 , the genes in networks A, B, and C were mostly involved in developmental diseases, while the genes in network A were involved in metabolic diseases. Many DEGs in the eutopic and ectopic tissues were also expressed in different tissues, such as the cerebellar tissue, craniofacial tissue, salivary glands, and the adrenal cortex. Moreover, DEGs in the eutopic tissues were also identified in Burkitt lymphoma Raji and BM-CD33 + Myeloid (Table 3 ). Network A genes were expressed in the cerebellar, craniofacial, salivary gland, lymphoma Burkitts Raji and BM-CD33 + Myeloid tissues. Genes in network B were expressed in uterine tumor tissue and embryo development. Network C genes were identified in embryo development (Table 3 ). On the basis of the GO annotation, most DEGs in eutopic and ectopic tissues were involved in the transcription and development or morphogenesis of various systems. “Anterior/posterior pattern specification”, “embryonic skeletal system morphogenesis”, “dopaminergic neuron differentiation”, and “multicellular organism development” were some of the most important pathways listed for DEGs in ectopic tissues. Eutopic tissue DEGs were also involved in the “multicellular organism development”, “neuron development”, and “dopaminergic neuron differentiation” pathways (Table S 3 , S 4 ). “Anterior/posterior pattern specification”, multicellular organism development and “skeletal system morphogenesis” were listed for all of the networks (A, B, and C). Nervous system-related pathways like dopaminergic neuron differentiation and neuron development were listed for networks A and B (Tables S 3 , S 4 ).

The Ethics Committee of the Clinical and Research Centers for Infertility, Yazd, Iran and the Research Ethics Committee of Royan Institute, Tehran, Iran (IR.ACECR.ROYAN.REC.1396.238) approved this cross-sectional study. All methods were carried out in accordance with the approved guidelines [ 12 ]. Participants signed an informed consent in accordance with the guidelines of the Declaration of Helsinki (2000 revision) and tissue samples were collected after receipt of participants’ consent. Study participants comprised 15 women with stages III or IV endometriosis who underwent laparoscopic surgery. The participants ranged in age from 24 to 38 years. Ectopic endometrial samples were collected from ectopic sites in the abdomen and eutopic endometrium samples were obtained by pipelle sampling. None of patients had evidence of any endometrial hyperplasia, visible endometrial hyperplasia or neoplasia, or any inflammatory diseases. The control group consisted of 15 patients (22 to 36 years old) who underwent diagnostic laparoscopy for secondary infertility. Control group participants had no evidence of endometriosis, polyps, myoma, or reproductive inflammatory diseases and had at least one child conceived by natural pregnancy. Normal endometrial samples were taken via pipelle sampling from these women during diagnostic laparoscopy procedures. Confounding factors for samples are listed in Table S1. All study participants did not have hormone therapy or an intrauterine device for three months before the surgical procedures, and all were in the follicular phase of their menstrual cycles. The obtained endometrial samples were placed immediately in an RNA protection reagent, RNAlater (Ambion, Austin, TX), frozen in liquid nitrogen, and stored at -80 °C. We pooled five tissue samples in the ectopic, eutopic, and normal groups as three biological repeats (each combined five tissue samples served as one biological repeat) and the following analyses were performed on the nine samples. RNA extraction of collected tissues was performed by using an RNeasy Microarray Tissue Mini Kit (Qiagen, cat. no: 73304) according to the manufacturer’s instructions. RNA quality and concentration were measured by using a Nanodrop 2000 spectrophotometer (Thermo Scientific). Then, cDNA synthesis was performed with an RT2 First Strand Kit (Qiagen, cat. no: 330404). The PCR array was conducted for the HOX family genes by using an RT2 Profiler PCR Array Human Homeobox Gene Kit (Qiagen, cat. no: PAHS-083Z) and an RT2 SYBR Green ROX qPCR Master Mix (Qiagen, cat. no: 330502) [ 13 ]. There were primer sets for 84 tests and the housekeeping genes on mentioned HOX PCR array kit, according to the manufacturer’s instructions. Cycling conditions for the StepOnePlus Real-time PCR system [ 14 ] were: an initial denaturation at 95 °C for 10 min followed by 40 cycles at 95 °C for 15 s, and 60 °C for one minute. Results were expressed as values of the cycle threshold [ 15 ] and then normalized to GAPDH, as the housekeeping gene (control gene). The 2 −∆∆Ct method was used to calculate fold changes in mRNA abundance. The statistical significance of difference in expression levels between groups was analyzed using the ANOVA test. P -values < 0.05 were considered to be statistically significant. Principal component analysis (PCA) [ 16 ] was carried out using Minitab 16 statistical software to simplify the large amount of data. PCA was applied on log2 of differentially expressed genes (DEGs) and non-DEGs data separately, to underlying cluster structures of endometrial samples. Hierarchical clustering [ 17 ] was implemented by using the correlation coefficient and complete linkage in Minitab 16 statistical software [ 18 ]. The methods were performed twice with normalized gene expression data from the DEGs and non-DEGs among the sample groups, in three clusters. The gene co-expression networks were constructed based on the normalized gene expression data. MATLAB was used to compute the Pearson correlation between each pair of genes, then the significant correlation pairs ( p  < 0.05) were imported into Cytoscape software (version 3.6.1) [ 19 ] for visualization and the circle algorithms were used [ 20 ]. The co-expression networks of the normal tissues were constructed. Gene function was annotated based on Gene Ontology (GO), the intervention in other diseases, and expression in other tissues by using DAVID ( https://david.ncifcrf.gov ) [ 21 ] to clarify the function of the DEGs and the mechanism of endometriosis.

Endometriosis is an estrogen-dependent gynecological disorder that affects 10%-15% of reproductive age women [ 1 , 2 ]. It is characterized by the presence of endometrial tissues located outside of the uterus. The symptoms include dysmenorrhea, pelvic pain, pelvic masses, and infertility [ 3 – 6 ]. Social relationships, sexuality, and mental health are affected by endometriosis [ 7 ]. There is no definitive cure for all patients who suffer from endometriosis. Therefore, developing advanced efficient therapies in addition to prevention of endometriosis development and recurrence are of utmost importance [ 8 ]. It is essential to demystify the molecular mechanisms that underlie endometriosis in order to develop advanced therapies for this disease. In this regard, HOX genes are remarkable because they regulate endometrial development and its receptivity during implantation [ 9 ]. These genes are essential for endometrial growth, differentiation, and receptivity because they mediate some functions of the sex steroids during each reproductive period [ 10 ]. Dysregulation in HOXD3 [ 10 ], HOXA11 [ 11 ] , and HOXB3 [ 10 ] have been reported in endometriosis patients. However, the changes in the expression levels of other HOX genes or HOX -cofactor genes in patients with endometriosis have not been elucidated and interactions between these genes remain to be determined. Access to expression profiles of these genes would provide the information needed to uncover new genes that could assist with endometriosis diagnosis and treatment. In this study, we report about the expression profiles of HOX genes, their cofactors, and the interactions between these genes.

Our study provides the first insight into the expression pattern of HOX clusters and their cofactors in endometriosis. The results suggest an important role for these genes in the pathophysiology of this disease. Some genes from this platform should be additionally investigated in the future because of their thousand-fold changes in endometriosis tissue biopsies compared to control tissue biopsies.

In this study, we have observed disruptions in expression patterns of HOX genes and their cofactors in eutopic and ectopic samples from endometriosis patients compared to samples derived from women without endometriosis. The data demonstrated that most of the evaluated genes had altered expression patterns in the eutopic and ectopic samples from endometriosis patients compared to the control samples. Genes involved in the HOXA-D clusters showed disturbed patterns of expression in tissues derived from women with endometriosis compared to the control tissues. In line with our findings, several studies have reported an association between HOX genes and endometriosis, as well as cancer [ 10 , 22 , 23 ]. Differential expressions of HOXB2, B3 and B4 [ 10 , 24 ], and the HOXD2 and D3 genes [ 10 ], which we detected in the current study, agreed with previous reports on endometriosis. Altered expressions of HOXB2, B3, B4, B7, and B9, which we observed in the endometriosis samples, correlated with previous reports for these genes in endometriosis and in progression and development of many cancers [ 24 – 26 ]. Here we reported misregulation of HOXC cluster genes in endometriosis. Their involvement with endometriosis, cancer, and the female reproductive tract have been shown previously [ 10 , 27 , 28 ]. Moreover, we reported contrariwise expression patterns for HOXD12 and HOXD13 (downregulated) and HOXD1 and D3 (upregulated) in this study on endometriotic lesions. Upregulated HOXD genes have a correlation with their expression profiles in some carcinomas [ 29 ]. HOXC12 , EN1 , SHOX , TLX1 , VAX1 , VSX1 , PAX3 , PITX3 , DLX1, and LBX1 showed intensive changes as they are upregulated in eutopic tissues of women with endometriosis compared to the control group. Accordingly, we reported that most of these genes were hypomethylated, which correlated with elevated expressions of these genes in endometriosis patients [ 30 ]. The expression level of HOXC12 is increased in breast tumor samples [ 31 ] and TLX1 shows misregulation in cancer [ 32 ]. Therefore, our finding re-enforces the previously reported association between cancer and endometriosis [ 33 , 34 ]. There is PAX3 expression in cytotrophoblast cells and decidua cells during early pregnancy [ 35 ], so its disturbed expression might be involved in implantation defects in endometriosis. EN1 , PITX3, and SHOX2 are involved in neuron formation in eutopic and ectopic tissues [ 36 , 37 ], and the altered expressions in these genes may be related to nerve formation and pain sensation in endometriosis patients [ 38 ]. Most DEG genes in network A showed significantly altered expressions in ectopic tissues. These genes are mainly involved in development and nerve function. In contrast, most genes categorized in networks B and C showed significantly altered expressions in eutopic tissue and are involved in developmental functions (Fig.  3 A). Our analysis indicated that ectopic and eutopic DEG genes are mainly expressed in the salivary gland and Burkett's lymphoma (Table 3 ). We detected pathways related to stem cell pluripotency for DEGs, which was in line with previous reports for this pathway in endometriosis [ 39 , 40 ]. Moreover, we identified a “transcriptional misregulation in cancer” pathway for the studied genes that correlated with a previously proven relationship between endometriosis and cancer [ 41 , 42 ]. HOX gene expressions are regulated by the hormones estradiol and progesterone. Endocrine regulation of HOX gene expressions is important for reproduction [ 27 , 43 ]. Because endometriosis is associated with disturbed hormonal function (estrogen hyperactivity and progesterone resistance) [ 44 ], it is possible that an altered expression profile of HOX in endometriosis tissues is a result of hormonal dysregulation. The results of this study provide preliminary evidence for the involvement of certain pathways (neural and immune system, as well as cancer) in endometriosis because the DEG genes identified by our study are reported to play a role in the above pathways. However, the role of these pathways in endometriosis should be investigated. It remains to be determined how neurogenesis or cancer pathways induce endometriosis or promote the development of this disease. Recently developed endometriosis organoids are suitable pre-clinical models to investigate molecular mechanisms that underlie this disease [ 30 , 45 ]. We suggest knock out the upregulated genes with extensive fold changes like HOXC12 in endometriosis organoids and over-express these genes in endometrial organoids developed from non-endometriosis patients to explore the role of these genes in the development of endometriosis. Endometriosis/endometrium organoids and stromal cell cultures provide a suitable platform to investigate the role of progesterone and estrogen in the disturbed expression profiles of HOX genes.

Additional file 1. Additional file 1.