Differential expression and functional analysis of micro RNAs in Papio anubis induced with endometriosis for early detection of the disease

1.   Introduction

Endometriosis is a common infertility-related gynecologic disorder defined by the presence of endometrial-like tissue outside the uterine cavity. The disease is estimated to affect up to 10% of women in their reproductive years worldwide [1]–[3] and infertility is a common outcome in 30–50% of these patients [4]. Women with endometriosis frequently suffer from symptoms including non-menstrual pelvic pain, painful menstrual cramps, and pain during intercourse, fatigue, and infertility [1],[5]. This negatively affects their quality of life [6]. Due to its chronic and asymptomatic character, lack of treatment, high prevalence, and significant morbidity associated with the disease, the global endometriosis-related healthcare burden has been estimated to be in the hundreds of billions of Euros each year [7].

To date, non-invasive approaches such as ultrasound and magnetic resonance imaging have not been effective in the diagnosis of endometriosis [8]. The current gold standard for the diagnosis of endometriosis is laparoscopic examination with histological confirmation of glands and/or stroma in the excised lesions [9]. The invasive nature and cost of surgery, misinterpretation of symptoms [10] coupled with the lack of molecular biomarker(s), results in a diagnostic delay of between 4–10 years from the onset of symptoms to definitive diagnosis [9]. The prolonged delay contributes to years of suffering, potential infertility, and disease progression in up to 50% of the affected women [11]. Non-invasive diagnosis of endometriosis would allow early diagnosis and treatment but so far this has not been achieved meaning it is a priority in endometriosis research [4].

Owing to ethical reasons and differences in progression of the disease, controlled invasive studies cannot be carried out in humans. An appropriate animal model is of paramount importance in the search for diagnostic biomarkers. A model of experimentally induced endometriosis in baboons with documented regular menstrual cycles has been developed by several groups [12]–[15]. Intraperitoneal inoculation with autologous menstrual endometrium results in the formation of endometriotic lesions in the baboons with histological and morphological characteristics similar to those in women [13]–[15]. Obtaining endometrial biopsy (EB) samples by either pipelle or curette is minimally invasive and has previously proven to be useful as a diagnostic tool for endometriosis in an outpatient setting [16],[17]. According to the widely accepted theory of retrograde menstruation [18], the menstrual endometrium is the source of ectopic endometriotic foci and therefore using the direct source of the disease is perhaps logical in our quest to identify biomarkers for endometriosis. Further, the finding that retrograde menstruation is present in 90% of women but not all of them suffer from endometriosis [1] suggests that molecular differences between eutopic endometrium from women with and without endometriosis may exist that lead to the development of the condition in certain women but not in others [1]. The endometrium of women with endometriosis has been documented to respond differently to ovarian hormones [19] and therefore, this offers the opportunity to discover new potential biomarkers in tissue biopsies obtained through a minimally invasive procedure [20].

MicroRNAs (miRNA) are single-stranded, noncoding, small RNA molecules that regulate gene expression by inhibiting mRNA translation or by facilitating cleavage of the target mRNA [21],[22]. MicroRNAs are attractive as biomarkers because of their lower complexity, no known post-processing modifications, simple detection and amplification methods, tissue-restricted expression profiles, and sequence conservation between humans and model organisms [23]. In addition, miRNAs are robustly expressed and conserved [24]. MicroRNAs are likely to regulate mRNAs and molecular networks that contribute to the key pathways proposed in the pathophysiology of endometriosis, like inflammation, angiogenesis, tissue repair, and extracellular matrix [25],[26]. Women with endometriosis have distinct circulating miRNA signatures, as identified in several studies [25],[27],[28]. A number of studies have reported differential expression of miRNAs between eutopic and ectopic endometrial tissues [24],[25],[29]–[33]. Previous studies have shown dysregulation of miRNAs in endometrium and lesions of endometriosis patients [20],[34],[35]. However, different studies have not reached a consensus on which particular miRNAs are most relevant in endometriosis, as dysregulated miRNAs reported in one study are only occasionally confirmed by others.

In this study, we analyzed miRNA profiles in eutopic endometrial tissues from baboons induced with endometriosis. This provided a more controlled assessment of the disease. We performed a genome-wide miRNA expression analysis in endometrium biopsy (EB) samples to identify the differentially expressed miRNAs. Additionally, their targets were predicted and their functions in relation to the endometriosis were determined.

2.   Materials and method

2.1.   Ethical statement

The study was approved by the Institutional Scientific Evaluation and Review Committee and the Animal Care and Use Committee of the Institute of Primate Research (IPR) Nairobi – Kenya.

2.2.   Study design

A case-cohort study involving Papio anubis was undertaken at Institute of Primate Research (IPR)-Nairobi. Three adult female olive baboons with an average mean weight of 15.2 kg were used in the study. The baboons were captured in the wild and maintained in quarantine for 3 months at IPR animal cages and were screened for common pathogens (bacterial and viral infections as well as parasites) and tuberculosis, simian T-lymphotropic virus-1, and simian immunodeficiency virus to ensure they were disease-free. They were housed indoors with natural lighting in group cages and fed on commercial food pellets with fruits and vegetable supplementation three times a week and water ad libitum.

2.2.1.   Induction of endometriosis

To confirm absence of endometriosis, the baboons were screened by video laparoscopy during the mid-luteal phase (approximately 25^th day of the cycle) and animals were then allowed to recover for one menstrual cycle as previously described [36]. Disease induction was experimentally done in baboons (n = 3) as previously described [15]. Briefly, on days 1–2 after the onset of the next menses, 1 gram of menstrual endometrium was harvested by transcervical uterine curettage from each animal and minced through an 18 – gauge needle. The menstrual endometrium was autologously seeded onto ectopic sites (uterosacral ligaments, uterovesical fold, pouch of douglas, ovaries) as described before [37] . For the surgical procedure, the baboons were anesthetized with a mixture of ketamine (Anesketin, 15 mg/kg; Eurovet NV/SA, Heusden-Zolder, Belgium) and xylazine (Bomazine 2%, 2 mg/kg; Bomac Laboratories Ltd, Auckland, New Zealand) administered intramuscularly for induction, and 1wer%–2% halothane (Halothane; Nicholas Piramal India Ltd, Andhra Pradesh, India) with N2O/O2 (70%/30%) for maintenance. After surgery, the animals received antibiotics for 1 week (Clamoxyl LA; Pfizer, Paris, France), and pain was controlled with ibuprofen (Ketofen; Merial, Lyon, France) for 3 days.

2.2.2.   Sample collection and preparation

500 milligrams of endometrial tissue samples were collected using a suction curette (Pipet Curet, CooperSurgical, USA) in 1ml sterile tube and processed within 1 hour after collection. The samples were collected from day zero (before disease induction), day 25 and 50 Post Induction (PI) through video laparoscopy procedures performed by veterinary doctors. A total of ten samples from the baboons (n = 3) were used for the study; three control samples at baseline, three diseased on day 25 and three on day 50 PI while one sample was a replicate. (Each sample was divided into two portions: one of which was fixed in 10% formalin and processed for histological evaluation (hematoxylin-eosin [H-E]); the second portion was frozen at −80°C for subsequent RNA extraction. The endometrial tissue biopsies were homogenized in TRIzol™ Reagent (Invitrogen, USA) before RNA extraction.

2.3.   RNA extraction

100 µl of total RNA was extracted from EB using the Direct-zol RNA MiniPrep kit (Zymo Research - USA) according to the manufacturer's protocol. The purity and concentration of RNA was determined by OD260/280 from a NanoDrop spectrophotometer (ND-1000, Thermo Fisher Scientific-USA). RNA with an OD260/280 between 1.8 and 2.0 was considered of good-quality and was included in further experiments.

2.4.   Library construction and Small RNA (sRNA) Sequencing

The RNA samples were subjected to library construction using SMARTer smRNA-Seq Kit (Clontech Laboratories, USA) according to the manufacturer's protocol [38]. Briefly, 6µl of total RNA was used as the input for RNA adapter ligation (using 3′ and 5′ RNA adapters). This was followed by reverse transcription and PCR amplification with bar-coded primers. To verify the size of PCR enriched fragments, template size distribution was checked by running on an Agilent Technologies 2100 Bioanalyzer using a DNA 1000 chip [39]. Further, the prepared libraries were quantified using qPCR according to the Illumina qPCR Quantification Protocol Guide [40]. Sequencing of the libraries was achieved under Illumina SBS Technology platform utilizing a proprietary reversible terminator-based method to detect single bases as they are incorporated into DNA template strands [41]. The platform was selected to minimize natural competition and base incorporation bias as well as reduce raw error rates compared to other technologies. Raw sequencing reads from the output were retrieved from Macrogen servers and submitted to Sequence Reads Archive (SRA) Bioprojects PRJNA625977 and PRJNA637777.

2.5.   Data analysis

2.5.1.   Preprocessing, mapping and quality control

Sequences in FASTQ format were exposed to quality control check using FastQC (Babraham Bioinformatics; UK; version 0.10.1). The raw sequence reads were filtered and only those with a Phred score of >20 and length between 15–50 nt (to include sequences with sustained 3′-adaptor sequences at the 3′-end) were retained. The adaptor sequences were trimmed using Cutadapt [42] with the following parameters: set adapter sequence -AAAAAAAAAA-, default adapter detection at 10% error rate and length filter set at ≥15 nt. Sequences that failed the quality cutoff were regarded as unusable and were not considered for subsequent analysis.

2.5.2.   miRNA identification and distribution of sRNA

Considering that some of the sequence reads may not fall under miRNA category but still belong to sRNA family as well as the possibility of miRNAs mapping to non-miRNA positions in the genome, we adopted the two sRNAbench mapping strategies [43]. These strategies include genome mapping approach to target the Papio anubis (P. anubis ) genome and library mapping strategy targeting mature miRNA in mirBase database [44]. The processed reads from 10 samples (control samples n = 3, day 25PI n = 3, day 50PI n = 3, replicate n = 1) were allotted in the genome mapping strategy. In this strategy, the reads were first aligned to the bowtie-indexed reference transcript of Papio anubis (Panu_3_0) in the Ensembl database. The transcript mapping was adopted to cutoff artifacts that map to non-coding genomic regions of P. anubis. For the library strategy, Bowtie tool incorporated in sRNA workbench was used in mapping sequence reads against the miRNA library (hairpin.fa + mature.fa) of human (hsa) deposited in miRBase v. 22 [44]. This was done as the miRNA library annotated for P. anubis at the time of our analysis was not available. Homo sapiens miRNA library are highly homologous to Papio anubis miRNAs (∼90% homology) as suggested in previous studies [45],[46]. Reads that were mapped to the reference Ensembl were given the coordinates for their exact location in the chromosome.

2.5.3.   Differential expression (DE) of mature miRNAs and prediction of novel miRNAs

Due to lack of consistency in recalling some miRNAs across the samples of the same treatment group, we relied on 38 miRNAs that were consistently expressed across all samples (control and diseased samples). After excluding miRNAs with less than 10 raw reads from the libraries, differential expression analysis for the selected miRNA (38 miRNAs) was done. For genome-mapped reads, we used the top 300 coding regions to identify differentially mapped (hotspots) non-microRNA sites among the P. anubis transcripts. To determine and visualize differentially expressed mature miRNAs and genomic hotspots, DeSeq and Heatmap packages in R (version 3.6.6) were used respectively. False discovery rate (FDR) cutoff value for differentially expressed miRNAs and genes between control and diseased samples was set at FDR < 0.05 while DeSeq expression threshold was adjusted as p = 0.05. Prediction of novel miRNAs was done using prediction tool incorporated in sRNAbench following the default parameters as set in the prediction manual [47]. A venn diagram software was used to establish the relationship between the DE miRNAs on day 25 and 50 PI.

2.5.4.   Prediction of miRNA targets and functional analysis of DE miRNAs and genes

To investigate the biological role of the twelve dysregulated miRNAs, the potential target genes were predicted using TargetScan [48] and DIANA's mirPATH 3.0 software [49]. Enrichment analysis including pathway and functional annotation of DE miRNA target genes (MTGs) and dysregulated genes was done using shinyGo and DIANA-mirPath v.3.0 web server [49]. To investigate the potential functional mechanism of the dysregulated miRNAs, we employed enrichment Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the DIANA-mirPath v.3.0 web server [49] with a FDR cutoff of p = 0.05.

3.   Results

3.1.   Identification of miRNAs and sRNA distribution in normal and diseased eutopic endometrium of P. anubis.

Sequencing of 10 EB small RNA libraries from normal and diseased samples yielded a total of 103,261,611 sequence reads with an average of 10,326,161 reads per sample (raw sequences). Removal of adaptors, contaminants and low-quality reads resulted in a total of 32,025,676 clean reads from the libraries with an average of 3,202,568 clean reads per sample (Supplementary Figure S1A). Majority of the sequence reads that mapped on the genome were of length 21–24 nt (Supplementary Figure S1B). The small RNAs exhibited a diverse size distribution of sequence reads that aligned to P. anubis genome where majority were miRNAs (Figure 1A). Apart from sRNAs, other short RNA fragments in the samples were also detected during the analysis (Figure 1A). The total high-quality reads that perfectly mapped to P. anubis genome was 77.43% (Figure 1B).

Figure 1.  Identification of EB miRNAs and other RNAs from P. anubis by Illumina SBS sequencing. (A) shows percentage read composition of each read length. In (A), the order from bottom upwards include portions of miRNA, snRNA, tRNA, ncRNA, rRNA, mRNA and other RNAs. (B) shows the percentage of total number of high quality reads that were mapped in both genome and library (miRBase) mode.

DownLoad: Full-Size Img PowerPoint

3.2.   Differential Expression analysis of miRNAs and genes in normal and diseased EB from olive baboons

3.2.1.   DE analysis of miRNAs

Comparison of individual microRNA expression levels between normal eutopic endometrium (X0dpi1-X0dpi3) and eutopic endometrium of baboons with endometriosis (X25dpi1-X25dpi3 and X50dpi1-X50dpi4) was done. Due to lack of consistency in recalling some miRNAs across the samples of the same treatment group (baseline, day 25& 50PI), we relied on 38 miRNAs which were consistently expressed across all the samples (diseased and control groups). Differential expression analysis for the selected 38miRNAs in both the diseased and control EB samples was done as shown in Figure 2A. There were twelve miRNAs that were significantly differentially expressed in diseased eutopic endometrium compared to control endometrium (Table 1). Seven miRNAs were significantly up-regulated while five miRNAs were significantly down-regulated in the diseased baboons compared to controls (Table 1 and Figure 2B). The Venn diagram in Figure 2B depicts differentially expressed miRNAs through each time point of disease and overlap between these transitions. Of these, six miRNAs (miR-143-3p, miR-145-5p, miR-199a-3p, miR-199a-5p, miR-214-3p and miR-125b-5p) and four miRNAs (miR-145-5p, miR-199a-3p, miR-214-3p and miR-10b-5p) were over-expressed on day 25 and 50PI respectively (Table 1, Figure 2B). Five miRNAs (miR-16-5p, let-7g-5p, miR-29b-3p, miR-378a-3p and miR-342-3p) and two miRNAs (miR-378a-3p and let-7g-5p) were down regulated on day 25 and 50PI respectively (Table 1, Figure 2B).

Table 1.  Expression profiles of miRNAs between normal and diseased eutopic EB from Olive baboons at day 25&50 PI identified by DESeq program (FDR < 0.05) and P = 0.05.

EB (day 25PI)
Name	p-Value	Fold change	Status (Up/Down regulated)
miR-16-5p	1.83E-05	0.117211287	Down
let-7g-5p	1.65E-05	0.113225127	Down
miR-29b-3p	2.73E-09	0.075663826	Down
miR-378a-3p	3.49E-12	0.142046673	Down
miR-342-3p	0.003417	0.418169826	Down
miR-143-3p	0.00016	2.987779	Up
miR-145-5p	4.01E-11	15.92506	Up
miR-199a-3p	4.19E-08	6.629909	Up
miR-199a-5p	0.000121	3.113243	Up
miR-214-3p	1.00E-08	8.490738	Up
miR-125b-5p	5.46E-05	4.194129	Up

---

 | Show Table

DownLoad: CSV

Figure 2.  Most abundantly expressed miRNAs in diseased and normal eutopic EB. (A) Representative heatmap for the top 38 expressed miRNAs. The map was generated using log2 of median RPM. Dendrograms on the left represents Pearson correlations from hierarchical clustering analysis on the miRNAs. X0dpi-X0dpi3 miRNAs from control samples while X25dpi1-X25dpi3 & X50dpi1-X50dpi3 are miRNAs expressed on day 25&50 PI respectively. The red indicated high miRNA expression level and the blue shows low miRNA expression level. (B) Venn diagram showing the relationship between down and up regulated miRNAs at two time points after induction of endometriosis. The Venn diagram includes both up-regulated (UR) and down-regulated (DR) miRNAs on day 25 and 50PI (p ≤ 0.05). Three miRNAs were common (overlapped) in the two disease time points while 3 miRNAs were unique on day 25P1 and 1 was unique on day 50PI. An overlap of 2 DR miRNAs on day 25 &50 PI is shown while 3 miRNAs were unique on day 25PI while none was unique on day 50PI.

DownLoad: Full-Size Img PowerPoint

3.2.2.   Differentially expressed genes based on RNA abundance (Mapping against Ensembl transcript database)

For genome-mapped reads, we used the top 300 recalled coding domains to identify differentially expressed (hotspots) non-microRNA sites within the Papio anubis transcripts. From these 300 transcripts, the top 40 transcripts were presented in the heatmap due to space restriction (Figure 3). We identified seventeen genes that were significantly differentially expressed on the two time points (Table 2). Eleven genes were significantly down regulated while six genes were upregulated in diseased compared to the normal subjects (Table 2).

Table 2.  Differentially expressed genes derived from P. anubis Ensembl transcripts based on RNA abundance identified by DESeq program (FDR < 0.05) and P = 0.05.

Hotspots transcripts	Pvalue	Fold Change	Up/Down regulated
ITGA5-201	1.38E-02	7.00E-02	Down
LTB	0.002517	0.120625	Down
ADCY6-202	5.54E-02	2.54E-01	Down
ADCY6-201	5.54E-02	2.54E-01	Down
ADCY6-203	5.72E-02	2.56E-01	Down
PCDHGA4-205	5.66E-01	1.75E+00	Down
CMIP	0.000341783	0.143671631	Down
RNF150	0.002357946	0.121402398	Down
HDGFL3	0.000283676	0.223186038	Down
MDN1	1.05E-01	2.29E+00	Down
FOXC2	0.002725407	0.128973737	Down
KMT2C-201	4.27E-03	3.21E+00	Up
SZT2-201	9.28E-03	3.06E+00	Up
SZT2-202	9.28E-03	3.06E+00	Up
TBX4-201	6.08E-02	2.15E+00	Up
DCLK3	3.66E-06	38.60476048	Up
CYP3A5	0.001855972	19.14041971	Up

---

 | Show Table

DownLoad: CSV

Figure 3.  Representative heatmap for the top 40 hotspots across Papio anubis genome. The map was generated using log2 of median RPM. Dendrograms on the left and above the heatmap represents Pearson correlations from hierarchical clustering analysis on the genes and samples respectively. Blue colour indicates low expression while red colour indicates high expression.

DownLoad: Full-Size Img PowerPoint

3.3.   Identification of Novel miRNAs

Some novel miRNAs turned out to have very low or even no expression in our data and were discarded as candidates for novel miRNAs. In total, 20 potential novel miRNAs were predicted from sRNA sequencing data of EB using sRNAbench (Supplementary Table S1). Seven novel miRNAs showed similarities with miR451 (M-P001-M-P007), one was similar to miR342 (M-P009) while another novel miRNA was similar to miR310 (M-P008) according to the miRBase database (Supplementary Table S1). Eleven miRNAs had no similarity with any miRNA in the miRbase and therefore were deposited to enrich the miRBase database.

3.4.   MiRNA target gene identification and functional analysis

3.4.1.   Gene Ontology analysis

GO enrichment analysis showing potential biological roles of the MTGs using shinygo at a p-value = 0.05. The DE miRNAs target genes are represented in Supplementary Table S2. There were 30 GO biological processes (BP) terms significantly enriched for MTGs (Figure 4A) while for the DE genes, the most significant enriched terms were twenty at a p-value of < 0.05 (Table 3). Two DE genes which are ITGA5-201 and FOXC2 were involved in most of the GO biological processes. The Relationship between the target genes of these miRNAs was established using a Venn diagram software. There were 142 target genes that were common in both the upregulated and down regulated MTGs while 1,374 genes were uniquely identified as targets for the down regulated miRNAs as shown in Figure 4B. For the up-regulated miRNAs, 5,826 genes were uniquely identified as their targets (Figure 4B).

3.4.2.   KEGG pathway enrichment analysis

Significant KEGG enrichment analysis was done using Fishers Exact Methodology Test incorporated in DIANA-mirPath software at a p-Value of 0.05. A detailed summary of the DE miRNAs, the KEGG pathways and the number of targeted genes is provided in Table 4. There were 67 pathways that were significantly enriched but we selected sixteen pathways based on their association with endometriosis. The sixteen pathways include endometrial cancer, adherens junction, progesterone-mediated oocyte maturation, TNF signaling pathway, ECM-receptor interaction, estrogen signaling pathway, signaling pathways regulating pluripotency of stem cells, focal adhesion, FoxO signaling pathway, proteoglycans in cancer, p53 signaling pathway, steroid biosynthesis, transcriptional mis-regulation in cancer, PI3K-Akt signaling pathway, TGF-beta signaling pathway and pathways in cancer. A heat map was done to display the relationship between the biological pathways and the DE miRNAs (Figure 5).

Table 3.  GO analysis for DE genes.

Enrichment FDR	Genes in list	Total genes	Functional Category	Genes
0.022142	2	26	Positive regulation of sprouting angiogenesis	ITGA5-201 FOXC2
0.028686	2	49	Regulation of sprouting angiogenesis	ITGA5-201 FOXC2
0.028686	2	51	Cell adhesion mediated by integrin	ITGA5-201 FOXC2
0.030893	3	338	Angiogenesis	ITGA5-201 FOXC2 TBX4-201
0.030893	2	97	Formation of primary germ layer	ITGA5-201 FOXC2
0.030893	3	379	Morphogenesis of an epithelium	ITGA5-201 FOXC2 TBX4-201
0.030893	2	83	Sprouting angiogenesis	ITGA5-201 FOXC2
0.030893	2	79	Regulation of tube size	FOXC2 ADCY6-201
0.030893	2	78	Regulation of blood vessel size	FOXC2 ADCY6-201
0.030893	2	93	Vascular process in circulatory system	FOXC2 ADCY6-201
0.035742	3	417	Blood vessel morphogenesis	ITGA5-201 FOXC2 TBX4-201
0.035898	2	125	Positive regulation of vasculature development	ITGA5-201 FOXC2
0.035898	3	487	Blood vessel development	ITGA5-201 FOXC2 TBX4-201
0.035898	2	139	Gastrulation	ITGA5-201 FOXC2
0.035898	3	508	Vasculature development	ITGA5-201 FOXC2 TBX4-201
0.035898	2	115	Positive regulation of angiogenesis	ITGA5-201 FOXC2
0.035898	3	470	Embryonic morphogenesis	ITGA5-201 FOXC2 TBX4-201
0.035898	3	466	Tissue morphogenesis	ITGA5-201 FOXC2 TBX4-201
0.035898	3	516	Cardiovascular system development	ITGA5-201 FOXC2 TBX4-201
0.035898	2	146	Cell-cell adhesion via plasma-membrane adhesion molecules	ITGA5-201 PCDHGA4-205

---

 | Show Table

DownLoad: CSV

Figure 4.  (A) A hierarchical clustering tree summarizing the correlation among significant pathways associated with MTGs. Pathways with many shared genes are clustered together. Bigger dots indicate more significant P-values. (B) Comparison between up and down-regulated miRNA target genes. An overlap of genes between the up-regulated and down regulated target genes involved 142 genes while 1,374 genes were uniquely identified as targets for down-regulated miRNAs. The target genes uniquely identified for the up-regulated were 5,826.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Functional analysis of dysregulated miRNAs based on Gene ontology biological process pathways. Clusters on the right represent the dysregulated miRNAs and the prediction tool used in identifying the targets. GO biological processes are clustered under the heatmap. This was done using DIANA's mirPATH 3.0 tool. Red colour shows the most significant BP involving each DE miRNA.

DownLoad: Full-Size Img PowerPoint

Table 4.  Enriched KEGG pathways involving the dysregulated miRNAs. Analysis done by DIANA miRPath v.3.0 software (p-value = 0.05).

KEGG pathway	p-value	miRNAs
Endometrial cancer	0.0010266851934	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Adherens junction	7.596765016e-11	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Progesterone-mediated oocyte maturation	0.0312943265702	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
TNF signaling pathway	0.0034735818690	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
ECM-receptor interaction	1.044183390e-10	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Estrogen signaling pathway	0.0025192410437	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Signaling pathways regulating pluripotency of stem cells	5.03291855e-06	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Focal adhesion	1.1704258648e-05	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
FoxO signaling pathway	0.0015352923167	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Proteoglycans in cancer	2.016100718e-21	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
p53 signaling pathway	0.0001167150178	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Steroid biosynthesis	0.0098103022005	let-7g-5p, miR-378a-3p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Transcriptional mis-regulation in cancer	0.0010266851934	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
PI3K-Akt signaling pathway	0.007019728167	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
MAPK signaling pathway	0.033769490039	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
TGF-beta signaling pathway	2.485180553e-08	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Pathways in cancer	1.150577629e-07	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p
Cell cycle	1.044183391e-10	let-7g-5p, miR-378a-3p, miR-125b-5p, miR-143-3p, miR-10b-5p, miR-16-5p, miR-214-3p,miR-145-5p, miR-29b-3p,miR-199a-3p,miR-199a-5p, miR-342-3p

---

 | Show Table

DownLoad: CSV

4.   Discussion

Differential expression analysis of miRNA between endometriotic lesions and eutopic endometrium from women with endometriosis has been previously reported [50] but few studies have focused on differences between eutopic endometrium from women with and without endometriosis [27],[50],[51]. In the present study, we evaluated the expression profiles of miRNAs in eutopic endometrium as potential biomarkers for early diagnosis of endometriosis. We experimentally induced endometriosis in olive baboons to provide a more controlled assessment for the disease since in most human studies, a proper control for the disease has been a limitation. Our findings reveal twelve miRNAs that were significantly dysregulated in eutopic endometrium of olive baboons induced with endometriosis compared with normal endometrium during early disease stage. Seven miRNAs namely miR-199a-3p, miR-145-5p, miR-214-3p, miR-143-3p, miR-125b-5p, miR-199a-5p and miR-10b-5p were significantly upregulated while five miRNAs which are miR-29b-3p, miR-16-5p, miR-342-3p, miR-378a-3p, let-7g-5p were significantly down regulated. The dysregulated miRNAs from our findings could be attributed to their roles in development and progression of endometriosis.

The up-regulation of miR-143 and miR-145 is in agreement with other studies that have shown similar results in patients with endometriosis [23],[25],[52]. These MiRNAs (miR-145 and miR-143) are important in proliferation, inflammation, apoptosis, invasion, growth and differentiation as previously reported [53]. Elevated levels of miR214 have been reported in previous studies on reproductive disorders [54],[55] and has been shown to induce cell survival and promote metastases by targeting PTEN gene which plays an important role in the development of endometriosis [56]. Our findings on increased expression levels of miR199a-3p and miR199a-5p have been replicated before [57] and their roles in pelvic adhesion, inflammation, angiogenesis, cell proliferation and lesion distribution during endometriosis have been reported. Similar findings on over expression of miR-125b-5p and its role in cell proliferation, angiogenesis and adhesion have previously been reported [25]. Our results on up regulation of miR-10b-5p and its apoptotic function was earlier reported [58],[59].

The decreased expression levels of miR-16-5p, miR-29b-3p, miR-378a-3p, miR-342-3p and let 7g-5p from our results have been replicated in previous studies [51],[60]–[63]. The roles of miR-16-5p in endometriosis development include apoptosis, cell proliferation, inflammation and angiogenesis [64] while miRNA 29b-3p is involved in cell differentiation. Angiogenesis, cell migration and invasion processes have been associated with miR-378a-3p. Micro RNA 342-3p has been linked to alteration of lipid metabolism leading to altered immunity during endometriosis while let 7g-5p plays a role in cell proliferation and survival.

Our findings from GO enrichment analysis show that the dysregulated miRNAs regulate several processes involved in the pathogenesis of endometriosis such as inflammatory and immune response, cell invasion, extracellular matrix remodeling, angiogenesis, cell proliferation, cell communication, signaling, epigenetic regulation, metastasis and apoptosis. These processes which are important in pathogenesis of endometriosis have been previously reported in other studies [26],[65]. The KEGG analysis of the dysregulated miRNAs from our study reveal important pathways that have been previously linked with endometriosis. We selected the pathways that had an association with endometriosis as shown in Table 4. These selected pathways regulate cell growth, adhesion, proliferation, differentiation, angiogenesis and apoptosis which are crucial in development and progression of endometriosis [66]. For instance, in patients with endometriosis, progesterone and eostrogen pathways are altered and therefore affecting several pathways such as progesterone-mediated oocyte maturation, estrogen signaling pathway, PI3K-AKT signaling pathway, MAPK signaling pathway and FOXO signaling pathway shown on Table 4 which have been previously supported [67].

We also reported seventeen differentially expressed genes based on the abundance of other small RNAs sequences in our EB samples. The functional analysis reveal that some of these DE genes are putative target genes for the DE miRNAs. The enriched results of GO analysis and KEGG pathways suggested that these dysregulated genes were significantly enriched in angiogenesis, cell and biological adhesion and oxidative stress which are pathological endometriosis processes. These genes have been previously associated with endometriosis [68]–[70]. For instance, Integrin alpha 5 (ITGA5) identified in our study is important in cell adhesion which is a process in establishment and progression of endometriosis. We demonstrate that the diseased eutopic endometrium compared to normal (control) endometrium exhibits molecular abnormalities. Our findings have also uncovered novel miRNAs never described before and therefore will be used to enrich the miRbase.

From our study, there was a high correlation between the identified DE miRNAs in the baboon model and altered gene expression in endometrium of women with endometriosis as we have discussed above. The selected time points for the disease induction (day 25 and 50PI) in our study represent early disease stages as this has previously been established in other studies [12],[37],[71]. Staging of endometriosis in the baboons was done using revised American Fertility Society (rAFS) scoring system [72] after modification to baboon size [37]. We collected samples prior to the induction of disease and then compared these control samples to diseased samples after disease induction. The use of each animal to serve as its own control reduces inter-animal variability in experimental designs and reduces the animal numbers required per study [15]. Although there were miRNAs that were dysregulated in both time points of the disease (day 25 and 50 PI), some were dysregulated only in a single disease time point. More studies with a large number of baboons to confirm this are needed. Aberrant gene expression and molecular abnormalities in eutopic endometrium during early disease stages can accurately be determined in baboons as opposed to humans [71]. The key pathways identified in our study in relation to endometriosis are important in understanding the mechanisms involved in the development of disease.

5.   Limitations

Our study had several limitations which include small sample size, inclusion of only minimal and mild stages of disease, and collection of control tissues only at a single time point before disease induction (baseline) and not at consecutive time points for exact matching. Further studies to validate these DE miRNAs in eutopic EB with large number of baboons with induced endometriosis are necessary. Despite these limitations, we identified DE miRNAs and genes which may be considered for further evaluation as potential biomarkers for minimally invasive diagnosis of endometriosis in patients.

6.   Conclusion

Our findings indicate that altered expression of specific miRNAs is associated with endometriosis and therefore further evaluation of these miRNAs as potential biomarkers for the disease diagnosis could lead to production of a diagnostic kit. The differential expression of these miRNAs and their putative molecular pathways constituted by their target genes reveal that these miRNAs could be involved in development and progression of endometriosis. Our KEGG analysis results reveal important pathways targeted by these dysregulated miRNAs in relation to endometriosis which could be further explored as therapeutic targets.