Genome-wide profiling of miRNA expression patterns in tubal endometriosis

MicroRNA (miRNA) expression profiles in tubal endometriosis (EM) are still poorly understood. In this study, we analyzed the differential expression of miRNAs and the related gene networks and signaling pathways in tubal EM. Four tubal epithelium samples from tubal EM patients and five normal tubal epithelium samples from uterine leiomyoma patients were collected for miRNA microarray. Bioinformatics analyses, including Ingenuity Pathway Analysis (IPA), Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were performed. Quantitative real-time polymerase chain reaction (qRT-PCR) validation of five miRNAs was performed in six tubal epithelium samples from tubal EM and six from control. A total of 17 significantly differentially expressed miRNAs and 4343 potential miRNA-target genes involved in tubal EM were identified (fold change >1.5 and FDR-adjusted P value <0.05). IPA indicated connections between miRNAs, target genes and other gynecological diseases like endometrial carcinoma. GO and KEGG analysis revealed that most of the identified genes were involved in the mTOR signaling pathway, SNARE interactions in vesicular transport and endocytosis. We constructed an miRNA-gene-disease network using target gene prediction. Functional analysis showed that the mTOR pathway was connected closely to tubal EM. Our results demonstrate for the first time the differentially expressed miRNAs and the related signal pathways involved in the pathogenesis of tubal EM which contribute to elucidating the pathogenic mechanism of tubal EM-related infertility. Reproduction (2019) 157 525–534


Introduction
Endometriosis (EM) is the presence of functional endometrial glands and stroma outside the uterine cavity. EM is seen in 5-10% of the general female population, with most affected individuals being of the reproductive age (Woodward et al. 2001). Patients may be asymptomatic or may present with chronic pelvic pain, infertility and dyspareunia (Woodward et al. 2001).
Tubal EM refers to functional endometrial glands and stromal implant involving the fallopian tube and is a part of pelvic EM. Only a few studies have focused on tubal EM. The incidence of tubal EM varies widely and accounts for about 0.29-4% of pelvic EM. (Xue & Zhang 2017). In an earlier study, EM lesions were found in the fallopian tubes of 60% (21/35) EM patients; therefore, tubal involvement in moderate and severe EM is common (Xia et al. 2018). Moreover, EM with tubal EM resulted in lower ciliary beat frequency, lower ciliated cell percentage, weaker muscular contractile activity and lower contraction frequency (Xia et al. 2018). Thus, tubal EM may cause infertility in women by affecting the function of the fallopian tube (Xia et al. 2018).
miRNAs are short non-coding RNA molecules (22-25 nucleotides long), which are transcribed in the nucleus as primary miRNA and undergo a series of maturation steps utilizing the sequential action of the endonucleases, Drosha and Dicer, to attain functional competence (Taganov et al. 2007). miRNAs have diverse functions, including morphogenesis, tissue maintenance, cell growth, differentiation, apoptosis and metabolism (Ouellet et al. 2006, Song & Tuan 2006. Aberrant miRNA expression is associated with a number of human diseases (Singh et al. 2008), including benign gynecological conditions, gynecological malignancies and fertility disorders of the human female reproductive tract. Furthermore, miRNAs can impose epigenetic effects on the pathogenesis of EM (Zondervan et al. 2018).
Previous studies have identified the miRNA expression profiles of eutopic and ectopic endometrium from women with EM (Ohlsson Teague et al. 2009) and differentially expressed miRNAs between women with EM and those without EM (Jia et al. 2013, Okamoto et al. 2015. By performing bioinformatics analyses on miRNA and mRNA microarray data, novel miRNAassociated molecular pathways that are likely to play an important role in the pathophysiology of this disorder were identified (Hull et al. 2008, Ohlsson Teague et al. 2009). Additionally, animal models of EM and microarray methods have identified numerous individual EM-associated transcripts that encode the proteins involved in broad signaling pathways. These pathways mediate inflammation, tissue remodeling, apoptosis, cellular proliferation, angiogenesis and wound healing in endometriotic tissues (Eyster et al. 2007, Flores et al. 2007, Hever et al. 2007, Borghese et al. 2008, Hull et al. 2008. Nevertheless, little is known about the expression profiles and functions of miRNAs in tubal EM.
In this study, the miRNA expression profiles in tubal EM were examined by microarray analysis. Differentially expressed miRNAs between tubal EM and controls were identified, followed by the identification of target genes using an integrated method. These differentially expressed miRNAs may play a potential role in the pathogenesis and development of tubal EM.

Clinical tissue specimens
Seven premenopausal patients with tubal EM (case group) and seven premenopausal patients suffering from uterine leiomyoma (control group) were recruited. Tubal EM and normal fallopian tubes in the control group were confirmed by histopathological examinations and were reviewed by at least two experienced pathologists. This study was approved by the Institutional Ethics Committee of the International Peace Maternity and Child Health Hospital in Shanghai, China (No. GKLW 2015-34) and written informed consents was obtained from all participants. From June 2016 to March 2017, tubal epithelium specimens were obtained from surgically resected fallopian tubes. Specimens were immediately preserved in liquid nitrogen and were transported from operation room to the lab for experiments within 10 min. Patient demographics and clinical pathological data were also collected.

RNA extraction and quality control
Total RNA was extracted and purified using mirVana™ miRNA Isolation Kit without phenol (Cat # AM1561, Ambion), following the manufacturer's instructions and checked for RNA integrity by RNA integrity number (RIN) using an Agilent Bioanalyzer 2100 (Agilent Technologies). The initial sample for chip experiment using total RNA was examined using a NanoDrop nd-2000 spectrophotometer and an Agilent Bioanalyzer 2100 (Agilent Technologies).

MiRNA microarray
miRNA in total RNA was labeled by miRNA Complete Labeling and Hyb Kit (Cat # 5190-0456, Agilent Technologies) following the manufacturer's instructions. Each slide was hybridized with 100 ng Cy3-labeled RNA in a hybridization oven (Cat # G2545A, Agilent Technologies) at 55°C and 0.06 g for 20 h. The Agilent miRNA Complete Labeling and Hyb Kit provides all the needed components for sample labeling and hybridization preparation. These include calf intestinal alkaline phosphatase (CIP), 10× calf intestinal phosphatase buffer, T4 RNA ligase, 10× T4 RNA ligase buffer, DMSO, nuclease-free water, cyanine 3-pCp, 10× GE-blocking agent and 2X Hi-RPM hybridization buffer. After hybridization, slides were washed in staining dishes (Cat # 121, Thermos Shandon, Waltham, MA, USA) with Gene-Expression Wash Buffer Kit (Cat #5188-5327, Agilent Technologies).
The slides were then scanned using an Agilent Microarray Scanner (Cat # G2565CA, Agilent Technologies) and the data were extracted using Agilent Feature Extraction Software (Agilent Technologies) with default settings.

Bio-informatic analyses to predict miRNA-target genes and network
To further investigate the functional roles of miRNAs, putative targets of miRNAs were predicted by TargetScan (http:// www.targetscan.org/), and for more accuracy, the Ingenuity Pathway Analysis software (IPA software, Qiagen) was also used for the prediction. Besides, IPA identified the miRNAtarget gene-disease networks, which connected tubal EM with other similar diseases. Genes were then subjected to enrichment analysis of Gene Ontology (GO) functions, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and DAVID (https://david.ncifcrf.gov). GO analysis was performed to explore the functional roles of target genes in terms of biological processes (BP) and molecular functions (MFs). Pathway analysis is a functional analysis for mapping genes to KEGG. The FDR-adjusted P value denotes the significance of the pathway correlated to the conditions. IPA showed how miRNAs participate in the mTOR signaling pathway.

Quantitative real-time PCR
The expression level of five miRNAs was measured by qRT-PCR in six paired tubal epitheliums (six tubal EM and six control women). Total RNA was isolated using TRIzol (Invitrogen), and cDNA was prepared using the bulge-loop™ miRNA qRT-PCR Starter kit. miRNA qPCR analysis was conducted to determine the miRNA expression using the SYBR Green Mix system (RiboBio, Guangzhou, China) in 384-well plates. The amplification thermal cycling conditions were an initial 95°C for 10 min followed by 40 cycles of 95°C for 2 s, 60°C for 20 s and 70°C for 10 s. After PCR, a dissociation curve was constructed at 70°C for 5 s for detection of PCR product specificity. The bulge-loop™ miRNA qRT-PCR Primer Sets (one RT primer and a pair of qPCR primers for each set) specific for miR-21-3p, miR-223-3p, miR-3663-3p, miR-450a-5p and miR-503-5p were designed by RiboBio. The relative expression https://rep.bioscientifica.com Reproduction (2019) 157 525-534 level of miRNAs was normalized to U6 RNA and calculated using the comparative Ct method (2 −∆∆CT ).

Statistical analysis
Raw data were normalized using the Quantile algorithm, included in the R package AgiMicroRna (López-Romero, P. BMC Genomics. 2011). The obtained array data were normalized to external and internal controls to remove systemrelated variations. The statistical significance of the difference was estimated using the Student 's t-test, fold change >1.5 and FDR-adjusted P values <0.05 were considered statistically significant. RT-qPCR data were analyzed using GraphPad Prism 7 (GraphPad Software Inc.).

Participant demographics and clinical pathological data
After initial sample examinations, nine qualified samples (four patients and five controls) were identified for subsequent miRNA chip experiments. The participants' demographics and clinical pathological data are summarized in Table 1. N302837-BL, N129922-BL, N305896-BL and N308632-BL were included in the case group and N305643-DZ, N305044-DZ, N306224-DZ, N303534-DZ and N306515-DZ were included in the control group. The participants varied in age from 42 to 47 years and were in secretory phase of their menstrual cycles, as estimated by the number of days since their last menstrual cycle and the endometrial histology. The P value of unpaired Student's t-test for age was 0.53 indicating that there was no significant difference in age between the case group and the control group.

Normalizing the raw data
Normalization is an essential step in microarray data preprocessing. It aims to remove non-biologic effects resulting from the experimental process so that biologic effects can be accurately identified (Qin & Satagopan 2009). After homogenization and logarithmic transformation of the original data from single fluorescence chip, boxplots were used to evaluate the data to visually compare overall signal distribution between multiple samples (Fig. 1A). The overall distribution and dispersion of the hybridization signal value data of each chip were estimated using boxplots. The boxplots showed that the distributions of miRNAs for the nine samples were similar after normalization (Fig. 1A).

Identification of differential miRNA profiles
A total of 2532 miRNAs were scanned during microarray analysis. After normalization of the raw data, fold change and Student's t-test were used to screen the differentially expressed genes; the selection criteria were as follows: fold change >1.5 and FDRadjusted P value <0.05; finally, 17 miRNAs were found significantly differentially expressed including 4 upregulated and 13 downregulated miRNAs in tubal EM (shown in Supplementary Table 1, see section on supplementary data given at the end of this article). To be more intuitive, differentially expressed miRNAs with statistical significance between the two groups were identified using hierarchical clustering (Fig. 1B).

Validation of differentially expressed miRNAs with qPCR
To confirm the miRNA-microarray results, three upregulated miRNAs and two downregulated miRNAs with fold change >2 were identified with qRT-PCR.
The results of qRT-PCR (Fig. 1C, D, E, F and G) showed that the expression of miRNAs was consistent with the miRNA-microarray results.

IPA of miRNAs
Based on an analysis of predicted gene targets using the IPA software, one upregulated miRNA hsa-miR-223-3p and two downregulated miRNAs hsa-miR-450a -5p and hsa-miR-503-5p were predicted to top 70 target genes involved in the pathogenesis of EM. These miRNAs were selected because their expression profiling showed a fold change >2 which was confirmed by qRT-PCR. Besides, IPA showed, that among the  Fig. 2A). Among the miR-450a-5p target genes, 22 were associated with endometrioid carcinoma, 22 with endometrial carcinoma, 24 with female genital tract adenocarcinoma, 27 with female genital tract cancer and 27 with tumorigenesis of the reproductive tract (Fig. 2B). Among the miR-503-5p target genes, 41 were associated with endometrioid carcinoma, 43 with endometrial carcinoma, 46 with female genital tract adenocarcinoma, 52 with tumorigenesis of the reproductive tract and 52 were associated with female genital tract cancer (Fig. 2B). More details of the miRNA-gene disease network are shown in Supplementary Table 2.

Integrating miRNAs to functional pathways
TargetScan predicted 4343 target genes on which we performed functional annotation including GO and KEGG analyses. The results of the GO analysis of these genes according to the values in the enrichment score under BP and MF are shown in Fig. 3A. In BP, the top two processes were regulation of nuclear-transcribed mRNA poly(A) tail shortening and positive regulation of nuclear-transcribed mRNA poly(A) tail shortening. In MF, there were type-2 fibroblast growth factor receptor binding and copper ion transmembrane transporter activity. Other GO categories were also enriched such as regulation of nuclear-transcribed mRNA catabolic process and positive regulation of nuclear-transcribed mRNA catabolic process (Fig. 3A). KEGG enrichment  analysis of differentially expressed genes can enrich pathways which are significantly different and is suitable for finding biological regulatory pathways with significant differences under experimental conditions. Count numbers higher than two and FDR-adjusted P value less than 0.05 were chosen as the cut-off criterion. Five hundred and ten KEGG pathways were identified in all the genes (data not shown). The top 30 KEGG pathways are shown in Fig. 3B. KEGG results showed that pathways in cancer were associated with 85 genes, which was the largest association. The interesting candidate pathways included the MAPK signaling pathway and Wnt signaling pathway, which has been reported to be involved in the pathogenesis of EM. Moreover, SNARE interactions in vesicular transport and mTOR signaling pathway were the top three pathways in terms of enriching factor. 3500 miRNA downstream target genes were predicted by IPA. DAVID was used to explore the KEGG pathway  Fig. 4A. We focused our attention on to mTOR signaling mechanisms because both TargetScan and IPA-predicted genes enriched in this pathway. Remarkably, four genes at different levels of the mTOR signaling pathway were targeted by miR-450a-5p and miR-503-5p ( Fig. 4B and C). Based on these results, we propose that the downregulated miRNAs may elicit the constitutive activation of the mTOR signaling pathway.

Discussion
We have seen numerous advances in the field of EM in the past few decades; however, only a few studies have focused on the pathogenesis, diagnosis and treatment of tubal EM. Till now, the gold standard of diagnosis of EM was laparoscopic visualization, ideally with histological verification (Zondervan et al. 2018). In EM, miRNAs play an emerging important role and many researches have explored the clinical applications of miRNA-based therapeutics (Panir et al. 2018). MicroRNA has also been studied as a potential diagnostic marker of EM (Nisenblat et al. 2016) since the presence of altered miRNA profiles in the plasma was reported in several types of solid tumors, such as ovarian cancer (Taylor & Gercel-Taylor 2008), and other diseases. Additionally, since a single miRNA targets multiple downstream genes, miRNA dysregulation can affect several pathways involved in the development of diseases. Unfortunately, the expression profile and the role of miRNAs in tubal EM is still unknown. We performed miRNA-microarray screening followed by qRT-PCR validation, and functional analysis to explore their roles in tubal EM. Bioinformatic analysis showed that specific miRNAs might promote the occurrence and development of tubal EM by taking part in SNARE interactions in vesicular transport, the mTOR signaling pathway and endocytosis. Tubal EM may lead to infertility, but the pathogenic mechanisms of it is still unclear. By regulating target genes and participating in mTOR signaling pathway, these differentially expressed miRNAs play potential role in the pathogenesis of it. Altered miRNAs and gene expression profiling affected micro-environment in fallopian tube and then may cause infertility, though the tubal epithelium showed no significant damage. Our study emphasizes miRNAs and related signaling pathways in the tubal epithelium of tubal EM which provides new clue to the pathogenesis of tubal EM.
In this study, we identified 17 differentially expressed miRNAs in the tubal epithelium from tubal EM. What we found was expected because miRNAs play a vital role in the regulation of diverse molecular pathways as a generic mechanism and as a core disease phenotype (Ben-Hamo & Efroni 2015), which are strongly linked with the pathogenesis of EM (Teague et al. 2010, Hawkins et al. 2011, Braza-Boils et al. 2014. Previous studies found that miRNA-21 was involved in tissue repair of EM and miRNA-223 was associated with angiogenesis of EM (Teague et al. 2010). miRNA-21 was overexpressed and was found to directly participate in the regulation and activation of a number of EM-related signaling pathways such as Wnt and TGF-β (Wang et al. 2012). miRNA-21 overexpression has also been found to contribute to cell proliferation by targeting PTEN in endometrioid endometrial cancer (Qin et al. 2012). It was found that the miR-503 expression in endometriotic cyst stromal cells was downregulated and that miR-503 inhibited cell proliferation, VEGF-A production and induced apoptosis and G0/G1 cell-cycle arrest of these cells (Hirakawa et al. 2016). A reduced expression of miR-503 has also been identified in endometrial cancer (Xu et al. 2013). These findings were consistent with our results that miRNA-21-3p and miRNA-223-3p were upregulated, whereas miRNA-503-3p was downregulated in tubal EM compared to the controls. These suggest that these differentially expressed miRNAs may play a potential role in the pathogenesis and development of tubal EM.
To explore how these miRNAs participate in tubal EM pathogenesis, an miRNA-target gene-disease network was generated, which identified the potential connections between miRNAs and their downstream genes. IPA identified 78 target genes (e.g. ALDOB and ZNF695) which were associated with endometrioid carcinoma, 79 (e.g. ARSK and ELL2) with female genital tract cancer, 79 (e.g. ARSK, CDK5RAP1 and HS6ST2) with tumorigenesis of reproductive tract, 86 (e.g. ARSK, ELL2 and TOX) with female genital tract adenocarcinoma, and 87 (e.g. ARSK, CERKL and ZNF484) with endometrial carcinoma. Most target genes were cancer-related. Previous studies showed that HS6ST2 was involved in ovarian cancer (Cole et al. 2014) and pancreatic cancer (Song et al. 2011). CDK5RAP1 was associated with tumor growth in breast cancer (Wang et al. 2015). Though tubal EM is not a malignant disease, some biological behaviors of it are similar to cancer, including invasion, adhesion, angiogenesis and proliferation. Therefore, we conclude that these target genes are associated with the clinical phenotype of tubal EM and the 17 differentially expressed miRNAs may regulate the pathogenesis of tubal EM by these target genes.
In this study, bioinformatic analysis and prediction of the downstream genes of miRNAs were conducted on miRNAs isolated from the fallopian tube epithelium from tubal EM patients and control patients. The five identified miRNAs were predicted to participate in a mTOR signaling pathway which is related to autophagy and cancer. Some studies have reported that in women with EM, this pathway is overexpressed and may significantly modulate survival, the proliferation of endometriotic cells and angiogenesis in endometriotic implants (Lee & Kim 2014). The mTOR pathway is associated with cell viability and proliferation and many kinases within this pathway may be overactive in endometriotic cells (McKinnon et al. 2016). Besides, it has been suggested that mTOR inhibitors can be used in the treatment of EM. Our results suggest that the downregulated miRNAs may participate in the pathogenesis of tubal EM by activating mTOR signaling pathway.
However, certain limitations in our study should be considered when interpreting the results -first, the number of samples was small, so, recruiting more participants is important for future studies; second, target genes related to biological functions, MFs and signaling pathways were obtained from existing databases. These findings must be further confirmed by mechanistic studies performed in tubal EM cell lines and animal experiments; third, the relationship between differentially expressed miRNAs and tubal reproductive functions needs further confirmation; fourth, patients with uterine leiomyoma were recruited as control group, which affects the accuracy of our results; therefore, we will select appropriate controls in further studies.
Overall, based on the microarray data, we found for the first time, an miRNA expression profiling and an associated miRNA-gene-disease network in tubal EM. We revealed that 17 miRNAs and 4343 genes were differentially expressed in the tubal EM patients compared to the control group. Most of the genes targeted by the five identified miRNAs appear to be cancer-related and anticipate primarily in the mTOR signaling pathway. Further studies on the pathogenesis of tubal EM are necessary to confirm the five identified miRNAs as pathogenic factors of tubal EM. Our results suggest that it is feasible to recruit more patients with tubal EM to study miRNAs and their target genes within tissues in order to provide novel idea to the pathogenesis of tubal EM-related infertility and to identify molecular targets for the diagnosis and treatment of tubal EM.

Supplementary data
This is linked to the online version of the paper at https://doi.org/10.1530/REP-18-0631.

Declaration of interest
The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Funding
This work was supported by Shanghai Jiaotong University (grant number YG2016ZD07).