Genome-wide expressions in autologous eutopic and ectopic endometrium of fertile women with endometriosis

Background In order to obtain a lead of the pathophysiology of endometriosis, genome-wide expressional analyses of eutopic and ectopic endometrium have earlier been reported, however, the effects of stages of severity and phases of menstrual cycle on expressional profiles have not been examined. The effect of genetic heterogeneity and fertility history on transcriptional activity was also not considered. In the present study, a genome-wide expression analysis of autologous, paired eutopic and ectopic endometrial samples obtained from fertile women (n = 18) suffering from moderate (stage 3; n = 8) or severe (stage 4; n = 10) ovarian endometriosis during proliferative (n = 13) and secretory (n = 5) phases of menstrual cycle was performed. Methods Individual pure RNA samples were subjected to Agilent’s Whole Human Genome 44K microarray experiments. Microarray data were validated (P < 0.01) by estimating transcript copy numbers by performing real time RT-PCR of seven (7) arbitrarily selected genes in all samples. The data obtained were subjected to differential expression (DE) and differential co-expression (DC) analyses followed by networks and enrichment analysis, and gene set enrichment analysis (GSEA). The reproducibility of prediction based on GSEA implementation of DC results was assessed by examining the relative expressions of twenty eight (28) selected genes in RNA samples obtained from fresh pool of eutopic and ectopic samples from confirmed ovarian endometriosis patients with stages 3 and 4 (n = 4/each) during proliferative and secretory (n = 4/each) phases. Results Higher clustering effect of pairing (cluster distance, cd = 0.1) in samples from same individuals on expressional arrays among eutopic and ectopic samples was observed as compared to that of clinical stages of severity (cd = 0.5) and phases of menstrual cycle (cd = 0.6). Post hoc analysis revealed anomaly in the expressional profiles of several genes associated with immunological, neuracrine and endocrine functions and gynecological cancers however with no overt oncogenic potential in endometriotic tissue. Dys-regulation of three (CLOCK, ESR1, and MYC) major transcription factors appeared to be significant causative factors in the pathogenesis of ovarian endometriosis. A novel cohort of twenty-eight (28) genes representing potential marker for ovarian endometriosis in fertile women was discovered. Conclusions Dysfunctional expression of immuno-neuro-endocrine behaviour in endometrium appeared critical to endometriosis. Although no overt oncogenic potential was evident, several genes associated with gynecological cancers were observed to be high in the expressional profiles in endometriotic tissue.


Background
Endometriosis is a complex disorder involving pathogenesis and clinical presentation of ectopically implanted endometrium [1]. It is generally assumed that elucidation of molecular expressional specificities of eutopic and ectopic endometrium may provide leads towards a better understanding of the pathophysiology of the disorder [2]. To this end, several studies exploring the differential expression of genes between autologous eutopic and ectopic endometrium from patients with endometriosis have been reported, however, with no specific comparison for stages of severity, fertility history and phases of menstrual cycle [3][4][5][6][7], except a recent report [8]. More over, it is notable that two types of endometriosis, namely ovarian endometriosis and peritoneal endometriosis reportedly show differential characteristics [4,9]. Furthermore, there is evidence to support the idea that deep infiltrating endometriosis also show differential pathophysiology as compared to ovarian and peritoneal endometriosis [10,11]. In the present study, we examined a genome-wide large-scale transcript survey of autologous, paired eutopic and ectopic endometrial samples obtained from fertile women suffering from moderate to severe ovarian endometriosis, and excluded cases of peritoneal endometriosis and deep infiltrating endometriosis. We assumed that the present model of subject selection would reduce the impact of biological noise derived from genetic and pathogenetic heterogeneity and subfertilityassociated variability on the transcriptional activity in the target tissue. We report here for the first time that clustering effect of expressional arrays among eutopic and ectopic samples was higher for genetic homogeneity (i.e. pairing of eutopic and ectopic samples from same individuals) than that of clinical stages of severity and phases of menstrual cycle. Based on the present transcriptomics data, we have also hypothesized that dysfunctional immuno-neuro-endocrine behaviour in endometrium was associated with the pathogenesis of endometriosis. Additionally, we did not observe an overt oncogenic potential in the expressional profiles in endometriotic tissue, however, several genes associated with gynecological cancers were highly expressed in the eutopic and ectopic endometrium. Finally, a novel cohort of 28 genes was identified, the expression of which carry potential marker value for endometriosis in fertile women. A flow diagram of the experimental design is shown in Figure 1.

Subjects and tissue samples
The present study was approved by the Ethics Committee on the Use of Human Subjects, All India Institute of Medical Sciences (AIIMS), New Delhi. The patients enrolled in the Department of Obstetrics and Gynecology -AIIMS and showing evidence of endometriotic lesions, adhesions and endometriotic cyst were selected to participate in the present study. All the patients were reportedly fertile and referred from the Pain Clinics, and had voluntarily agreed to donate their samples after understanding the purpose of the proposed study. Signed informed consent was obtained from each participant of this study. As shown in Figure 1, twentysix (26) normally cycling and proven fertile women (age: 24-45 y) with history of pregnancy and with at least one living biological offspring, and body mass indices within normal ranges (20-22 k/m 2 ) having ovarian endometriosis were selected for the present study. Confirmation of ovarian endometriosis and exclusion of other types of endometriosis was achieved from reports of pelvic imaging based on ultrasound, MRI and/or diagnostic laparoscopy as described elsewhere [8]. Severity stages 3 and 4 of the disease condition were defined at the time of surgical laparoscopy [8] according to rASRM protocol [12]. Selected subjects (n = 18; shown as 'E' in Additional file 1: Table S1) contributed their eutopic (shown as ' A' in Figure 2) and ectopic (shown as 'B' in Figure 2) samples during proliferative (days 9-14) phase (n = 17) and secretory (days 17-24) phase (n = 8) of menstrual cycle as described elsewhere [8]. Additional paired samples collected from different group of subjects (n = 8; shown as 'Ep' in Additional file 1: Table S1) with confirmed ovarian endometriosis as described above and with classified menstrual (proliferative: n = 4; secretory: n = 4) phases and severity stages 3 (n = 4) and 4 (n = 4) were employed for validating the prediction as described below. A small piece from each specimen was processed for chemical fixation in neutral buffered formaldehyde (4%, w/v) for subsequent confirmation of phase of cycle, state of pathology and cell types from eutopic and ectopic samples, and the residual portions were transported on ice to the laboratory within 10 minutes of collection for further processing for RNA extraction.

Experimental procedure
The methodological details of RNA extraction followed by the estimation of its yield and purity using standard electrophoretic and spectrophometric protocols and its RIN score using the Agilent 2100 Bioanalyzer, RNA 6000 Nano LabChip kit and Agilent 2100 Expert Software (Agilent Technologies, Santa Clara, CA, USA) have been given elsewhere [8,13]. Individual RNA samples from eutopic and ectopic tissue samples (n = 18) from confirmed stages 3 (n = 8) and 4 (n = 10) collected during proliferative (n = 13) and secretory (n = 5) phases and having RIN scores >8.0 were subjected to whole transcriptome array experiment using the Agilent Whole Human Genome 60-mer 4X44K microarray according to the manufacturer's recommendations. Thus, seven (7) samples could not be used either for insufficient RNA yield or RIN scores (see Additional file 1: Table S1 for the subject details of the selected samples). Hybridized arrays were scanned with Agilent's G2505B microarray scanner system and the raw data were imported into GeneSpring 11.5.0 software (Agilent Technologies, Santa Clara, CA, USA) for further analysis. Pearson's correlation coefficients done to assess the reliability of data obtained from two separate hybridization runs for same RNA preparation for four (4) eutopic and ectopic samples confirmed the reproducibility assurance (P < 0.01) among hybridizations. Analysis of the data retrieved  from separate chips with the same RNA samples yielded QC statistics highly concordant with that of the manufacturer, and it revealed more than 95% confidence level.

Data analysis
Unsupervised and supervised hierarchical clustering analysis (HCA), and non-hierarchical K-mean cluster analysis of expression arrays were performed with the help of GeneSpring software 11.5.0. Analysis of variance followed by pair-wise differential (>3-fold at P < 0.01) expression (DE) for specific genes between eutopic and ectopic samples, as well as, between proliferative and secretory phases, and between clinical stage 3 and stage 4 of severity for eutopic endometrium, and for ectopic  endometrium, respectively were done using multiple comparison tests as described elsewhere [14].

Post-hoc analysis
Networks and enrichment analysis were done using gene lists obtained from the above analyses and based on a priori setting of a cut-off threshold (pFDR(p) = 0.05) with the help of the GeneSpring11.5.0 software and Metacore platform (GeneGo, St. Joseph, MI, USA). The K-mean clusters were further used for differential co-expression (DC) analyses and analyzed in terms of Gene Ontology (GO) enriched categories using GeneSpring11.5.0 software. Gene Set Enrichment Analysis (GSEA) version 3.7 was applied to each of the K-mean clusters independently to examine at FDR ≤ 0.25 for not less than 10 genes for a set with a maximum of 1000 permutation whether preannotated BROAD gene sets [15]: C1 (cytogenetic sets), C2 (functional sets), C3 (regulatory sets), C4 (cancer neighborhood sets), and C5 (gene ontology sets) could identify any interesting information in the DC sets [16].

Quantification of candidate gene expression by real time RT-PCR
In order to validate the microarray data, relative expression of arbitrarily chosen seven (7) selected genes (ATX, DDHD1, DYNLT1, FTH1, LAMR1, MIER2, and WDR87) in eutopic and ectopic samples collected from all patients were performed using Taqman multiplexing technology on iCycler iQTm real time RT-PCR detection system (BioRad, Hercules, CA, USA). GAPDH was selected as an endogenous control based on its observed expressional consistency in arrays on data analysis. Primers and probes were designed on Beacon Designer software7 (Labware Scientific Inc., Milipitas, CA, USA) and obtained from Qiagen (Cologne, Germany) (see Additional file 2: Table S2 for the details). The ratio of estimated efficiency of the primers for the selected genes and GAPDH was~1.0. An optimized kit (QuantiTect multiplex PCR kit, Qiagen, Cologne, Germany) was used to synthesize cDNA from respective RNA (5 μg) samples. Relative expression ratios between groups were calculated by using 2 -ΔΔCt method [17]. Quantification of copy numbers for target transcripts in complex RNA samples was obtained as described elsewhere [18]. Comparison between fold change data obtained from real time RT-PCR and microarray image analysis for selected seven (7) genes revealed a high degree of concordance and pattern similarity in expression profile. Concordance correlation test between real time RT-PCR based quantitative data and microarray data for the seven (7) genes showed a high degree of correlation (P < 0.01) [8].
In order to test the reproducibility of prediction derived from analysis of microarray results, the relative expressions of twenty eight (28) selected genes in individual RNA samples obtained from eight (8) additional subjects giving paired eutopic and ectopic samples with confirmed endometriosis stages 3 (n = 4) and 4 (n = 4) during proliferative (n = 4) and secretory (n = 4) were examined using real time RT-PCR technology. The subject details are shown in Additional file 1: Table S1. The genes selected based on GSEA implementation of DC results were employed to test the predictability function of the expression of those genes. The details of RNA methodologies are given above. GAPDH was selected as an endogenous control based on its observed expressional consistency in arrays on data analysis. All primers were designed on the Beacon Designer software7.0 (Labware Scientific Inc., Milipitas, CA, USA) based on SYBR green chemistry and obtained from Qiagen (Cologne, Germany). QuantiTect Reverse Transcription kit for cDNA synthesis and QuantiFast SYBR green PCR kit for PCR amplification from Qiagen (Cologne, Germany) were used according to the protocol given by the manufacturer. The estimates of relative expression ratios between groups and copy numbers for target transcripts in complex RNA samples were obtained as described above.

Results
The data sets are available at NCBI-GEO website [19].
A distribution histogram of the number of probes and genes for different ranges of expression in autologous, paired eutopic and ectopic samples obtained from eighteen (18) fertile women with confirmed ovarian endometriosis is shown in Figure 2A. Total numbers and per cent estimates of probes/genes expressed in eutopic and ectopic samples in optimized scale are shown in Table 1. On average,~75% and~50% of expressed genes showed marked signal in eutopic and ectopic samples, respectively.   Table S3. Note that the areas in the Venn distribution analysis are not drawn to scale. 0.2) than that of phases of cycle (cd: 0.3), but not than that of the clinical stages of severity (cd: 0.1).

Differential expression (DE)
Additional file 3: Table S3 gives the list of the genes along with their differential expression (DE) patterns under different categories based on expressional arrays in autologous, paired eutopic and ectopic samples obtained from 18 fertile women with ovarian endometriosis. Figure 3 shows the number of genes with DE in different categories of comparison and the lists of common genes in it. Table 2 highlights the enriched categories of pathways for the common genes from above-mentioned DE analysis between eutopic and ectopic endometrium. It appeared that different signaling pathways associated with immune response, several neuronal processes, and ERBB family signaling pathways were commonly selected. A summary of DE analysis of the non-common genes showing differential display under different categories and their enrichment analysis are shown in Table 3. Collectively, it appeared that informational flow for a wide array of pathways involving cellular signaling, apoptosis and survival, cytoskeleton remodeling, chemotaxis, cell adhesion, immune response and several neurophysiological processes were affected.
K-mean clusters and differential co-expression (DC) As shown in Figure 4, K-mean cluster analysis identified four clusters of expression patterns and profiles based on normalized hybridization signals for all expressed genes in all samples. The genes in cluster 1 (K1) did not show any specific expression pattern, while other three clusters showed overt patterns for menstrual cycle phases and severity stages. A large number of genes belonging to cluster 2 (K2) showed over-expression in severity stage 4 secretory phase endometrium (Fig. 4B). The co-expressed genes in cluster 3 (K3) and cluster 4 (K4) showed very similar patterns with an overall higher expression in stage 3 as compared to stage 4 endometrium samples irrespective of cycle phases. Table 4 shows the pathways-based enrichment analysis of groups of genes in four (4) K-mean clusters revealing differential co-expression (DC) profiles between paired eutopic and ectopic endometrial tissues. It essentially substantiated the observation obtained from DE analysis that transcriptomic signals related to cell cycle, signal transduction, cytoskeleton remodeling, apoptosis and survival, chemotaxis, cell adhesion, and immune response were affected in the pathogenesis process of endometriosis.
Gene-set enrichment analysis (GSEA) Table 5 provides a summary of the results of GSEA implementation on co-expressed genes with differential display (DC) in the four K-mean clusters. In K4, one (1) cytoband i.e. C1 set and two (2) gene ontology i.e. C5 sets were selected. More over, three (3) DC gene setsone each in K1, K2 and K4, respectivelywere selected under BROAD regulatory gene motif sets, C3. It is notable that two (2) selected regulatory motif sets belonging to K1 and K2 were significantly (p < 0.0001) associated with ectopic sample as evident from their negative normalized enrichment scores (NES). Further, four (4) DC gene setstwo (2) each in K1 and K4, respectively -  Figure 2, b A, c B and d C.     were selected under BROAD cancer gene neighborhood sets, C4. Table 5 also shows the major gene families selected in GSEA and names of the genes showing differential display in the comparison between eutopic and ectopic endometrium.
Expressional cohort of marker genes Table 6 provides the list of selected twenty eight (28) genes that appeared significant from combined analysis of GSEA-selected gene sets followed by DC analysis and from DE analysis of microarray data of 18 paired samples. Table 6 also shows that the validity of the prediction value of the expressional cohort based on quantitative analysis in a different set of autologous, paired eutopic and ectopic samples obtained from a separate group of 8 subjects was markedly high.

Discussion
The awareness that whole genome expression array analysis may yield high dimension knowledge towards deciphering patho-etiology of complex diseases [20,21] has prompted several groups of investigators to employ this approach to examine the transcriptomics basis of endometriosis using eutopic and ectopic samples [3][4][5][6][7][8][9]. Although significant and interesting observations have emerged from these reports, these studies did not include the possible impact of one or more of the factors like the demographic characteristics, position of endometriosis, fertility history, severity stages and phases of menstrual cycle influencing the genomic expression in eutopic and ectopic tissues [2,8,22]. In the present study, we have examined the whole genome transcriptomics of autologous, paired eutopic and ectopic samples obtained from fertile Indian women with ovarian endometriosis of known clinical severity and phases of menstrual cycle but with no history of previous treatment for endometriosis at the time of tissue collection. We analyzed the expression profiles to delineate the impact of stages of severity and phases of cycle in eutopic and ectopic samples. We observed that clustering effect on expression arrays was maximum in paired samples, followed by stages of clinical severity and positional cue. The phase of menstrual cycle exhibited minimal clustering effect on expressional profiles in the experimental samples. Generally, we observed that eutopic tissue yielded a normal frequency distribution histogram of gene expressions for different levels of expression and that an overall higher numbers of genes in eutopic endometrium expressed higher transcriptomic signals as compared to ectopic samples; ectopic samples yielded a truncated frequency distribution histogram. Furthermore, higher numbers of genes bearing expression levels at the high and at the low ends of the frequency distribution were observed in ectopic tissues as compared to eutopic tissues. We believe that implementation of appropriate computational models based on Shannon's noise-signal entities and probability of size of loss of signals may yield in the future new leads about the global genomic expression pattern in the ectopic tissue [23,24]. It is notable in this regard that: (i) a large number (~0.7 K) of genes were silenced in the ectopic tissue at stage 4 condition as compared to stage 3, and (ii) expressional clustering cohesion was very high (cd: 0.07) between the eutopic and ectopic endometrium in stage 4 disease condition. Taken together, it is suggestive of high degree of pathognomonicity in stage 4 eutopic endometrium [8,25]. Major highlights in the previous studies on large scale expressional array analysis were to explore the genespecific DE in paired analysis between eutopic and ectopic endometrium with an assumption that a 2-fold change at P < 0.05 between two groups of tissue samples was sufficiently significant for further analysis. As pointed out elsewhere, this may give rise to different sets of biases and inadequacies in interpretation and discovery [26,27]. To circumvent these acknowledged insufficiencies, we have employed a 3-fold change at P < 0.01 as the pre-set filter for DE of individual genes between groups followed by pathway networks based enrichment analysis, and for DC analysis of K-mean based expressional clusters followed by gene set enrichment analysis (GSEA) model [16] to interpret the present transcriptomics data.
Post-hoc analysis of expressional signals in eutopic and ectopic tissues under different sets of categorical comparison revealed that several signaling pathways related to immune response were commonly affected in eutopic and ectopic endometrium. The results from the present study support the observation made by Zhao et al. based on GSEA of archival transcriptomics data sets of endometriosis that the main canonical pathways putatively involved in the process of endometriosis were related to that of immune and inflammatory diseases [27]. Additionally, we hypothesize from the results of the present study that functional connectivity between overexpression of CLOCK and inflammatory disorder, as well as, between over-expression of genes associated with lipid metabolism and inflammation at the local tissue level are operative in the pathogenesis of endometriosis [28][29][30].
We also observed that neuronal processes involving nNOS signaling pathways and GABA synthesis and metabolism were commonly expressed in both tissue types and a large number of genes involving several signaling pathways (corticoliberin, opioid receptors, serotonin receptors, melatonin, dopamine receptor, neuronal cell adhesion, NGF) associated with neurophysiological processes were up-regulated in stage 3 ectopic endometrium. Earlier the possible involvement of neuroendocrine processes in the pathogenesis of endometriosis has been implicated [31][32][33][34]. Pathways and networks based enrichment analysis of DE of individual genes and DC of gene cohorts in four clusters in the present study revealed that expression of genes in pathways directly and indirectly associated with cell apoptosis and survival, cytoskeleton remodeling, chemotaxis and cell adhesion were differentially affected in eutopic and ectopic samples. Involvement of these pathways in endometriosis has earlier been reported by several groups based on different experimental models [35][36][37][38].
The pathogenesis of endometriosis has also been associated with excessive production of estrogens by up-regulated expression of aromatase and 17β-HSD type 1, and suppression of 17β-HSD types 2 and 4, collectively resulting in an increased ratio of estradiol-17β to estrone in ectopic tissues [39][40][41][42]. It has however been reported by others that mRNA and protein expressions of aromatase were minimal in ectopic tissues [42,43]. Our transcriptomics data also failed to identify any overt change in the expression of genes for aromatase (CYP19A1) and 17β-HSD (HSD17B1-B17) in Table 4 Estimates and description of top scored enriched categories of differentially regulated (eutopic-to-ectopic) co-expressed (DC) genes Cluster identity (Total number of genes) [Number of genes with differential display: Up-regulated/Down-regulated] Enriched category (p-value) Gene symbols of major candidate genes (Vector of differential)  endometriosis. However, our observation that genes (NR5A1, STAR) for steroidogenic factor (SF)1 and steroidogenic acute regulatory protein (StAR), which are known to be significant regulators of steroidogenesis [44,45] were highly expressed in the ectopic endometrium substantiates previous reports [46][47][48]. An over-expression of ERBB3 in proliferative phase eutopic endometrium and secretory phase ectopic endometrium was seen in the present study; it has been associated with ligand-independent activation of estrogen receptors (ESR1 and ESR2) in target tissues [49]. However, ERBB3 was up-regulated in eutopic tissues as compared to ectopic tissues in pooled analysis. The activation of ESRs and relative down-regulation of progesterone receptor (PGR) in the secretory phase ectopic endometrium is suggestive of relative suppression of progesterone action in the ectopic endometrium. Indeed, the phenomenon of progesterone resistance in the ectopic endometrium has earlier been documented by other groups [8,30,50].
neuro-endocrine responses may bear vulnerability to give rise to endometriotic lesion if deposited ectopically. Subsequently, the ectopically placed endometrial tissue with positional input and under fluctuating levels of sex steroid hormones in vulnerable subjects may develop differential expression repertoire related to cell survival, adhesion, migration and growth resulting in endometriotic lesion [51,52]. More over, it appears that a pathwaysnetwork of several transcription factors including CLOCK-ESR1-MYC may be involved at the transcriptomic level towards pathoetiology of ovarian endometriosis  Based on transcript numbers from RT-PCR analysis of RNA samples from a different set of autologous paired eutopic and ectopic tissues (n = 8). *Genes selected in gene families identified in GSEA gene sets (see Table 5).
( Fig. 5). Further studies are warranted to test this hypothesis. Endometriosis, by definition, is a benign disease, however, there are a few reports indicating risk of malignant transformation in endometriosis [53][54][55]. In the present study, we observed a general suppression in the expression of genes associated with cell cycle and DNA damage repair in both eutopic and ectopic endometrium in fertile women with endometriosis. Interestingly, there is a recent report indicating a better survival rate for women with endometriosis for all malignancies combined, and specifically for ovarian and breast cancer, while it was poorer in malignant melanoma [56]. While the present results revealing the lack of overt oncogenic potential in endometriotic tissue concur with some of the earlier reports [6][7][8], genes (CHEK1, ERBB family, laminin gamma and Ki-67) associated with gynecological cancers [57][58][59][60] were highly expressed in autologous, paired eutopic and ectopic tissues. Thus, the possibility of inducement of oncogenic transformation through critical phase transition [61] in the course of endometriosis disease progression cannot be ruled out, especially in the high risk population [62,63].
Finally, we have identified for the first time a cohort of twenty-eight (28) genes with high degree of predictability index for ovarian endometriosis in fertile women. We believe this cohort of genes can be used for further study to discover patho-physiology of ovarian endometriosis. Figure 5 Knowledge-based construction of the pathways-network of transcription factors putatively associated with pathogenesis of endometriosis. The transcription factors were identified from GSEA implementation on co-expressed genes. It is notable that CLOCK, ESR1, and MYC (shown inside blue dotted rectangle) are differentially co-expressed in endometriosis as shown is Table 5.