Exploration of the pathological mechanism of endometriosis in rats using high-throughput sequencing

Background : This study aimed to explore the pathological mechanism of endometriosis (EMs) in rats by high-throughput sequencing. Methods : Rat EMs model, sham operation group and drug treatment group were established. Uterine tissue was collected for construction of sequencing library and high-throughput sequencing. Data quality was examined. KEGG pathway enrichment analysis was carried out. Results : Percentages of both high-quality sequence Reads and high-quality sequence bases accounted for more than 98%, suggesting that the data quality was acceptable. Total sequences mapped to reference genome (Total Mapped) accounted for more than 90% of total sequences used for mapping (Clean Reads), suggesting that the mapping results of the sequencing data were acceptable. There were 440 differentially expressed genes (DEGs) in the drug treatment group compared with the sham operation group, 382 DEGs in the drug treatment group compared with the model group, and 503 DEGs in the sham operation group compared with the model group. We screened genes ENSRNOG00000023079 and ENSRNOG00000012175 related to vascular endothelial growth factor pathway. The DEGs were mainly enriched in the signaling pathways such as phagosome, natural killer cell mediated cytotoxicity, Janus kinase-signal transducers and activators of transcription signaling pathway, hematopoietic cell lineage, cytokine-cytokine receptor interaction, regulation of actin cytoskeleton and extracellular matrix-receptor interaction. Conclusions : EMs might begin with the inflammatory response of the ectopic endometrium. key this process. The ectopic adhered to the abdominal wall with the help of the inflammation reaction, generated blood vessels, and finally transformed into growing lesions.


Abstract
Background : This study aimed to explore the pathological mechanism of endometriosis (EMs) in rats by high-throughput sequencing. Methods : Rat EMs model, sham operation group and drug treatment group were established. Uterine tissue was collected for construction of sequencing library and high-throughput sequencing. Data quality was examined. KEGG pathway enrichment analysis was carried out. Results : Percentages of both high-quality sequence Reads and high-quality sequence bases accounted for more than 98%, suggesting that the data quality was acceptable. Total sequences mapped to reference genome (Total Mapped) accounted for more than 90% of total sequences used for mapping (Clean Reads), suggesting that the mapping results of the sequencing data were acceptable. There were 440 differentially expressed genes (DEGs) in the drug treatment group compared with the sham operation group, 382 DEGs in the drug treatment group compared with the model group, and 503 DEGs in the sham operation group compared with the model group. We screened genes ENSRNOG00000023079 and ENSRNOG00000012175 related to vascular endothelial growth factor pathway. The DEGs were mainly enriched in the signaling pathways such as phagosome, natural killer cell mediated cytotoxicity, Janus kinase-signal transducers and activators of transcription signaling pathway, hematopoietic cell lineage, cytokine-cytokine receptor interaction, regulation of actin cytoskeleton and extracellular matrix-receptor interaction. Conclusions : EMs might begin with the inflammatory response of the ectopic endometrium.
Phagocytes played a key role in this process. The ectopic endometrium adhered to the abdominal wall with the help of the inflammation reaction, generated blood vessels, and finally transformed into growing lesions. Background 3 Endometriosis (EMs) is a disease caused by the presence of endometrial tissue with growth ability and function in other parts of the body besides the mucous membrane of uterine cavity. The name of EMs was first described by Von Rokitansky in 1860, but its etiology and pathogenesis have not yet been fully elucidated [1,2]. EMs is characterized by secondary and progressive exacerbation of dysmenorrhea. Distending pain of the lower abdomen and anus is the predominant symptom. It can begin 1-2 days before menstruation and disappear after menstruation. The pain begins at the lumbosacral and abdominal regions and radiates to the vagina, perineum, anus or thigh. The degree of pain is related to the location of the ectopic focus [3]. Epidemiological studies show that the incidence of EMs is on the rise and it is a common and frequently-occurring gynecological disease [4]. The incidence of EMs among women of childbearing age is more than 15%, which seriously affects women's physical and mental health, work and reproduction [4].
80% of the EMs patients have dysmenorrhea and 50% are complicated with sterility [1].
The traditional EMs treatment is to use drugs such as gestrinone to inhibit estrogen [5].
Most of these drugs are hormone compounds, which are not suitable for long-term use [5].
Studies have shown that the development of endometrial debris to EMs must undergo three processes: attachment, aggregation and angiogenesis [6]. With the deepening of research, it has been recognized that a variety of signaling pathways are involved, such as extracellular signal-regulated kinase 1/2 (ERK1/2)-mitogen-activated protein kinase (MAPK) pathway to promote the proliferation and adhesion of ectopic endometrial cells, wnt/β-catenin pathway to enhance the invasion and adhesion of ectopic endometrium, nuclear factor-κB (NF-κB) and transforming growth factor-β (TGF-β) pathways to promote angiogenesis and tissue remodeling [7][8][9][10]. Researchers are committed to investigating how to suppress the development and progression of EMs. Endostatin and blockers of vascular endothelial growth factor (VEGF) pathway can treat EMs by inhibiting 4 angiogenesis [11][12][13].
Therefore, in this study, we systematically analyzed the pathological mechanism of EMs using high-throughput sequencing in rat EMs model, which might help to reveal the mechanism of the development and progression of EMs.

Establishment of rat EMs model
The rats were anesthetized by intraperitoneally injecting 1% pentobarbital sodium at 45mg/kg. The cornua uteri on the side of the rats themselves was taken by operation and cut into small pieces. The endometrium was facing the abdominal cavity. The uterine tissue was sewn on the inner wall of the abdominal cavity. Then the abdominal cavity was closed by suture. Six rats were modeled by this method. Another three rats were taken as sham operation group. The rats in the sham operation group underwent laparotomy and abdominal suture only. A week later, three EMs rats were randomly selected as drug treatment group. Gestrinone was intragastrically administered on a daily dose of 60mg/kg once a day for 7 consecutive days. The remaining rats including another three model rats (the model group) and three rats in the sham operation group were intragastrically administrated with normal saline once a day for 7 consecutive days. Then the abdominal cavity of the rat was reopened and the uterine tissue on the inner wall of the abdominal cavity was taken. In the sham operation group, the corresponding abdominal wall was taken. The obtained tissues were used for subsequent extraction of RNA. Rats were anesthetized by intraperitoneal injection of 10% chloral hydrate at 3 ml/kg and then sacrificed by injection of 2 ml air through caudal vein.

Construction of sequencing library
Total RNA from the obtained tissues was extracted using TRIzol reagent (Invitrogen,

Data analysis
Data quality was examined using SepPrep 12.1 software (www.symantec.com) and Sickle 0.994 software (bioinformatics.ucdavis.edu). After quality control, the data were mapped to reference human genomic data using Bowtie2 [14]. KEGG pathway enrichment analysis was carried out by KEGG orthology based annotation system 3.0 (KOBAS 3.0) and KEGG annotation was performed [15][16].

Establishment of rat EMs model
The process of establishing the rat EMs model was shown in Fig. 1. Fig. 1a showed the ligation and taking of the cornua uteri. Fig. 1b showed the obtained endometrium. Fig. 1c showed that the endometrium was sutured to the abdominal cavity. Fig. 1d showed the wound closure in rats after operation.

Sequencing data reduction
Sequencing raw data contained some connectors and low-quality Reads which would greatly interfere with subsequent information analysis, so sequencing data needed to be sorted out. The criteria for data filtering were as follows: 1) Cutadapt was used to remove the 3' end connector, and the removed part had at least 10 bp overlap (AGATCGGAAG) with known connector, allowing 20% base mismatch; 2) Removal of the Read whose base recognition accuracy was less than 99%. The results of the filtered data were shown in Table 1. Both the percentage of high-quality sequence Reads to sequencing Reads and the percentage of high-quality sequence bases to sequencing bases in each group were more than 98%, which indicated that the data quality was acceptable.

Sequencing data mapping
The reference genome index was established through Bowtie2, and then the filtered Reads were mapped to the reference genome using Tophat2. When mapping using Tophat2, the mismatch between Reads and the reference genome sequence was not more than two, indicating the successful mapping. The mapping results between the sequencing data and the reference genome were shown in Table 2. Total sequences mapped to the reference genome (Total Mapped) accounted for more than 90% of total sequences used for mapping The three biological repetitions (rep1, rep2 and rep3) of the model group had little difference. In the drug treatment group, little difference was also found between the biological repetition 1 (rep1) and 3 (rep3) and the PC2 value of biological repetition 2 (rep2) was markedly higher than those of biological repetition 1 (rep1) and 3 (rep3). The difference between biological repetition 1 (rep1) and 2 (rep2) of the sham operation group was extremely small and their PC2 values were greater than those of the biological repetition 1 (rep1) and 3 (rep3) of the drug treatment group. In the sham operation group, the PC2 value of biological repetition 3 (rep3) was sharply higher than those of biological repetition 1 (rep1) and 2 (rep2), and was also obviously higher than that of biological repetition 2 (rep2) in the drug treatment group. Comprehensively, the gene expression among the three groups formed clusters respectively, which indicated that there were significant differences in the gene expression among the three groups. Table 3 summarized the differentially expressed genes (DEGs) obtained from the pair-wise comparison of the three groups. The differential expression of genes among the pair-wise comparisons of the three groups was shown in Fig. 3. There were 440 DEGs in the drug treatment group compared with the sham operation group. There were 382 DEGs in the 8 drug treatment group compared with the model group. There were 503 DEGs in the sham operation group compared with the model group. In addition, in the intersection of "the drug treatment group versus the model group" and "the sham operation group versus the model group", we screened the genes related to VEGF signaling pathway from the 86 DEGs were after excluding 4 DEGs in the intersection of the three pair-wise comparisons and found the genes ENSRNOG00000023079 and ENSRNOG00000012175. The corresponding false discovery rates (FDR) were 0.808848037 and 0.977323743, respectively.

Discussion
In this study, the rat EMs model and its corresponding sham operation group were 9 constructed. Three EMs model rats were intragastrically administrated with gestrinone to form the drug treatment group. The sham operation group served as control. The model group served as the up-regulation group of EMs related genes, while the drug treatment group served as the down-regulation group of EMs related genes compared with the model group.
Through high throughput sequencing, we found 440 DEGs in the drug treatment group compared with the sham operation group, 382 DEGs in the drug treatment group compared with the model group, and 503 DEGs in the sham operation group compared with the model group. Furthermore, we performed enrichment analysis of KEGG signaling pathway for DEGs obtained from the pair-wise comparison of the three groups and found that the DEGs were mainly enriched in the signaling pathways such as phagosome, natural killer cell mediated cytotoxicity, JAK-STAT signaling pathway, hematopoietic cell lineage, cytokine-cytokine receptor interaction, regulation of actin cytoskeleton and ECM-receptor interaction. It was suggested that EMs had a strong relationship with immune response and inflammatory response. A previous study had found that non-specific immunity, especially the change of macrophage phagocytosis, was one of the basic factors [17].
Immune factors secreted by macrophages such as prostaglandin E 2 (PGE 2) and interleukin 6 (IL-6) interweaved with VEGF to form a complex signal regulatory network and participated jointly in the angiogenesis of EMs [17]. The cytokine TGF-β is regulated by the JAK-STAT signaling pathway. STAT3 is highly active in the lesions of EMs patients, and TGF-β/Smad signaling pathway plays a key role in tissue fibrosis and has the effect of promoting ectopic endometrial adhesion [18][19][20]. Furthermore, we screened the genes related to VEGF signaling pathway and found the genes ENSRNOG00000023079 and ENSRNOG00000012175. However, the VEGF signaling pathway was not prominent in the enrichment of KEGG signaling pathway. The corresponding FDR of these two genes was also high, probably because the model itself was in the early stage of EMs and the angiogenesis process had not reached its peak yet.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. sham: the sham operation group   Figure 1 The process of establishing the rat endometriosis model. a: the ligation and taking of the cornua uteri. b: the obtained endometrium. c: the endometrium was sutured to the abdominal cavity. d: the wound closure in rats after operation